PROGRAM wave1d_40x40 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 F ORCE F INITE E LEMENT W AVE A NISOTROPIC * 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,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 COMMON SHARE(13000) MTOT=13000 C C ASSIGN TAPE NO.S 5 AND 6 TO INPUT AND OUTPUT FILES RESPECTIVELY C LR=5 LW=6 OPEN(LR,FILE='wave1d_4x40.DAT',STATUS='OLD',ERR=888) OPEN(LW,FILE='wave1d_4x40.out',STATUS='UNKNOWN',ERR=888) C C ASSIGN TAPE NO.s 10, 11 and 12 to PLTWD, QF3D and QF103D RESPECTIVELY C OPEN(10,FILE='pltwave1d_4x40.dat',STATUS='UNKNOWN',ERR=888) OPEN(11,FILE='node003d.dat',STATUS='UNKNOWN',ERR=888) OPEN(12,FILE='node103d.dat',STATUS='UNKNOWN',ERR=888) 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(4,FILE='temp4',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') 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,NLOAD,NDISP,KEY c 20 FORMAT(/,5I6,/) READ(LR,*)NE,NDS,NLOAD,NDISP,KEY C WRITE(LW,22)NE,NDS,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' NUMBER OF NODE LOADS .......................',I6,/, 4' NO. OF PRESCRIBED NODE DISPLACEMENTS .......',I6,/, 6' OUTPUT FILE PRINT KEY ......................',I6) 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,42)THETA READ(LR,*)THETA c 40 FORMAT(3(6X,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,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 WRITE(LW,64)ITYPE 64 FORMAT(/,10X,'CONSTRAINT TYPE = ',I1, */,23X,' ( = 1, PLANE STRAIN)', */,23X,' ( = 2, PLANE STRESS)',//) 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,FT,IPRNT READ(LR,*)DT READ(LR,*)ALT READ(LR,*)DELT READ(LR,*)FT READ(LR,*)IPRNT c 68 FORMAT(4(1X,E11.4,/),1X,I11,/) WRITE(LW,69)DT,ALT,DELT,FT,IPRNT 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, */,' FINAL TIME, FT=',E11.4, */,' PRINT INTERVAL, IPRNT=',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,AMP READ(LR,*)WF READ(LR,*)WS READ(LR,*)RHO READ(LR,*)NWL READ(LR,*)AMP c 71 FORMAT(3(1X,E11.4,/),1X,I11,/,1X,E11.4,/) 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 WRITE(LW,72)WF,WS,RHO,NWL,AMP,WC,WP,WL,ESZ 72 FORMAT(' PERIODIC WAVE PROPERTIS AT BOUNDARY',/ * ,' WAVE FREQUENCY, WF=',E11.4,/ * ,' WAVE SPEED, WS=',E11.4,/ * ,' MATERIAL DENSITY, RHO=',E11.4,/ * ,' ELEMENTS/WAVE LENGTH, NWL=',I3,/ * ,' PEAK FORCE AMPLT., AMP=',E11.4,/ * ,' WAVE TIME-CONSTANT, WC=',E11.4,/ * ,' WAVE PERIOD, WP=',E11.4,/ * ,' WAVE LENGTH, WF=',E11.4,/ * ,' ELEMENT SIZE, ESZ=',E11.4,//) C C WRITE INFO. FOR PLOTTING WAVE PROPAGATION ON TAPE4 C XSCALP=0.25 YSCALP=0.25 XORIG=0.25 YORIG=2.50 WRITE(10,11)(TITLE(I),I=1,20) WRITE(10,24)NDS,NE,XSCALP,YSCALP,XORIG,YORIG 24 FORMAT(2I5,4F8.3) 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 25 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,ESZ) C C WRITE INFO. FOR PLOTTING WAVE PROPAGATION ON TAPE4 C SCALD=1.0/ESZ NINC=FT/DT/IPRNT+1 WRITE(10,26)NINC,SCALD 26 FORMAT(I2,E10.3) 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 C ARRAYS USED IN SUBROUTINE PLOAD 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 NTOT1=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 73 FORMAT(' TOTAL SIZE OF ARRAYS STORED IN COMMON',/, 1' AVAILBLE(',I6,') REQUIRED BY SUBROUTINE PLOAD (',I6,')',/) WRITE(LW,999) C IF(NTOT1.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.0 ICNT=0 C C CALCULATE P RESCRIDED L OADS FOR EACH TIME STEP C ID(NDS,2) MAXA(NEQ1) A(NCOFS) V(NEQ) MA(NCOFS) 100 CALL PLOAD(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),SHARE(N5), *SHARE(N6),SHARE(N7),SHARE(N8),NDS,NEQ1,NCOFS,NEQ,MBAND,NLOAD,T,DT) C U(NEQ) UD(NEQ) UDD(NEQ) 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,NEQ,T,ICNT) C U(NEQ) UD(NEQ) UDD(NEQ) C C I F R E Q U E S T E D C A L C U L A T E S T R E S S E S C C D(NDS,2) IF((KEY.EQ.6).OR.(KEY.EQ.7))CALL STR(SHARE(N3),NE,NDS,DD) C C RECOVER CPU TIME AND WRITE TO OUTPUT FILE C call time(tt) WRITE(LW,333)T,tt,ICNT 333 FORMAT(' T=',E11.4,' CPU TIME=',A8,' ICNT=',I2) C C I N C R E M E N T T I M E S T E P C T=T+DT ICNT=ICNT+1 IF(ICNT.EQ.IPRNT)ICNT=0 IF(T.LE.FT)GO TO 100 C C E N D O F T I M E L O O P C 888 STOP END C*********************************************************************** SUBROUTINE NEDATA(ID,H,Z,NDS,NE,ESZ) 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,THETA C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),H(NDS),Z(NDS) 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 C FOR GRAPHICAL POST PROCESSING C WRITE(10,21)I,ID(I,1),ID(I,2),H(I),Z(I) 21 FORMAT(I6,2I3,2F12.5) 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 FOR GRAPHICAL POST PROCESSING C WRITE(10,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 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 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 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=',I4) 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=',I6,/) C IF((KEY.EQ.4).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 + A0*MA(I) * 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,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 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 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 MASS MATRICES ON OUTPUT FILE C 2. RIGHT-HAND-SIDE LOAD VECTOR ON OUTPUT FILE C 3. GLOBAL MASS MATRIX 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),MA(1) 75 FORMAT(1H1,//,' GLOBAL STIFFNESS ARRAY', *//,' COLUMN NO. 1', */,' 1 RST DIAGONAL POSITION', */,7X,' A( 1)=',E11.4,7X,' MA( 1)=',E11.4,/) 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),I,MA(I) 76 FORMAT(8X,'A(',I4,')=',E14.7,8X,'MA(',I4,')=',E14.7) GO TO 84 77 IF(I.NE.NCOFS)WRITE(LW,78)IDIAG,I,A(I),I,MA(I),IDIAG 78 FORMAT(/,' COLUMN NO.',I4,/,8X,'A(',I4,')=',E14.7,8X,'MA(',I4, *')=',E14.7,' DIAGONAL NO.',I4) C IF(I.EQ.NCOFS)WRITE(LW,79)I,A(I),I,MA(I) 79 FORMAT(//,' LAST STORED VALUES',//, *8X,'A(',I4,')=',E14.7,8X,'MA(',I4,')=',E14.7) 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 C ADD THE GLOBAL MASS MATRIX TO THE STIFFNESS MATRIX FOR DYNAMIC SOLUTION C 999 DO 110 I=1,NCOFS A(I)=A(I)+A0*MA(I) 110 CONTINUE C C STORE GLOBAL STIFFNESSES AND LOADS 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 K=1,2 DO 40 J=1,8 H(K,J)=0.0 DO 40 I=1,3 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 PLOAD(ID,MAXA,A,V,MA,U,UD,UDD,NDS,NEQ1,NCOFS,NEQ,MBAND, * NLOAD,T,DT) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M PLOAD * C * 1. ADDS THE PRESCRIBED BOUNDARY FORCES TO THE LOAD VECTOR. * C * FOR THE FIRST TIME STEP THE ACCELERATIONS ARE CALCULATED * C * BY USING THE COLSOL EQUATION SOLVER. * 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 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) DIMENSION NOD(20),IDIRN(20),FLOAD(20) 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 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 LOAD IS PRESCRIBED C 2. IDIRN, DIRECTION OF LOAD C 3. FLOAD, PERCENTAGE OF LOAD AT NODE NOD(ILOAD) C (NOTE: THAT THE AMPLITUDE OF ALL NODAL LOADS WAS PREVIOUSLY C READ FROM THE INPUT FILE AS THE VARIABLE AMP) C IF(T.GT.0.0)GO TO 16 DO 10 ILOAD=1,NLOAD 10 READ(LR,*)NOD(ILOAD),IDIRN(ILOAD),FLOAD(ILOAD) c 10 READ(LR,15)NOD(ILOAD),IDIRN(ILOAD),FLOAD(ILOAD) c 15 FORMAT(2I6,F8.2) C 16 continue c 16 WRITE(LW,17)T c 17 FORMAT(1X,80(1H*),/,' TIME=',E11.4) C C RECOVER EQN NO. (IEQN) C FOR THE NODE WHERE A BOUNDARY LOAD IS PRESCRIBED C DO 30 ILOAD=1,NLOAD IEQN=ID(NOD(ILOAD),IDIRN(ILOAD)) C C WRITE(LW,31)ILOAD,NOD(ILOAD),IDIRN(ILOAD),IEQN,AMP,IEQN,V(IEQN) C 31 FORMAT(/,5X,' ILOAD=',I4,' NODE NO.=',I4,' LOAD DIRECTION=', C *I3,/,6X,'EQN NO.=',I5,14X,'FORCE PEAK AMPLITUDE=',E10.3,//, C *' OLD RIGHT-HAND-SIDE VECTOR V(',I4,')=',E11.4) C C ADD FORCE TO GLOBAL V(I) ARRAY C V(IEQN)=V(IEQN)+FLOAD(ILOAD)*BFORCE(T) C C WRITE(LW,51)ILOAD,FLOAD(ILOAD),BFORCE(T),IEQN,V(IEQN) C 51 FORMAT(5X,' FLOAD(',I2,')=',E10.3,' BFORCE=',E11.4,/ C *,' NEW RIGHT-HAND-SIDE VECTOR V(',I4,')=',E11.4,/) C 30 CONTINUE C C FOR THE FIRST TIME STEP CALCULATE THE INITIAL ACCELERATIONS FROM C THE INITIAL LOADS WHERE THE INITIAL DISPLACEMENTS ARE ZERO. C IF(T.EQ.0.0)RETURN IF(T.GT.DT)GO TO 36 WRITE(4)(V(I),I=1,NEQ) CALL COLSOL(MAXA,MA,V,NEQ,NCOFS,NEQ1,1) CALL COLSOL(MAXA,MA,V,NEQ,NCOFS,NEQ1,2) DO 35 IEQN=1,NEQ UDD(IEQN)=V(IEQN) 35 CONTINUE C WRITE(LW,34)(IEQN,UDD(IEQN),IEQN=1,NEQ) C 34 FORMAT(' UDD(',I3,')=',E10.4,' UDD(',I3,')=',E10.4 C * ,' UDD(',I3,')=',E10.4,' UDD(',I3,')=',E10.4) REWIND 7 REWIND 4 READ(7)(MA(I),I=1,NCOFS) READ(4)(V(I),I=1,NEQ) C C CALCULATE THE LOAD VECTOR INCLUDING THE INERTIAL TERMS OF NEWMARK APPROX. C FOR EACH EQUATION WHERE A DEGREE OF FREEDOM EXISTS C 36 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 VEC=A0*U(IN)+A2*UD(IN)+A3*UDD(IN) FORCE=FORCE+MA(IM)*VEC C WRITE(LW,22)IN,VEC,IM,MA(IM),FORCE C 22 FORMAT(' VECTOR(',I4,')=',E11.4,' MASS(',I4,')=',E11.4, C *' FORCE=',E11.4) 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 VEC=A0*U(IN)+A2*UD(IN)+A3*UDD(IN) FORCE=FORCE+MA(IM)*VEC C WRITE(LW,22)IN,VEC,IM,MA(IM),FORCE 50 IN=IN+1 60 CONTINUE C 65 V(IEQN)=V(IEQN)+FORCE C C WRITE(LW,61)IEQN,V(IEQN) C 61 FORMAT(' V(',I3,')=',E11.4,/,1X,80(1H*)) 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-2) 40,150,150 40 DO 140 N=1,NN KN=MAXA(N) KL=KN+1 KU=MAXA(N+1)-1 KH=KU-KL IF (KH) 110,90,50 50 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 IF (ND) 80,80,60 60 KK=IC IF (ND.LT.IC)KK=ND C=0.0 DO 70 L=1,KK 70 C=C+A(KI+L)*A(KLT+L) A(KLT)=A(KLT)-C 80 K=K+1 90 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) 100 A(KK)=C A(KN)=A(KN)-B 110 IF (A(KN)) 120,120,140 120 WRITE (IOUT,2000) N,A(KN) STOP 140 CONTINUE RETURN C C *** REDUCE RIGHT-HAND-SIDE LOAD VECTOR C 150 DO 180 N=1,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 IF (KU-KL) 180,160,160 160 K=N C=0.0 DO 170 KK=KL,KU K=K-1 170 C=C+A(KK)*V(K) V(N)=V(N)-C 180 CONTINUE C C *** BACK SUBSTITUTE C DO 200 N=1,NN K=MAXA(N) 200 V(N)=V(N)/A(K) IF (NN.EQ.1) RETURN N=NN DO 230 L=2,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 IF (KU-KL) 230,210,210 210 K=N DO 220 KK=KL,KU K=K-1 220 V(K)=V(K)-A(KK)*V(N) 230 N=N-1 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,NEQ,T,ICNT) 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 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).EQ.0)GO TO 10 D(I,J)=V(IEQN) IEQN=IEQN+1 10 CONTINUE WRITE(11,11)T,D(3,1),D(3,2) WRITE(12,11)T,D(103,1),D(103,2) 11 FORMAT(' TIME=',E10.3,' X-DISPL.=',E11.4,' Y-DISPL.=',E11.4) 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 WRITE(10,30)T C WRITE(LW,34)(IEQN,UDD(IEQN),IEQN=1,NEQ) C 34 FORMAT(' UDD(',I3,')=',E10.4,' UDD(',I3,')=',E10.4 C * ,' UDD(',I3,')=',E10.4,' UDD(',I3,')=',E10.4) C WRITE(LW,35)(IEQN,UD(IEQN),IEQN=1,NEQ) C 35 FORMAT(' UD(',I3,')=',E10.4,' UD(',I3,')=',E10.4 C * ,' UD(',I3,')=',E10.4,' UD(',I3,')=',E10.4) C WRITE(LW,36)(IEQN,U(IEQN),IEQN=1,NEQ) C 36 FORMAT(' U(',I3,')=',E10.4,' U(',I3,')=',E10.4 C * ,' U(',I3,')=',E10.4,' U(',I3,')=',E10.4) C C WRITE DISPLACEMENTS ONTO TAPE4 FOR GRAPHICAL POST PROCESSING C IF(ICNT.GT.0)GO TO 60 WRITE(10,30)T 30 FORMAT(' TIME=',E11.4) DO 40 IND=1,NDS 40 WRITE(10,50)IND,(D(IND,J),J=1,2) 50 FORMAT(1X,I6,2(1X,E11.4)) ENDFILE 4 C 60 RETURN END C*********************************************************************** SUBROUTINE STR(D,NE,NDS,DD) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATES STR ESSES FROM STRAINS FROM DISPLACEMENTS * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION D(NDS,2) DIMENSION B(3,8),DD(3,3),IND(4),E(3),S(3) C C R E W I N D S T O R E D E L E M E N T D A T A C C TAPE 1 ND1,ND2,ND3,ND4,THIC REWIND 1 C C TAPE 2 NELS,LXS,LXY,((B(I,J),I=1,3),J=1,8) REWIND 2 C C WRITE HEADER FOR STRESSES ON OUTPUT FILE CONDITIONAL ON KEY C IF((KEY.EQ.6).OR.(KEY.EQ.7))WRITE(LW,5) 5 FORMAT(1H1,//,25X,'** STRESSES **',/, * ' ELE. LX LY',7X,'SH',13X,'SZ',13X,'SHZ') C C L O O P T H R U E L E M E N T S C DO 70 NEL=1,NE C C READ ELEMENT DATA READ(1)(IND(I),I=1,4),DUM C C CHOSE NO. OF QUASS POINTS C NINT=2 C C LOOP THRU (NINT)X(NINT) QUASS POINTS C DO 30 LX=1,NINT DO 30 LY=1,NINT C READ(2)NELS,LXS,LYS,((B(I,J),I=1,3),J=1,8) C C COMPARE DO LOOP PARAMETERS WITH DATA ON TAPE C IF(NEL.EQ.NELS)GO TO 7 GO TO 9 7 IF(LX.EQ.LXS)GO TO 8 GO TO 9 8 IF(LY.EQ.LYS)GO TO 10 9 WRITE(LW,15)NEL,NELS,LX,LXS,LY,LYS 15 FORMAT(//,' *** ERROR *** DATA ON TAPE-2 AND DO LOOPS',/, * ' DOES NOT COMPARE IDENTICALLY',/, * ' NEL=',I6,' NELS=',I6,/, * ' LX=',I2,' LXS=',I2,/, * ' LY=',I2,' LYS=',I2) STOP C C C CALCULATE STRAINS FROM DISPLACEMENTS C 10 E(1)=B(1,1)*D(IND(1),1)+B(1,2)*D(IND(1),2) * +B(1,3)*D(IND(2),1)+B(1,4)*D(IND(2),2) * +B(1,5)*D(IND(3),1)+B(1,6)*D(IND(3),2) * +B(1,7)*D(IND(4),1)+B(1,8)*D(IND(4),2) C E(2)=B(2,1)*D(IND(1),1)+B(2,2)*D(IND(1),2) * +B(2,3)*D(IND(2),1)+B(2,4)*D(IND(2),2) * +B(2,5)*D(IND(3),1)+B(2,6)*D(IND(3),2) * +B(2,7)*D(IND(4),1)+B(2,8)*D(IND(4),2) C E(3)=B(3,1)*D(IND(1),1)+B(3,2)*D(IND(1),2) * +B(3,3)*D(IND(2),1)+B(3,4)*D(IND(2),2) * +B(3,5)*D(IND(3),1)+B(3,6)*D(IND(3),2) * +B(3,7)*D(IND(4),1)+B(3,8)*D(IND(4),2) C C CALCULATE STRESSES FROM STRAINS C DO 20 J=1,3 S(J)=0.0 DO 20 I=1,3 20 S(J)=S(J)+DD(J,I)*E(I) C C WRITE STRESSES ON OUTPUT FILE C WRITE(LW,25)NEL,LX,LY,S(1),S(2),S(3) 25 FORMAT(2X,I5,2I3,3(1X,E14.7)) C 30 CONTINUE C C E N D Q U A S S P O I N T L O O P C 70 CONTINUE C C E N D E L E M E N T L O O P C RETURN END C******************************************************************************* FUNCTION BFORCE(T) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * FUNCTION BFORCE * C * CALCULATES THE MAGNITUDE OF THE FORCE PRESCRIBED AT THE BOUNDARY. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK3/A0,A2,A3,A6,A7,WC,RHO,AMP C C DEFINE THE MAGNITUDE OF THE BOUNDARY FORCE C IF(T.EQ.0.0)BFORCE=0.0 IF(T.GT.0.0)BFORCE=AMP C RETURN END