PROGRAM wave2d_30x60 c PROGRAM TESTZ C C23456789012345678901234567890123456789012345678901234567890123456789012 C 1 2 3 4 5 6 7 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * M A I N P R O G R A M * C * NAME: ( Q4FEWA ) * C * Q UADRILATERAL 4-NODED F INITE E LEMENT W AVE A NISOTROPIC * C * * C * MODIFIED TESTG FOR GENERATION AND TRACKING OF ONE PULSE. * C * BOUNDARY SIMULATING PULSE REMAINS FIXED DISPLACEMENT. * C * * C * CORE STORAGE MINIMIZED * C * 1. STORES ELEMENT-NODE DATA ON SCRATCH FILES (TAPE 1,3,7-9). * C * 2. DIMENSIONED VARIABLES SHARE HIGH SPEED CORE STORAGE. * C * 3. USES THE SKYLINE STORAGE SCHEME, DEVELOPED BY BATHE AND * C * WILSON ( NUMERICAL METHODS IN FINITE ELEMENT ANALYSIS), TO * C * STORE THE STIFFNESS ARRAY AS A SINGLE DIMENSIONED ARRAY * C * SUCH THAT A COLUMN ORIENTED GAUSS ELIMINATION TECHNIQUE IS * C * USED IN SUBROUTINE COLSOL TO SOLVE THE SET OF LINEAR EQNS. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C REAL NU12,NU13,NU23 character tt*8 C COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E1,E2,E3,G12,G13,G23,NU12,NU13,NU23,AL1,AL2,AL3,THETA COMMON/RDK3/A0,A2,A3,A6,A7,WC,RHO,AMP C DIMENSION TITLE(20),DD(3,3) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * THE FOLLOWING TWO CARDS DETERMINE THE MAXIMUM HIGH SPEED CORE * C * STORAGE SHARED BY DIMENSIONED VARIABLES IN VARIOUS SUBROUTINES. * C * TO CHANGE STORAGE, CHANGE BOTH MTOT AND SHARE(MTOT). * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C c COMMON SHARE(3222231) c MTOT=3222231 C COMMON SHARE(491557) MTOT=491557 C C ASSIGN TAPE NO.S 5 AND 6 TO INPUT AND OUTPUT FILES RESPECTIVELY C LR=5 LW=6 open(lr,file='wave2d_30x60.DAT',status='old',err=8888) open(lw,file='wave2d_30x60.out',status='unknown',err=8888) c c Assign tape no.s 4 and 12 to pltfl.dat and pltbnd.dat respectively c open(4,file='pltfl.dat',status='unknown',err=8888) open(12,file='pltbnd.dat',status='unknown',err=8888) c c Assign other tape no.s for unformatted scratch files c OPEN(1,FILE='temp1',STATUS='UNKNOWN',FORM='unformatted') OPEN(2,FILE='temp2',STATUS='UNKNOWN',FORM='unformatted') OPEN(3,FILE='temp3',STATUS='UNKNOWN',FORM='unformatted') OPEN(7,FILE='temp7',STATUS='UNKNOWN',FORM='unformatted') OPEN(8,FILE='temp8',STATUS='UNKNOWN',FORM='unformatted') OPEN(9,FILE='temp9',STATUS='UNKNOWN',FORM='unformatted') OPEN(10,FILE='temp10',STATUS='UNKNOWN',FORM='unformatted') OPEN(11,FILE='temp11',STATUS='UNKNOWN',FORM='unformatted') C C H E A D I N G C WRITE(LW,99) 99 FORMAT(1H1) WRITE(LW,999) 999 FORMAT(/,80(1H*),/) READ(LR,10)(TITLE(I),I=1,20) 10 FORMAT(20A4) WRITE(LW,11)(TITLE(I),I=1,20) 11 FORMAT(1X,20A4) WRITE(LW,999) C C RECOVER CPU TIME AND WRITE TO OUTPUT FILE C call time(tt) WRITE(LW,331)tt 331 FORMAT(' CPU TIME1=',A8) C C G R I D P A R A M E T E R S C c READ(LR,20)NE,NDS,NBN,NLOAD,NDISP,KEY c 20 FORMAT(/,3I6,3I3,/) READ(LR,*)NE,NDS,NBN,NLOAD,NDISP,KEY C WRITE(LW,22)NE,NDS,NBN,NLOAD,NDISP,KEY 22 FORMAT(//,5X,'G R I D P A R A M E T E R S',/, 1' NUMBER OF ELEMENTS .........................',I6,/, 2' NUMBER OF NODES ............................',I6,/, 3' NO. OF NODES ALONG CHOOSEN BOUNDARY ........',I6,/, 4' NUMBER OF NODE LOADS .......................',I3,/, 5' NO. OF PRESCRIBED NODE DISPLACEMENTS .......',I3,/, 6' OUTPUT FILE PRINT KEY ......................',I3) C C M A T E R I A L P R O P E R I T I E S C WRITE(LW,30) 30 FORMAT(//,10X,'M A T E R I A L P R O P E R T I E S',/) C c READ(LR,40)E1,E2,E3 READ(LR,*)E1,E2,E3 c READ(LR,40)G12,G13,G23 READ(LR,*)G12,G13,G23 c READ(LR,40)NU12,NU13,NU23 READ(LR,*)NU12,NU13,NU23 c READ(LR,41)AL1,AL2,AL3 READ(LR,*)AL1,AL2,AL3 c READ(LR,42)THETA READ(LR,*)THETA c 40 FORMAT(3(6X,E10.3)) c 41 FORMAT(3(8X,E10.3)) c 42 FORMAT(8X,E10.3,/) WRITE(LW,43)E1,E2,E3 43 FORMAT(/,' YOUNGS MODULI E1=',E10.3,/, * 14X,' E2=',E10.3,/, * 14X,' E3=',E10.3,/) WRITE(LW,44)G12,G13,G23 44 FORMAT(/,' SHEAR MODULI G12=',E10.3,/, * 13X,' G13=',E10.3,/, * 13X,' G23=',E10.3,/) WRITE(LW,45)NU12,NU13,NU23 45 FORMAT(/,' POISSONS RATIOS NU12=',E10.3,/, * 17X,'NU13=',E10.3,/, * 17X,'NU23=',E10.3,/) WRITE(LW,51)AL1,AL2,AL3 51 FORMAT(/,' THERMAL EXPANSION COEFFICIENTS ALPHA1=',E10.3,/, * 32X,'ALPHA2=',E10.3,/, * 32X,'ALPHA3=',E10.3,/) WRITE(LW,52)THETA 52 FORMAT(/,' FIBER ORIENTATION FROM H-AXIS, THETA(DEG)=',F6.2) C C L O A D I N F O R M A T I O N C O T H E R T H E N L O A D S A T N O D E P O I N T S C WRITE(LW,60) 60 FORMAT(1H1,10X,'L O A D I N F O R M A T I O N',/,2X, *'O T H E R T H E N L O A D S A T N O D E P O I N T S', *//) C c READ(LR,61)ITYPE READ(LR,*)ITYPE c 61 FORMAT(I6) c READ(LR,61)LTYPE READ(LR,*)LTYPE c READ(LR,62)TLOAD READ(LR,*)TLOAD c 62 FORMAT(F12.0,/) C WRITE(LW,64)LTYPE,ITYPE 64 FORMAT(/,10X,'LOAD TYPE = ',I1, */,23X,' ( = 1, MECHANICAL LOADS IN THE H-Z PLANE ONLY)', */,23X,' ( = 2, CHANGE IN THERMAL CONDITIONS SUPERPOSED)',/ */,10X,'CONSTRAINT TYPE = ',I1, */,23X,' ( = 1, PLANE STRAIN)', */,23X,' ( = 2, PLANE STRESS)',//) IF(LTYPE.EQ.2)WRITE(LW,67)TLOAD 67 FORMAT(' NONMECHANICAL LOADS: TEMPERATURE CHANGE, TLOAD=',E11.4) C C T I M E I N T E G R A T I O N I N F O R M A T I O N C C NEWMARK CONSTANTS AND PERIODIC WAVE PROPERTIES ON BOUNDARY C c READ(LR,68)DT,ALT,DELT,STG,FTG,IPRG,STB,FTB,IPRB READ(LR,*)DT READ(LR,*)ALT READ(LR,*)DELT READ(LR,*)STG READ(LR,*)FTG READ(LR,*)IPRG READ(LR,*)STB READ(LR,*)FTB READ(LR,*)IPRB c 68 FORMAT(5(1X,E11.4,/),1X,I11,/,2(1X,E11.4,/),1X,I11,/) WRITE(LW,69)DT,ALT,DELT,STG,FTG,IPRG,STB,FTB,IPRB 69 FORMAT(' T I M E I N T E G R A T I O N I N F O R M A T I O N', *//,' NEWMARK TIME-INTEGRATION CONSTANTS', */,' TIME STEP, DT=',E11.4, */,' CONSTANT, ALT=',E11.4, */,' CONSTANT, DELT=',E11.4,/ */,' WHOLE GRID:', */,' START PRINT TIME, STG=',E11.4, */,' FINAL TIME, FTG=',E11.4, */,' PRINT INTERVAL, IPRG=',I11,/ */,' BOUNDARY OF GRID:', */,' START PRINT TIME, STB=',E11.4, */,' FINAL TIME, FTG=',E11.4, */,' PRINT INTERVAL, IPRB=',I11,//) C C CALCULATE NEWMARK COEFFICIENTS AND WRITE TO OUTPUT FILE C A0=1./(ALT*DT**2) A2=1./(ALT*DT) A3=1./(2.*ALT)-1. A6=DT*(1.-DELT) A7=DT*DELT WRITE(LW,70)A0,A2,A3,A6,A7 70 FORMAT(' NEWMARK COEFF.S',/ * ,' A0=',E11.4,/ * ,' A2=',E11.4,/ * ,' A3=',E11.4,/ * ,' A6=',E11.4,/ * ,' A7=',E11.4,//) C C PERIODIC WAVE PROPERTIES ON BOUNDARY C c READ(LR,71)WF,WS,RHO,NWL READ(LR,*)WF READ(LR,*)WS READ(LR,*)RHO READ(LR,*)NWL c 71 FORMAT(3(1X,E11.4,/),1X,I11,/) C C CALCULATE OTHER WAVE PROPERTIES AND WRITE RESULTS TO OUTPUT FILE C PI=ACOS(-1.) WL=WS/WF WP=1./WF WC=2.*PI*WF ESZ=WL/NWL AMP=WL/(2.*PI) WRITE(LW,72)WF,WS,NWL,RHO,WC,WP,WL,ESZ,AMP 72 FORMAT(' PERIODIC WAVE PROPERTIS AT BOUNDARY',/ * ,' WAVE FREQUENCY, WF=',E11.4,/ * ,' WAVE SPEED, WS=',E11.4,/ * ,' ELEMENTS/WAVE LENGTH, NWL=',I3,/ * ,' MATERIAL DENSITY, RHO=',E11.4,/ * ,' WAVE TIME-CONSTANT, WC=',E11.4,/ * ,' WAVE PERIOD, WP=',E11.4,/ * ,' WAVE LENGTH, WF=',E11.4,/ * ,' ELEMENT SIZE, ESZ=',E11.4,/ * ,' PEAK DISPL. AMPLT., AMP=',E11.4,//) C C WRITE INFO. FOR PLOTTING WAVE PROPAGATION ON TAPE4 AND TAPE12 C XSCALP=0.0556 YSCALP=0.0556 C XSCALP=0.1667 C YSCALP=0.1667 XORIG=0.50 YORIG=0.50 WRITE(4,11)(TITLE(I),I=1,20) WRITE(4,24)NDS,NE,XSCALP,YSCALP,XORIG,YORIG 24 FORMAT(2I6,4F8.4) XBND=0.5 YBND=3.5 WRITE(12,11)(TITLE(I),I=1,20) WRITE(12,25)NDS,NBN,XSCALP,YSCALP,XBND,YBND 25 FORMAT(2I6,4F8.4) C C N O D E - E L E M E N T I N F O R M A T I O N C C ADDRESSES N1=1 N2=N1+2*NDS N3=N2+NDS C C NODE BOUNDARY CONDITION CODES ID(NDS,2) C NODE HORIZONTAL CO-ORDINATES H(NDS) C NODE VERTICAL CO-ORDINATES Z(NDS) C C ID(NDS,2) H(NDS) Z(NDS) CALL NEDATA(SHARE(N1),SHARE(N2),SHARE(N3),NDS,NE,NBN,ESZ,LTYPE, *TLOAD) C C WRITE INFO. FOR PLOTTING WAVE PROPAGATION ON TAPE4 AND TAPE12 C SCALG=1.0/ESZ SCALB=15.0*SCALG NBINC=(FTB-STB)/DT/IPRB+1 WRITE(12,26)NBINC,SCALB 26 FORMAT(I3,E10.3) ENDFILE 12 NGINC=(FTG-STG)/DT/IPRG+1 WRITE(4,26)NGINC,SCALG ENDFILE 4 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * C * F. E. M. S O L U T I O N * C * * * * * * * * * * * * * * * * * * * * * * * * * * * C C A S S I G N E Q U A T I O N N U M B E R S TO C E A C H D E G R E E O F F R E E D O M C STORED IN ARRAY ID(NDS,2) C C ID(NDS,2) CALL EQNUM(SHARE(N1),NDS,NEQ) C C C A L C U L A T E A D D R E S S E S U S E D T O S T O R E C T H E G L O B A L S T I F F N E S S A R R A Y C C GLOBAL ARRAY SIZE, NCOFS (NO. OF COEFFICIENTS) C GLOBAL ADDRESSES OF DIAGONAL TERMS, MAXA(NEQ1) C GLOBAL HALF BANDWIDTH INCLUDING DIAGONAL, MBAND C C CALCULATE NO. EQN.S PLUS ONE, NEQ1 USED TO C STORE DIAGONAL TERMS OF GLOBAL STIFFNESS ARRAY NEQ1=NEQ+1 C C ID(NDS,2) MAXA(NEQ1) CALL ADDRS(SHARE(N1),SHARE(N2),NE,NDS,NEQ1,MBAND,NCOFS) C C C A L C U L A T E T O T A L S I Z E O F A R R A Y S C T O B E S T O R E D I N C O M M O N C C ARRAYS USED IN SUBROUTINE QUAD4 AND SHARED BY MAIN PROG. COMMON C C ARRAY NAMES COMMON ADDRESSES C ID( NDS,2 ) = SHARE(N1) N1=1 C MAXA( NEQ1 ) = SHARE(N2) N2=N1+2*NDS C A( NCOFS ) = SHARE(N3) N3=N2+NEQ1 C V( NEQ ) = SHARE(N4) N4=N3+NCOFS C MA( NCOFS ) = SHARE(N5) N5=N4+NEQ C NTOT1=2*NDS+NEQ1+2.*NCOFS+NEQ C C ARRAYS USED IN SUBROUTINE PDVA AND SHARED BY MAIN PROG. COMMON C C ARRAY NAMES COMMON ADDRESSES C ID( NDS,2 ) = SHARE(N1) N1=1 C MAXA( NEQ1 ) = SHARE(N2) N2=N1+NDS*2 C A( NCOFS ) = SHARE(N3) N3=N2+NEQ1 C V( NEQ ) = SHARE(N4) N4=N3+NCOFS C MA( NCOFS ) = SHARE(N5) N5=N4+NEQ C U( NEQ ) = SHARE(N6) N6=N5+NCOFS C UD( NEQ ) = SHARE(N7) N7=N6+NEQ C UDD( NEQ ) = SHARE(N8) N8=N7+NEQ C NTOT2=2*NDS+NEQ1+4*NEQ+2*NCOFS C C BEFORE PROCEEDING COMPARE EACH SUBR. ARRAY SIZE WITH AVAILABLE MEMORY C C WRITE(LW,99) WRITE(LW,999) WRITE(LW,73)MTOT,NTOT1,NTOT2 73 FORMAT(' TOTAL SIZE OF ARRAYS STORED IN COMMON',/, 1' AVAILBLE(',I9,') REQUIRED BY SUBROUTINE QUAD4 (',I9,')',/, 233X,'SUBROUTINE PDVA (',I9,')') WRITE(LW,999) C IF(NTOT1.GT.MTOT)STOP IF(NTOT2.GT.MTOT)STOP C C C A L C U L A T E V A R I O U S E L E M E N T (QUAD4) A R R A Y S C CONNECTIVITY STIFFNESS FORCE MASS C LM( NDOFE ) AK( NDOFE,NDOFE ) F( NDOFE ) MK( NDOFE,NDOFE ) C C ---------- AND ---------- C C A S S E M B L E E L E M E N T A R R A Y S I N T O C G L O B A L A R R A Y S C STIFFNESS LOAD MASS C A( NCOFS ) V( NEQ ) MA( NCOFS ) C C ADDRESSES N3=N2+NEQ1 N4=N3+NCOFS N5=N4+NEQ N6=N5+NCOFS N7=N6+NEQ N8=N7+NEQ C C ID(NDS,2) MAXA(NEQ1) A(NCOFS) V(NEQ) MA(NCOFS) CALL QUAD4(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),SHARE(N5) * ,NDS,NEQ1,NCOFS,NEQ,ITYPE,NE,DD) call time(tt) WRITE(LW,332)tt 332 FORMAT(' CPU TIME2=',A8) C C B E G I N T I M E L O O P O N T R A N S I E N T S O L U T I O N C T=0.0E+00 ICNTG=0 ICNTB=0 NTINCG=FTG/DT+1 NTINCB=FTB/DT+1 ISTINCG=STG/DT+1 ISTINCB=STB/DT+1 DO 100 ITINC=1,NTINCG IF(ITINC.EQ.ISTINCG)ICNTG=1 IF(ITINC.EQ.ISTINCB)ICNTB=1 C C CALCULATE P RESCRIDED D ISPLACEMENTS, V ELOCITIES AND A CCELERATIONS C FOR EACH TIME STEP C ID(NDS,2) MAXA(NEQ1) A(NCOFS) V(NEQ) MA(NCOFS) CALL PDVA(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),SHARE(N5), *SHARE(N6),SHARE(N7),SHARE(N8) C U(NEQ) UD(NEQ) UDD(NEQ) * ,NDS,NEQ1,NCOFS,NEQ,MBAND,ESZ,NDISP,T) C C********** S O L V E S Y S T E M O F E Q N. S ************* C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C CALL COLSOL TWICE; 1ST TO TRIANGULARIZE A(NCOFS) C 2ND TO SOLVE FOR DISPLACEMENTS C LOADS STORED IN V(NEQ) REPLACED BY DISPLACEMENTS C C MAXA(NEQ1) A(NCOFS) V(NEQ) CALL COLSOL(SHARE(N2),SHARE(N3),SHARE(N4),NEQ,NCOFS,NEQ1,1) CALL COLSOL(SHARE(N2),SHARE(N3),SHARE(N4),NEQ,NCOFS,NEQ1,2) C C************* S Y S T E M O F E Q N. S S O L V E D ************ C C I F D I S P L A C E M E N T S W H E R E P R E S C R I B E D C R E C O V E R C O R R E S P O N D I N G N O D E F O R C E S C C C R E C O V E R D I S P L A C E M E N T S D(NDS,2) C F R O M COLSOL S O L T I O N S S T O R E D I N A R R A Y V(NEQ) C A N D C A L C U L A T E N E W V E L O C I T I E S A N D C A C C E L E R A T I O N S A T E A C H T I M E S T E P C C ID(NDS,2) D(NDS,2) V(NEQ) TUDD(NEQ) CALL DISP(SHARE(N1),SHARE(N3),SHARE(N4),SHARE(N5) * ,SHARE(N6),SHARE(N7),SHARE(N8),NDS,NBN,NEQ,T,ITINC,ICNTG,ICNTB) C U(NEQ) UD(NEQ) UDD(NEQ) C C RECOVER CPU TIME AND WRITE TO OUTPUT FILE C call time(tt) WRITE(LW,333)T,tt,ICNTG,ICNTB 333 FORMAT(' T=',E11.4,' CPU TIME=',A8,' ICNTG=',I2,' ICNTB=',I2) C C I N C R E M E N T T I M E S T E P C T=T+DT IF(ITINC.GE.ISTINCG)THEN IF(ICNTG.EQ.IPRG)ICNTG=0 ICNTG=ICNTG+1 ENDIF IF(ITINC.GE.ISTINCB)THEN IF(ICNTB.EQ.IPRB)ICNTB=0 ICNTB=ICNTB+1 ENDIF 100 CONTINUE C C E N D O F T I M E L O O P C 8888 STOP END C*********************************************************************** SUBROUTINE NEDATA(ID,H,Z,NDS,NE,NBN,ESZ,LTYPE,TLOAD) C CALLING PROGRAM C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * READS NODE - ELEMENT DATA FROM TAPE UNIT LR. * C * STORES PARAMETERS FOR REFERENCE ON TAPES (1)AND(3) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C REAL NU12,NU13,NU23 C COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E1,E2,E3,NU12,NU13,NU23,G12,G13,G23,AL1,AL2,AL3,THETA COMMON/RDK4/IBN(200) C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),H(NDS),Z(NDS) C C CALCULATE A UNIFORM THERMAL STRAIN FOR ALL ELEMENTS IN THE H-Z PLANE C TS1=0.0 TS2=0.0 TS3=0.0 IF(LTYPE.EQ.2)TS1=AL1*TLOAD IF(LTYPE.EQ.2)TS2=AL2*TLOAD IF(LTYPE.EQ.2)TS3=AL3*TLOAD C C R E A D N O D E D A T A C C WRITE HEADER ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,10) 10 FORMAT(1H1,5X,'N O D E I N F O R M A T I O N',//, 14X,'(FIXED=1',/,6X,'FREE=0)',/,' NODE MOTION CO-ORDINATES',/, 2' NO. H Z ',6X,'H',8X,'Z') C C L O O P T H R U A L L N O D E S C DO 40 I=1,NDS C C READ NODE DATA FROM INPUT FILE C READ(LR,20)ID(I,1),ID(I,2),H(I),Z(I) 20 FORMAT(6X,2I3,2F12.0) C C WRITE UNDEFORMED UNSCALED GRID COORD.S ON TAPE4 AND TAPE12 C WRITE(4,21)I,ID(I,1),ID(I,2),H(I),Z(I) 21 FORMAT(I6,2I3,2F12.5) WRITE(12,21)I,ID(I,1),ID(I,2),H(I),Z(I) C C WRITE NODE DATA ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,30)I,ID(I,1),ID(I,2) * ,H(I),Z(I) 30 FORMAT(1X,I4,2(1X,I1,1X),1X,2F9.5) C 40 CONTINUE C C E N D N O D E L O O P C C C R E A D - S T O R E E L E M E N T D A T A C C C A L C U L A T E : RESCALED NODE CO-ORDINATES C ASSOCIATED WITH EACH ELEMENT C C C A L C U L A T E - S T O R E 1. ELEMENT NODE POINT COORDINATES C C C WRITE HEADER ON OUTPUT FILE C 45 IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,50) 50 FORMAT(1H1,5X,'E L E M E N T I N F O R M A T I O N',//, 1' ELMT--------NODES--------',/, 2' NO. 1ST 2ND 3RD 4TH THICKNESS') C C L O O P T H R U A L L E L E M E N T S C DO 90 K=1,NE C C READ ELEMENT DATA FROM INPUT FILE C READ(LR,60)ND1,ND2,ND3,ND4,THIC 60 FORMAT(6X,4I6,F12.0) C C WRITE GRID TOPOLOGY ON TAPE4 C WRITE(4,61)K,ND1,ND2,ND3,ND4 61 FORMAT(5I6) C C STORE NODES ASSOCIATED WITH EACH ELEMENT ON TAPE 1 C WRITE(1)ND1,ND2,ND3,ND4,THIC,TS1,TS2,TS3 C C NODE CO-ORDINATES ASSOCIATED WITH EACH ELEMENT RESCALED C H1=H(ND1)*ESZ H2=H(ND2)*ESZ H3=H(ND3)*ESZ H4=H(ND4)*ESZ Z1=Z(ND1)*ESZ Z2=Z(ND2)*ESZ Z3=Z(ND3)*ESZ Z4=Z(ND4)*ESZ C C STORE ELEMENT NODE POINT COORD.S ON TAPE 3 C WRITE(3)H1,H2,H3,H4,Z1,Z2,Z3,Z4 C C WRITE ELEMENT DATA ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,70)K,ND1,ND2,ND3,ND4,THIC 70 FORMAT(5(1X,I4),E12.5) C 90 CONTINUE C C READ LINE OF NODES WHERE DISPL.S ARE STORED C READ(LR,130)IDUM DO 120 IBND=1,NBN READ(LR,100)IBN(IBND) 100 FORMAT(I6) 120 CONTINUE READ(LR,130)IDUM 130 FORMAT(A4) C C E N D E L E M E N T L O O P C RETURN END C*********************************************************************** SUBROUTINE EQNUM(ID,NDS,NEQ) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * GENERATES EQUATION NUMBERS BY REPLACING THE FREE DISPL * C * BOUNDARY CONDITION CODE (ID=0) WITH AN EQUATION NUMBER. * C * EQUATION NUMBERS ASSOCIATED WITH EACH FREE NODE DEGREE * C * OF FREEDOM (DOF) ARE STORED IN ARRAY ID(NDS,2) WHICH IS * C * SHARED WITH MAIN PROGRAM CONNOM SHARE(MTOT) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLE SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2) C C INITIALIZE VARIABLES DEFINING EQUATION NUMBERS C IEQN=0 C C WRITE HEADER ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,10) 10 FORMAT(1H1,13X,'EQUATION NUMBERS',//,5(12X,'NODE MOTION '),/, * 5(' # NODE HORIZ VERT')) C C SEQUENTIALLY TEST THE DISPLACEMENT BOUNDARY CONDITION C CODES( FIXED=1 OR FREE=0 ) FOR EACH NODE. THERE ARE C TWO POSSIBLE DEGREES OF FREEDOM ASSOCIATED WITH EACH NODE. C ID(I,1)=0, I-TH NODE IS FREE TO MOVE IN THE HORIZ. H-DIRECTION C ID(I,2)=0, I-TH NODE IS FREE TO MOVE IN THE VERT. Z-DIRECTION C IF A NODE IS FREE TO MOVE THEN AN EQUATION NUMBER (IEQN) IS C ASSOCIATED WITH THAT NODE DEGREE OF FREEDOM. C DO 40 I=1,NDS DO 30 J=1,2 IF(ID(I,J).EQ.1)GO TO 20 IEQN=IEQN+1 ID(I,J)=IEQN GO TO 30 20 ID(I,J)=0 30 CONTINUE 40 CONTINUE C C WRITE EQN. NO.S ON OUTPUT FILE C IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(LW,3)(I,(ID(I,J),J=1,2),I=1,NDS) 3 FORMAT(5(' #',I5,6X,2I6)) C C LAST VALUE OF IEQN IS EQUAL TO THE TOTAL NO. OF EQUATIONS C NEQ=IEQN C C WRITE TOTAL NUMBER OF EQUATIONS DUE TO NODAL D.O.F. S C ON OUTPUT FILE C WRITE(LW,45)NEQ 45 FORMAT(//,' TOTAL NO. OF EQN.S DUE TO NODAL DEGREES OF FREEDOM', *I5,//) C RETURN END C*********************************************************************** SUBROUTINE ADDRS(ID,MAXA,NE,NDS,NEQ1,MBAND,NCOFS) C C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * 1. DETERMINES COLUMN HEIGHTS, STORED TEMPORARILY * C * IN ARRAY MAXA(NEQ). * C * 2. DETERMINES ADDRESSES = ( I-TH SUBSCRIPT IN A(I) ) FOR * C * THE DIAGONAL TERMS IN THE GLOBAL ARRAY A(I) * C * ADDRESSES STORED IN MAXA(NEQ+1). * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES ID AND MAXA SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),MAXA(NEQ1),LM(8) C C CLEAR LM AND MAXA ARRAYS C DO 5 IEQN=1,NEQ1 5 MAXA(IEQN)=0 DO 10 I=1,8 10 LM(I)=0 C C USING EQN. NO.S STORED IN ID(NDS,2) ARRAY TOGETHER WITH C NODE NO.S (ND1,ND2,ND3,ND4) STORED ON TAPE 1 C CALCULATE CONNECTIVITY LM ARRAYS FOR EACH ELEMENT. C REWIND 1 C C L O O P T H R U A L L E L E M E N T S C DO 40 K=1,NE C C READ NODE NO.S ASSOCIATED WITH K-TH ELEMENT C READ(1)ND1,ND2,ND3,ND4,THIC,TS1,TS2,TS3 C C ASSIGN EQUATION NUMBERS TO ELEMENT CONNECTIVITY ARRAY C C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 1ST NODE)=ID(ND1,1)=LM(1) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 1ST NODE)=ID(ND1,2)=LM(2) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 2ND NODE)=ID(ND2,1)=LM(3) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 2ND NODE)=ID(ND2,2)=LM(4) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 3RD NODE)=ID(ND3,1)=LM(5) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 3RD NODE)=ID(ND3,2)=LM(6) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 4TH NODE)=ID(ND4,1)=LM(7) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 4TH NODE)=ID(ND4,2)=LM(8) C LM(1)=ID(ND1,1) LM(2)=ID(ND1,2) LM(3)=ID(ND2,1) LM(4)=ID(ND2,2) LM(5)=ID(ND3,1) LM(6)=ID(ND3,2) LM(7)=ID(ND4,1) LM(8)=ID(ND4,2) C C CHOOSE NO. OF DEGREES OF FREEDOM PER ELEMENT, NDOFE C NDOFE=8 C C CALCULATE THE COLUMN HEIGHTS ABOVE THE DIAGONAL C OF THE ASSEMBLED GLOBAL STIFFNESS ARRAY C (COLUMN HEIGHTS TEMPORARILY STORED IN MAXA ARRAY) C C THIS IS DONE BY: C COMPARING VALUES IN LM ARRAYS FOR ELEMENTS C WHICH SHARE A DEGREE OF FREEDOM(DOF) C COLUMN HEIGHT = MAXA( SHARED DOF ) - LOWEST DOF IN SHARED LM ARRAY C DO 30 I=1,NDOFE L=LM(I) IF(L.EQ.0)GO TO 30 DO 20 J=1,NDOFE M=LM(J) IF(M.EQ.0)GO TO 20 IF(M.GT.L)GO TO 20 N=L-M+1 IF(N.GT.MAXA(L))MAXA(L)=N 20 CONTINUE 30 CONTINUE 40 CONTINUE C C E N D E L E M E N T L O O P C C CALCUALTE WIDTH OF HALFBAND PLUS DIAGONAL, MBAND C WRITE COLUMN HEIGHTS AND MBAND ON OUTPUT FILE C MBAND=0 NEQ=NEQ1-1 DO 55 I=1,NEQ IF(MAXA(I).GT.MBAND)MBAND=MAXA(I) 55 CONTINUE IF((KEY.EQ.4).OR.(KEY.EQ.7))WRITE(LW,50)(I,MAXA(I),I=1,NEQ1) 50 FORMAT(4(' # COLUMN=',I6,' HEIGHT=',I6,' #')) WRITE(LW,56)MBAND 56 FORMAT(' WIDTH OF HALFBAND PLUS DIAGONAL, MBAND=',I9) C C USING COLUMN HEIGHTS STORED IN MAXA( NEQ ) C DETERMINE ADDRESSES = ( I-TH SUBSCRIPT IN A(I) ) C FOR DIAGONAL TERMS IN THE GLOBAL ARRAY A(I) C AND STORE VALUES IN MAXA( NEQ + 1 ). C C FIRST COLUMN DIAGONAL ADDRESS BY DEFINITION IS ONE C NN=1 MAXA(1)=1 MAXM1=1 C C ADD COLUMN HEIGHTS TO PREVIOUS DIAGONAL ADDRESS C DO 60 I=2,NEQ1 NN=NN+MAXM1 MAXM1=MAXA(I) MAXA(I)=NN 60 CONTINUE C C WRITE GLOBAL ARRAY SIZE AND DIAGONAL ADDRESSES ON THE OUTPUT FILE C NCOFS=MAXA(NEQ1) WRITE(LW,70)NCOFS 70 FORMAT(/,' NO. OF COEFF.S IN GLOBAL STIFFNESS ARRAY, NCOFS=',I9,/) C IF((KEY.EQ.6).OR.(KEY.EQ.7))WRITE(LW,75)(I,MAXA(I),I=1,NEQ) 75 FORMAT(3(' # COLUMN=',I6,' DIAGONAL ADDRESS=',I6,' #')) C RETURN END ************************************************************************ SUBROUTINE QUAD4(ID,MAXA,A,V,MA,NDS,NEQ1,NCOFS,NEQ,ITYPE,NE,DD) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M - QUAD4 * C * * C * CALCULATES THE STIFFNESS AND MASS MATRICES FOR ALL FOUR-NODED * C * QUADRILATERAL ELEMENTS IN PLANE STRESS OR PLANE STRAIN CONDITIONS * * C * * C * - - INPUT VARIABLES - - * C * NE = NUMBER OF ELEMENTS * C * NDS = NUMBER OF NODES * C * NCOFS = NUMBER OF COEFFICIENTS IN GLOBAL ARRAYS * C * NEQ = NUMBER OF EQUATIONS * C * NEQ1 = NEQ+1 * C * C = NODE POINT COORDINATES * C * ITYPE = ELEMENT TYPE * C * EQ.1 = PLANE STRAIN * C * EQ.2 = PLANE STRESS * C * NINT = GAUSS NUMERICAL INTEGRATION ORDER (MAX. 4) * C * THIC = THICKNESS OF ELEMENT * C * E1 , E2 , E3 = YOUNG'S MODULUI * C * G12, G13, G23 = SHEAR MODULUI * C * NU12,NU13,NU23 = POISSON'S RATIOS * C * RHO = MATERIAL DENSITY * C * * C * - - CALCULATES - - * C * AK(NDOFE,NDOFE) = ELEMENT STIFFNESS MATRIX * C * MK(NDOFE,NDOFE) = ELEMENT MASS MATRIX * C * F(NDOFE) = ELEMENT LOAD MATRIX * C * NDOFE = NO. OF DEGREES OF FREEDOM * C * * C * - - OUTPUT (CALCULATES AND RETURNS) * C * A(NCOFS) = GLOBAL STIFFNESS MATRIX * C * MA(NCOFS) = GLOBAL MASS MATRIX * C * V(NEQ) = GLOBAL LOAD VECTOR * C * DD(3,3) = MATERIAL ELASTIC STIFFNESS MATRIX * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT REAL*8(A-H,O-Z) C SQRT(X)=DSQRT(X) C ABS(X)=DABS(X) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON CDC * C * EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IMB OR UNIVAC * C * MACHINES. FOR SINGLE OR DOUBLE PRECISION ARITHMETIC ACTIVATE * C * DEACTIVATE OR ADJUST ABOVE CARDS * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C REAL NU12,NU13,NU23,NU21,NU31,NU32,MK,MA,MASS COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E1,E2,E3,G12,G13,G23,NU12,NU13,NU23,AL1,AL2,AL3,THETA COMMON/RDK3/A0,A2,A3,A6,A7,WC,RHO,AMP DIMENSION ID(NDS,2),MAXA(NEQ1),A(NCOFS),V(NEQ),MA(NCOFS) DIMENSION XG(4,4),WGT(4,4),DD(3,3),C(2,4),H(2,8),B(3,8),DB(3) DIMENSION AK(8,8),F(8),LM(8),MK(8,8) C C ZERO ALL ELEMENTS OF XG NAD WXG ARRAYS TO 15 SIGNIFICANT FIGURES C DO 5 I=1,4 DO 5 J=1,4 XG(I,J)=0.000000000000000 WGT(I,J)=0.000000000000000 5 CONTINUE C C MATRIX XG STORES GAUSS - LEGENDRE SAMPLING POINTS C XG(1,2)=-0.577350269189626 XG(2,2)=-XG(1,2) XG(1,3)=-0.774596669241483 XG(3,3)=-XG(1,3) XG(1,4)=-0.861136311594053 XG(2,4)=-0.339981043584856 XG(3,4)=-XG(2,4) XG(4,4)=-XG(1,4) C C MATRIX WGT STORESS GAUSS - LEGENDRE WEIGHTING FACTORS C WGT(1,1)=2.000000000000000 WGT(1,2)=1.000000000000000 WGT(2,2)=WGT(1,2) WGT(1,3)=0.555555555555556 WGT(2,3)=0.888888888888889 WGT(3,3)=WGT(1,3) WGT(1,4)=0.347854845137454 WGT(2,4)=0.652145154862546 WGT(3,4)=WGT(2,4) WGT(4,4)=WGT(1,4) C C O B T A I N S T R E S S - S T R A I N L A W C F O R O R T H O R H O M B I C M A T E R I A L C C CALCULATE RECIPROCAL POISSION'S RATIOS C NU21=E2*NU12/E1 NU31=E3*NU13/E1 NU32=E3*NU23/E2 C C CALCULATE CONTRACTED STIFFNESSES FROM ENGINEERING ELASTIC CONSTANTS C DEL=1.-NU12*NU21-NU23*NU32-NU31*NU13-2.*NU21*NU32*NU13 C11=(1.-NU23*NU32)*E1/DEL C22=(1.-NU13*NU31)*E2/DEL C33=(1.-NU12*NU21)*E3/DEL C12=(NU21+NU31*NU23)*E1/DEL C13=(NU31+NU21*NU32)*E1/DEL C23=(NU32+NU12*NU31)*E2/DEL C44=G23 C55=G13 C66=G12 C C CALCULATE STIFFNESSES ROTATED THETA-DEGREES FROM FIBER H-AXIS C PI=ACOS(-1.) THRAD=THETA*PI/180. CS=COS(THRAD) SN=SIN(THRAD) C2=CS*CS C4=C2*C2 S2=SN*SN S4=S2*S2 CB11=C11*C4+2.*C2*S2*(C12+2.*C66)+C22*S4 CB12=C2*S2*(C11+C22-4.*C66)+C12*(S4+C4) CB13=C13*C2+C23*S2 CB16=-CS*SN*(C11*C2-C22*S2-(C12+2.*C66)*(C2-S2)) CB22=C11*S4+2.*C2*S2*(C12+2.*C66)+C22*C4 CB23=C13*S2+C23*C2 CB26=-CS*SN*(C11*S2-C22*C2+(C12+2.*C66)*(C2-S2)) CB33=C33 CB36=CS*SN*(C23-C13) CB66=(C11+C22-2.*C12)*C2*S2+C66*(C2-S2)**2 C C FOR PLAIN STRAIN AYALYSIS C DD(1,1)=CB11 DD(1,2)=CB12 DD(1,3)=CB16 DD(2,1)=CB12 DD(2,2)=CB22 DD(2,3)=CB26 DD(3,1)=CB16 DD(3,2)=CB26 DD(3,3)=CB66 IF(ITYPE.EQ.1)THIC=1. IF(ITYPE.EQ.1)GO TO 10 C C FOR PLANE STRESS ANALYSIS C DD(1,1)=CB11-CB13**2/CB33 DD(2,2)=CB22-CB23**2/CB33 DD(3,3)=CB66-CB36**2/CB33 DD(1,2)=CB12-CB13*CB23/CB33 DD(2,1)=DD(1,2) DD(2,3)=CB26-CB23*CB36/CB33 DD(3,2)=DD(2,3) DD(1,3)=CB16-CB13*CB36/CB33 DD(3,1)=DD(1,3) C C WRITE STIFFNESS MATRIX ON OUTPUT FILE C WRITE(LW,6) 6 FORMAT(//,' ELASTIC STIFFNESS MATRIX FOR EACH ELEMENT',/) DO 7 I=1,3 7 WRITE(LW,8)(1,J,DD(I,J),J=1,3) 8 FORMAT(3(' DD(',I1,',',I1,')=',E11.4)) C C CLEAR GLOBAL ARRAYS C 10 DO 1 I=1,NCOFS MA(I)=0.0 1 A(I)=0.0 DO 2 I=1,NEQ 2 V(I)=0.0 C C REWIND TAPES 1 AND 3 C REWIND 1 REWIND 3 C C L O O P T H R U A L L E L E M E N T S C DO 100 NEL=1,NE C C WRITE(LW,21)NEL C 21 FORMAT(1X,80(1H*),/,' ELEMENT NO.=',I4) C C READ PARAMETERS FOR EACH ELEMENT C READ(1)ND1,ND2,ND3,ND4,THIC,TS1,TS2,TS3 READ(3)C(1,1),C(1,2),C(1,3),C(1,4),C(2,1),C(2,2),C(2,3),C(2,4) C C CREATE CONNECTIVITY ARRAY LM(NDOFE) FOR EACH ELEMENT C NDOFE=8 LM(1)=ID(ND1,1) LM(2)=ID(ND1,2) LM(3)=ID(ND2,1) LM(4)=ID(ND2,2) LM(5)=ID(ND3,1) LM(6)=ID(ND3,2) LM(7)=ID(ND4,1) LM(8)=ID(ND4,2) C C C A L C U L A T E E L E M E N T S T I F F N E S S C C ZERO ELEMENT ARRAYS C 20 DO 30 I=1,NDOFE F(I)=0.0 DO 30 J=1,NDOFE MK(I,J)=0.0 30 AK(I,J)=0.0 C NINT=2 DO 80 LX=1,NINT RI=XG(LX,NINT) DO 80 LY=1,NINT SI=XG(LY,NINT) C C EVALUATE DERIVATIVE OPERATOR B FOR 4-NODED QUADRILATERAL ELEMENTS C CALL STDM(C,B,H,DET,RI,SI,NEL) C C STORE B(I,J) ARRAY ON TAPE2 C 25 WRITE(2)NEL,LX,LY,((B(I,J),I=1,3),J=1,8) C C WITH B(I,J) ARRAY, CALCULATE ELEMENTAL STIFFNESS MATRIX C WT=WGT(LX,NINT)*WGT(LY,NINT)*THIC*DET C DO 70 J=1,NDOFE C DO 40 K=1,3 DB(K)=0.0 DO 40 L=1,3 DB(K)=DB(K)+DD(K,L)*B(L,J) 40 CONTINUE C DO 60 I=J,NDOFE C STIFF=0.0 DO 50 L=1,3 STIFF=STIFF+B(L,I)*DB(L) 50 CONTINUE C MASS=0.0 DO 51 L=1,2 MASS=MASS+H(L,I)*H(L,J) 51 CONTINUE C AK(I,J)=AK(I,J)+STIFF*WT MK(I,J)=MK(I,J)+MASS*WT*RHO C 60 CONTINUE C F(J)=F(J)+(B(1,J)*(DD(1,1)*TS1+DD(1,2)*TS2+DD(1,3)*TS3) * +B(2,J)*(DD(1,2)*TS1+DD(2,2)*TS2+DD(2,3)*TS3) * +B(3,J)*(DD(1,3)*TS1+DD(2,3)*TS2+DD(3,3)*TS3)) 70 CONTINUE 80 CONTINUE C C SYMMETRY C DO 90 J=1,NDOFE DO 90 I=1,NDOFE MK(J,I)=MK(I,J) 90 AK(J,I)=AK(I,J) C C ASSEMBLE EACH ELEMENT INTO THE GLOBAL MATRICES C 4 CALL ASSEMB(NDOFE,AK,F,LM,MK,MAXA,A,V,MA,NEQ1,NCOFS,NEQ) C 100 CONTINUE C C * * * * * * * * * * * * * * * * * * * C * E N D E L E M E N T L O O P * C * * * * * * * * * * * * * * * * * * * C C WRITE ASSEMBLED 1. GLOBAL STIFFNESS ARRAY ON OUTPUT FILE C 2. RIGHT-HAND-SIDE LOAD VECTOR ON OUTPUT FILE C IF((KEY.EQ.4).OR.(KEY.EQ.7))GO TO 74 GO TO 999 C C WRITE HEADER AND FIRST DIAGONAL TERM ON OUTPUT FILE C 74 WRITE(LW,75)A(1) 75 FORMAT(1H1,//,' GLOBAL STIFFNESS ARRAY',//,' COLUMN NO. 1', */,7X,' A( 1)=',E11.4,' 1 RST DIAGONAL POSITION') C C START LOOP AT 2 ND DIAGONAL TERM C IDIAG=2 C C LOOP THRU A(I) TERMS C DO 84 I=2,NCOFS IF(I.EQ.MAXA(IDIAG))GO TO 77 IF(A(I).EQ.0.0)GO TO 84 WRITE(LW,76)I,A(I) 76 FORMAT(8X,'A(',I4,')=',E14.7) GO TO 84 77 IF(I.NE.NCOFS)WRITE(LW,78)IDIAG,I,A(I),IDIAG 78 FORMAT(/,' COLUMN NO.',I4,/,8X,'A(',I4,')=',E14.7, *' DIAGONAL NO.',I4) C IF(I.EQ.NCOFS)WRITE(LW,79)I,A(I) 79 FORMAT(//,8X,'A(',I4,')=',E14.7,' LAST STORED VALUE IN A(I)') C IDIAG=IDIAG+1 C 84 CONTINUE C C WRITE RIGHT-HAND-SIDE LOAD VECTOR ON OUTPUT FILE C WRITE(LW,81) 81 FORMAT(1H1,//,' RIGHT-HAND-SIDE LOAD VECTOR',//) C INON=0 DO 82 I=1,NEQ IF(V(I).EQ.0.0)GO TO 82 INON=INON+1 WRITE(LW,83)I,V(I) 82 CONTINUE 83 FORMAT(7X,'V(',I4,')=',E14.7) IF(INON.EQ.0)WRITE(LW,85) 85 FORMAT(' ALL ELEMENTS OF THE LOAD VECTOR WERE ZERO') C C T H I S E N D S W R I T I N G G L O B A L A R R A Y S C C STORE GLOBAL STIFFNESS ARRAY ON TAPE3 C 999 REWIND 3 WRITE(3)(A(I),I=1,NCOFS) C C ADD THE GLOBAL MASS MATRIX TO THE STIFFNESS MATRIX FOR DYNAMIC SOLUTION C DO 110 I=1,NCOFS A(I)=A(I)+A0*MA(I) 110 CONTINUE C C STORE DYNAMIC GLOBAL ARRAYS ON TAPES 7,8,AND 9 FOR REF. C WRITE(7)(MA(I),I=1,NCOFS) WRITE(8)(A(I),I=1,NCOFS) WRITE(9)(V(I),I=1,NEQ) C RETURN END C*********************************************************************** SUBROUTINE STDM(XX,B,H,DET,R,S,NEL) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * 1. EVALUATE THE STRAIN-DISPLACEMENT TRANSFORMATION MATRIX B(3,8)* C * 2. EVALUATE THE DISPLACEMENT INTERPOLATION MATRIX H(2,8) * C * AT POINT (R,S) FOR A QUADRILATERAL 4-NODED ELEMENT * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C IMPLICIT REAL*8(A-H,O-Z) C COMMON/RDK1/LW,LR,KEY DIMENSION XX(2,8),B(3,8),Z(4),P(2,4),XJ(2,2),XJI(2,2),H(2,8) C C EVALUATE INTERPOLATION COEFF.S AND JACOBIAN TRANSFORMATION MATRICES C CALL JACOB(NEL,R,S,XX,Z,P,XJ,XJI,DET) C C CLEAR B(I,J) AND H(I,J) ARRAYS C DO 40 I=1,3 DO 40 J=1,8 DO 40 K=1,2 H(K,J)=0.0 B(I,J)=0.0 40 CONTINUE C C CONSTRUCT THE NONZERO TERMS OF H(2,8) FROM Z(4) C DO 45 K=1,4 J=1+2*(K-1) H(1,J)=Z(K) H(2,J+1)=Z(K) 45 CONTINUE C C EVALUATE GLOBAL DERIVATIVE OPERATOR B C K2=0 DO 60 K=1,4 K2=K2+2 B(1,K2-1)=0. B(1,K2)=0. B(2,K2-1)=0. B(2,K2)=0. DO 50 I=1,2 B(1,K2-1)=B(1,K2-1)+XJI(1,I)*P(I,K) 50 B(2,K2)=B(2,K2)+XJI(2,I)*P(I,K) B(3,K2)=B(1,K2-1) 60 B(3,K2-1)=B(2,K2) C RETURN END C*********************************************************************** SUBROUTINE JACOB(NEL,R,S,XX,H,P,XJ,XJI,DET) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATES INTERPOLATION POLYNOMIALS AND JACOBIAN TRANSFOR- * C * TION MATRICES FOR A 4-NODED QUADRILATERIAL ELEMENT. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY DIMENSION XX(2,4),H(4),P(2,4),XJ(2,2),XJI(2,2) C RP=1.0+R SP=1.0+S RM=1.0-R SM=1.0-S C C INTERPOLATION FUNCTIONS C H(1)=0.25*RP*SP H(2)=0.25*RM*SP H(3)=0.25*RM*SM H(4)=0.25*RP*SM C C NATURAL COORDINATE DERIVATIVES OF THE INTERPOLATION FUNCTIONS C C 1. WITH RESPECT TO R C P(1,1)=0.25*SP P(1,2)=-P(1,1) P(1,3)=-0.25*SM P(1,4)=-P(1,3) C C 2. WITH RESPECT TO S C P(2,1)=0.25*RP P(2,2)=0.25*RM P(2,3)=-P(2,2) P(2,4)=-P(2,1) C C EVALUATE THE JACOBIAN MATRIX AT POINT (R,S) C DO 30 I=1,2 DO 30 J=1,2 DUM=0.0 DO 20 K=1,4 20 DUM=DUM+P(I,K)*XX(J,K) 30 XJ(I,J)=DUM C C COMPUTE THE DETERMINANT OF THE JACOBIAN MATRIX AT POINT (R,S) C DET=XJ(1,1)*XJ(2,2)-XJ(2,1)*XJ(1,2) IF(DET.GT.0.00000001)GO TO 40 WRITE(LW,2000)NEL STOP C C COMPUTE INVERSE OF THE JACOBIAN MATRIX C 40 DUM=1./DET XJI(1,1)=+XJ(2,2)*DUM XJI(1,2)=-XJ(1,2)*DUM XJI(2,1)=-XJ(2,1)*DUM XJI(2,2)=+XJ(1,1)*DUM C RETURN C 2000 FORMAT(10H0*** ERROR, 1 52H ZERO OR NEGATIVE JACOBIAN DETERMINANT FOR ELEMNET (,I4, 2 1H) ) C END C*********************************************************************** SUBROUTINE ASSEMB(NDOFE,AK,F,LM,MK,MAXA,A,V,MA,NEQ1,NCOFS,NEQ) C CALLING PROGRAM - ELEMNT C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * ASSEMBLES : 1. EACH ELEMENT STIFFNESS ARRAY AK(I,J) * C * ONTO THE GLOBAL ARRAY A(I). * C * 2. EACH ELEMENT FORCE VECTOR F(I) ONTO THE * C * RIGHT-HAND-SIDE LOAD VECTOR V(I). * C * 3. EACH ELEMENT MASS ARRAY MK(I,J) * C * ONTO THE GLOBAL ARRAY MA(I). * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C REAL MA,MK DIMENSION AK(8,8),F(8),LM(8),MK(8,8) DIMENSION MAXA(NEQ1),A(NCOFS),V(NEQ),MA(NCOFS) C DO 20 I=1,NDOFE L=LM(I) IF(L.EQ.0)GO TO 20 NN=MAXA(L) DO 10 J=I,NDOFE M=LM(J) IF(M.EQ.0)GO TO 10 MM=MAXA(M) KK=M-L MM=MM+KK IF(KK.LT.0)MM=NN-KK C A(MM)=A(MM)+AK(I,J) MA(MM)=MA(MM)+MK(I,J) 10 CONTINUE C V(L)=V(L)+F(I) 20 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE PDVA(ID,MAXA,A,V,MA,U,UD,UDD,NDS,NEQ1,NCOFS,NEQ,MBAND, * ESZ,NDISP,T) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * 1. ADDS A HUGE SPRING STIFFNESS (HUGEK) TO THE ASSEMBLED * C * STIFFNESS (A(MAXA(IEQN))) AND LOAD (V(IEQN)) COMPONENTS * C * SO AS TO SIMULATE A PRESCRIBED NODAL DISPLACEMENT. * C * 2. UPDATES LOAD VECTOR WHERE INERTIAL TERMS ARE INCLUDED * C * FROM THE NEWMARK TIME INTEGRATION APPROXIMATION. * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C REAL MA COMMON/RDK1/LW,LR,KEY COMMON/RDK3/A0,A2,A3,A6,A7,WC,RHO,AMP C DIMENSION NOD(50),IDIRN(50) C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),MAXA(NEQ1),A(NCOFS),V(NEQ),MA(NCOFS) DIMENSION U(NEQ),UD(NEQ),UDD(NEQ) C C INTIALIZE DISP. VEL. AND ACCEL. TO ZERO AT TIME=0.0 C IF(T.GT.0.0)GO TO 5 DO 4 IEQN=1,NEQ U(IEQN)=0.0 UD(IEQN)=0.0 UDD(IEQN)=0.0 4 CONTINUE C C CALCULATE BANDWIDTH PROPERTIES C 5 MB=MBAND-1 LIMIT=NEQ1-MBAND C C READ GLOBAL ARRAYS FROM TAPES C REWIND 7 REWIND 8 REWIND 9 READ(7)(MA(I),I=1,NCOFS) READ(8)(A(I),I=1,NCOFS) READ(9)(V(I),I=1,NEQ) C C USING EQN NO.S STORED IN ID(NDS,2) ARRAY C RECOVER EQN NO. (IEQN) CORRESPONDING TO NODE DISPLACEMENT. C THEN 1. ADD LARGE SPRING STIFFNESS (HUGEK) TO DIAGONAL C TERM ( A(MAXA(IEQN)) ) OF THE GLOBAL STIFFNESS ARRAY. C 2. ADD LARGE SPRING FORCE ( HUGEK*AMP) TO THE C RIGHT-HAND-SIDE LOAD VECTOR ( V(IEQN) ). C HUGEK=1.0E+20 C C WRITE HEADER ON OUTPUT FILE C C WRITE(LW,10)NDISP C 10 FORMAT(1H1,//,' NUMBER OF PRESCRIBED DISPL.S=',I4,/) C C L O O P T H R U P R E S C R I B E D D I S P L. S C C C READ 1. NOD , NODE NO. WHERE I-TH DISPL. IS PRESCRIBED C 2. IDIRN, DIRECTION OF I-TH DISPL. C 3. DISP , MAGNITUDE OF EACH DISPL. C IF(T.GT.0.0)GO TO 20 DO 10 IDISP=1,NDISP 10 READ(LR,*)NOD(IDISP),IDIRN(IDISP) c 10 READ(LR,15)NOD(IDISP),IDIRN(IDISP) c 15 FORMAT(2I6) C 20 CONTINUE C WRITE(LW,71)T C 71 FORMAT(' TIME=',E11.4) C C CALCULATE CONSTANTS FOR QUASSIAN DISPL. DISTRIBUTIONS C SIGMA=ESZ*(NDISP-4)/2.355 XG=-ESZ*(NDISP-1)/2.0 C C CALCULATE CUTOFF TIME FOR TWO COMPLETE CYCLES OF PRESCRIBED BOUNDARY DISP. C TCUT=4.*ACOS(-1.)/WC IF(T.GT.TCUT)AMP=0.0 C C LOOP THRU PRESCRIBED NODAL DISPL.S C DO 30 IDISP=1,NDISP C C RECOVER EQN NO. (IEQN) AND DIAGONAL ADDRESSES MAXA(IEQN) C FOR THE NODE WHERE A DISPLACEMENT IS PRESCRIBED C IEQN=ID(NOD(IDISP),IDIRN(IDISP)) II=MAXA(IEQN) C C WRITE(LW,31)IDISP,NOD(IDISP),IDIRN(IDISP),IEQN,AMP,II,A(II),IEQN C *,V(IEQN) C 31 FORMAT(/,5X,' IDISP=',I4,' NODE NO.=',I4,' DISPL. DIRECTION=', C *I3,/,6X,'EQN NO.=',I5,14X,'DISPL. PEAK AMPLITUDE=',E10.3,//, C *5X,'OLD DIAGONAL A(',I4,')=',E11.4, C *' OLD RIGHT-HAND-SIDE VECTOR V(',I4,')=',E11.4) C C CALCULATE PRESCRIBED BOUNDARY DISPLACEMENTS FOR THE CURRENT TIME STEP C CONST=AMP*EXP(-XG**2/(2.*SIGMA**2)) U(IEQN)=CONST*SIN(WC*T) XG=XG+ESZ C C ADD SPRING STIFFNESS AND SPRING FORCE TO GLOBAL ARRAYS C A(II)=A(II)+HUGEK V(IEQN)=V(IEQN)+HUGEK*U(IEQN) C C WRITE(LW,51)II,A(II),IEQN,V(IEQN) C 51 FORMAT(5X,'NEW DIAGONAL A(',I4,')=',E11.4, C *' NEW RIGHT-HAND-SIDE VECTOR V(',I4,')=',E11.4,/) C 30 CONTINUE C C IF THIS IS THE FIRST TIME STEP CALCULATE THE INITIAL ACCELERATIONS FROM C THE PRESCRIBED BOUNDARY DISPLACEMENTS OTHERWISE CONTINUE C IF(T.GT.0.0)GO TO 35 WRITE(10)(A(I),I=1,NCOFS) WRITE(11)(V(I),I=1,NEQ) CALL ACCEL(U,UDD,MAXA,A,V,MA,NEQ,NEQ1,NCOFS) REWIND 7 REWIND 10 REWIND 11 READ(7)(MA(I),I=1,NCOFS) READ(10)(A(I),I=1,NCOFS) READ(11)(V(I),I=1,NEQ) C C DO 41 IEQN=1,NEQ C 41 WRITE(LW,42)IEQN,U(IEQN),UD(IEQN),UDD(IEQN) C 42 FORMAT(1X,I6,3(1X,E11.4)) C C CALCULATE THE LOAD VECTOR INCLUDING THE INERTIAL TERMS OF NEWMARK APPROX. C FOR EACH EQUATION WHERE A DEGREE OF FREEDOM EXISTS C 35 DO 70 IEQN=1,NEQ II=MAXA(IEQN) C C INITIALIZE VARIABLES USED TO POSITION SUMMATION OF TERMS C AT LEFT SIDE OF EQN NO. (IEQN) C C N , NO. OF COMPONENTS TO THE LEFT AND INCLUDING DIAGONAL C IN, SUBSCRIPT OF FIRST U, UD, UDD ARRAYS IN EQN NO. (IEQN) C IM, SUBSCRIPT OF FIRST MASS COMPONENT IN EQN NO. (IEQN) C N=MAXA(IEQN+1)-II IN=IEQN-N+1 IM=MAXA(IEQN+1)-1 C C SUM UP PRODUCTS OF (U, UD, UDD) AND MASSES C FOR TERMS TO THE LEFT AND INCLUDING DIAGONAL C FORCE=0.0 C DO 40 K=1,N FORCE=FORCE+MA(IM)*(A0*U(IN)+A2*UD(IN)+A3*UDD(IN)) IM=IM-1 IN=IN+1 40 CONTINUE C C SUM UP REMAINING PRODUCT OF TERMS TO RIGHT OF DIAGONAL. C BUT FIRST RECALCULATE TRUNCATED BANDWIDTH (MB) TO THE RIGHT C OF THE DIAGONAL IF THE EQN NO. (IEQN) IS GREATER THAN (LIMIT) C IF(IEQN.GT.LIMIT)MB=MBAND-1-IEQN+LIMIT IF(MB.EQ.0)GO TO 65 C DO 60 IB=1,MB IM=MAXA(IN)+IB IF(IM.GE.MAXA(IN+1))GO TO 50 FORCE=FORCE+MA(IM)*(A0*U(IN)+A2*UD(IN)+A3*UDD(IN)) 50 IN=IN+1 60 CONTINUE C 65 V(IEQN)=V(IEQN)+FORCE C C WRITE(LW,61)IEQN,V(IEQN) C 61 FORMAT(1X,I6,1X,E11.4) C 70 CONTINUE C C END LOOP ON EQUATION NUMBERS C RETURN END C*********************************************************************** SUBROUTINE COLSOL(MAXA,A,V,NN,NWK,NNM,KKK) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * C * THIS SUBROUTINE IS CALLED TWICE * C * * C * THE FIRST TIME * C * * TO PERFORM L D L TRANSPOSE OPERATION * C * ON THE GLOBAL STIFFNESS MATRIX * C * * C * THE SECOND TIME * C * TO SOLVE THE SYSTEM OF EQUATIONS * C * * C * PROGRAM * C * TO SOLVE FINITE ELEMENT STATIC EQUILIBRIUM EQUATIONS IN * C * CORE, USING COMPACTED STORAGE AND COLUMN REDUCTION SCHEME * C * - - INPUT VARIABLES - - * C * A(NWK) = STIFFNESS MATRIX STORED IN COMPACTED FORM * C * V(NN) = RIGHT-HAND-SIDE LOAD VECTOR * C * MAXA(NNM) = VECTOR CONTAINING ADDRESSES OF DIAGONAL * C * ELEMENTS OF STIFFNESS MATRIX IN A * C * NN = NEQ = NUMBER OF EQUATIONS * C * NWK (NOT USED)= NUMBER OF ELEMENTS BELOW SKYLINE OF MATRIX * C * NNM (NOT USED)= NN + 1 * C * KKK = INPUT FLAG * C * EQ. 1 TRIANLGULARIZATION OF STIFFNESS MATRIX * C * EQ. 2 REDUCTION AND BACK-SUBSTITUTION OF LOAD VECTOR * C * IOUT = NUMBER OF OUTPUT DEVICE * C * - - OUTPUT - - * C * A(NWK) = D AND L - FACTORS OF STIFFNESS MATRIX * C * V(NN) = DISPLACEMENT VECTOR * C * * * * * * * * * * * * * * * * * * * * * * * * C C IMPLICIT REAL*8(A-H,0-Z) C C * * * * * * * * * * * * * * * * * * * * * * * * C * THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON * C * CDC EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM * C * OR UNIVAC MACHINES. ACTIVATE, DEACTIVATE OR ADJUST ABOVE * C * CARD FOR SINGLE OR DOUBLE PRECISION ARITHMETIC. * C * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/IOUT,LR,KEY C DIMENSION MAXA(NNM),A(NWK),V(NN) C C *** PERFORM L*D*L(T) FACTORIZATION OF STIFFNESS MATRIX C IF(KKK.EQ.1)THEN DO 140 N=1,NN KN=MAXA(N) KL=KN+1 KU=MAXA(N+1)-1 KH=KU-KL C IF(KH .GT.0) THEN K=N-KH IC=0 KLT=KU DO 80 J=1,KH IC=IC+1 KLT=KLT-1 KI=MAXA(K) ND=MAXA(K+1)-KI-1 C IF(ND.GT.0)THEN KK=IC C IF(ND.LT.IC)THEN KK=ND ENDIF C C=0.0 DO 70 L=1,KK C=C+A(KI+L)*A(KLT+L) 70 CONTINUE A(KLT)=A(KLT)-C ENDIF C K=K+1 80 CONTINUE C ENDIF IF(KH.GE.0) THEN K=N B=0.0 DO 100 KK=KL,KU K=K-1 KI=MAXA(K) C=A(KK)/A(KI) B=B+C*A(KK) A(KK)=C 100 CONTINUE A(KN)=A(KN)-B C ENDIF C IF(A(KN).LE.0)THEN WRITE (IOUT,2000) N,A(KN) STOP ENDIF C 140 CONTINUE RETURN C ENDIF C C *** REDUCE RIGHT-HAND-SIDE LOAD VECTOR C DO 180 N=1,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 C IF((KU-KL).GE.0)THEN K=N C=0.0 DO 170 KK=KL,KU K=K-1 C=C+A(KK)*V(K) 170 CONTINUE V(N)=V(N)-C ENDIF C 180 CONTINUE C C *** BACK SUBSTITUTE C DO 200 N=1,NN K=MAXA(N) V(N)=V(N)/A(K) 200 CONTINUE C IF (NN.EQ.1) RETURN C N=NN DO 230 L=2,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 C IF((KU-KL).GE.0)THEN K=N DO 220 KK=KL,KU K=K-1 V(K)=V(K)-A(KK)*V(N) 220 CONTINUE ENDIF C N=N-1 230 CONTINUE RETURN C 2000 FORMAT (//48H STOP - STIFFNESS MATRIX NOT POSITIVE DEFINITE ,/ */ 32H NONPOSITIVE PIVOT FOR EQUATION ,I4// * 10H PIVOT = ,E20.12) END C*********************************************************************** SUBROUTINE DISP(ID,D,V,TUDD,U,UD,UDD,NDS,NBN,NEQ,T,ITINC * ,ICNTG,ICNTB) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * RECOVERS NODAL DISPL.S STORED IN V(NEQ) BY COLSOL. * C * AND ASSIGNS NODE DISPL.S TO ARRAY D(NDS,2). * C * AND CALCULATES NEW VEL.S AND ACCEL.S FOR EACH TIME STEP * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY COMMON/RDK3/A0,A2,A3,A6,A7,WC,RHO,AMP COMMON/RDK4/IBN(200) C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),D(NDS,2),V(NEQ),TUDD(NEQ) DIMENSION U(NEQ),UD(NEQ),UDD(NEQ) C C RECOVER NODAL DISPLACEMENTS C IEQN=1 DO 10 I=1,NDS DO 10 J=1,2 D(I,J)=0.0 IF(ID(I,J).GT.0)THEN D(I,J)=V(IEQN) IEQN=IEQN+1 ENDIF 10 CONTINUE C C UPDATE NEW VELOCITIES AND ACCELERATIONS FOR EACH TIME STEP C DO 20 IEQN=1,NEQ TUDD(IEQN)=A0*(V(IEQN)-U(IEQN))-A2*UD(IEQN)-A3*UDD(IEQN) UD(IEQN)=UD(IEQN)+A6*UDD(IEQN)+A7*TUDD(IEQN) UDD(IEQN)=TUDD(IEQN) U(IEQN)=V(IEQN) 20 CONTINUE C C WRITE WHOLE FIELD DISPL.S ONTO TAPE4 FOR POST PROCESSING C IF(ICNTG.EQ.1)THEN WRITE(4,30)T 30 FORMAT(' TIME=',E11.4) DO 40 IND=1,NDS 40 WRITE(4,50)IND,(D(IND,J),J=1,2) 50 FORMAT(1X,I6,2(1X,E11.4)) ENDFILE 4 ENDIF C C WRITE BOUNDARY DISPLACEMENTS ONTO TAPE12 FOR POST PROCESSING C IF(ICNTB.EQ.1)THEN WRITE(12,30)T DO 60 I=1,NBN IND=IBN(I) 60 WRITE(12,50)IND,(D(IND,J),J=1,2) ENDFILE 12 ENDIF C RETURN END C************************************************************************ SUBROUTINE ACCEL(U,UDD,MAXA,A,V,MA,NEQ,NEQ1,NCOFS) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * CALLED BY SUBROUTINE PDVA * C * THIS SUBROUTINE CALCULATES THE INITIAL ACCELERATIONS FROM THE * C * PRESCRIBED DISPLACEMENTS AT THE BOUNDARY FOR THE FIRST TIME * C * TIME STEP ONLY * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C REAL MA COMMON/RDK1/LW,LR,KEY C DIMENSION MAXA(NEQ1),A(NCOFS),V(NEQ),MA(NCOFS) DIMENSION U(NEQ),UDD(NEQ) C C CLEAR LOAD VECTOR C DO 10 I=1,NEQ V(I)=0.0 10 CONTINUE C C READ THE GLOBAL STIFFNESS ARRAY FROM TAPE3 C REWIND 3 READ(3)(A(I),I=1,NCOFS) C C CALCULATE LOAD VECTOR FOR THE PRESCRIBED BOUNDARY DISPLACEMENTS C DO 70 IEQN=1,NEQ II=MAXA(IEQN) C C INITIALIZE VARIABLES USED TO POSITION SUMMATION OF TERMS C AT LEFT SIDE OF EQN NO. (IEQN) C C N , NO. OF COMPONENTS TO THE LEFT AND INCLUDING DIAGONAL C IN, SUBSCRIPT OF FIRST U ARRAY IN EQN NO. (IEQN) C IM, SUBSCRIPT OF FIRST STIFFNESS COMPONENT IN EQN NO. (IEQN) C N=MAXA(IEQN+1)-II IN=IEQN-N+1 IM=MAXA(IEQN+1)-1 C C SUM UP PRODUCTS OF DISPLACEMENTS (U) AND STIFFNESS (A) C FOR TERMS TO THE LEFT AND INCLUDING DIAGONAL C FORCE=0.0 C DO 40 K=1,N FORCE=FORCE+A(IM)*U(IN) IM=IM-1 IN=IN+1 40 CONTINUE C C SUM UP REMAINING PRODUCT OF TERMS TO RIGHT OF DIAGONAL. C BUT FIRST RECALCULATE TRUNCATED BANDWIDTH (MB) TO THE RIGHT C OF THE DIAGONAL IF THE EQN NO. (IEQN) IS GREATER THAN (LIMIT) C IF(IEQN.GT.LIMIT)MB=MBAND-1-IEQN+LIMIT IF(MB.EQ.0)GO TO 65 C DO 60 IB=1,MB IM=MAXA(IN)+IB IF(IM.GE.MAXA(IN+1))GO TO 50 FORCE=FORCE+A(IM)*U(IN) 50 IN=IN+1 60 CONTINUE C 65 V(IEQN)=V(IEQN)+FORCE C 70 CONTINUE C C END LOOP ON EQUATION NUMBERS C C CALCULATE THE INITIAL ACCELERATIONS C CALL COLSOL(MAXA,MA,V,NEQ,NCOFS,NEQ1,1) CALL COLSOL(MAXA,MA,V,NEQ,NCOFS,NEQ1,2) C C TRANSFER SOLUTION TO THE UDD ARRAY C DO 80 IEQN=1,NEQ UDD(IEQN)=V(IEQN) 80 CONTINUE C RETURN END