program q8fema c PROGRAM QFEMA(INPUT,OUTPUT,TAPE1,TAPE3,TAPE4,TAPE5,TAPE6 c * ,TAPE7,TAPE8,TAPE9) 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: ( QFEMA ) * C * Q UADRILATERAL F INITE E LEMENT M ODEL 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 C COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E1,E2,E3,G12,G13,G23,NU12,NU13,NU23,AL1,AL2,AL3,THETA C DIMENSION TITLE(20),D(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(35000) MTOT=35000 c c openfiles for input and output c open(5,FILE='q8fema.dat',STATUS='OLD',ERR=888) OPEN(6,FILE='q8fema.out',STATUS='UNKNOWN',ERR=888) C C ASSIGN TAPE NO.S 5 AND 6 TO INPUT AND OUTPUT FILES RESPECTIVELY C LR=5 LW=6 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 G R I D P A R A M E T E R S C READ(LR,20)NE,NDS,NLOAD,NDISP,KEY 20 FORMAT(/,5I6) READ(LR,21)SMH,SMZ 21 FORMAT(2F12.0,/) C WRITE(LW,22)NE,NDS,NLOAD,NDISP,SMH,SMZ,KEY 22 FORMAT(//,5X,'G R I D P A R A M E T E R S',/, 1' NUMBER OF ELEMENTS .........................',I4,/, 2' NUMBER OF NODES ............................',I4,/, 3' NUMBER OF NODE LOADS .......................',I4,/, 4' NO. OF PRESCRIBED NODE DISPLACEMENTS .......',I4,/, 5' SCALE FOR HORIZONTAL COORD.S H(I) .........',F12.5,/, 6' SCALE FOR VERTICAL COORD.S Z(I) .........',F12.5,/, 7' OUTPUT FILE PRINT KEY ......................',I4) 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 READ(LR,40)E1,E2,E3 READ(LR,40)G12,G13,G23 READ(LR,40)NU12,NU13,NU23 READ(LR,41)AL1,AL2,AL3 READ(LR,42)THETA 40 FORMAT(3(6X,E10.3)) 41 FORMAT(3(8X,E10.3)) 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,/, * 14X,'G13=',E10.3,/, * 14X,'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=',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 READ(LR,61)ITYPE 61 FORMAT(I6) READ(LR,61)LTYPE READ(LR,62)TLOAD 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 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,SMH,SMZ * ,LTYPE,TLOAD) 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 ELEMENT ARRAY SIZE, NDOFE (DEGREES OF FREEDOM PER ELEMENT) 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 ELEMNT 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 NTOT1=2*NDS+NEQ1+NCOFS+NEQ C C ARRAYS USED IN SUBROUTINE STR AND SHARED BY MAIN PROG. COMMON C C ARRAY NAMES COMMON ADDRESSES C ID( NDS,2 ) = SHARE(N1) N1=1 SPACE ALLOCATED FOR ID ARRAY C PRIOR TO SUBPROG. STR C U( NDS,2 ) = SHARE(N2) N2=N1+NDS*2 C EH(NE,8,8) = SHARE(N3) N3=N2+NDS*2 C EZ(NE,8,8) = SHARE(N4) N4=N3+NE*8,8 C EHZ(NE,8,8) = SHARE(N5) N5=N4+NE*8,8 C NTOT2=4*NDS+3*NE*64 C C CHOOSE LARGEST SHARED ARRAY SIZE C NTOT=NTOT1 IF(NTOT2.GT.NTOT1)NTOT=NTOT2 C WRITE(LW,99) WRITE(LW,999) WRITE(LW,70)MTOT,NTOT1,NTOT2 70 FORMAT(' TOTAL SIZE OF ARRAYS STORED IN COMMON',/, 1' AVAILBLE(',I6,') REQUIRED BY SUBROUTINE ELEMNT (',I6,')',/, 233X,'SUBROUTINE STR (',I6,')') WRITE(LW,999) C IF(NTOT.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 A R R A Y S C CONNECTIVITY STIFFNESS FORCE C LM( NDOFE ) AK( NDOFE,NDOFE ) F( NDOFE ) C C ---------- AND ---------- 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 C A( NCOFS ) V( NEQ ) C C ADDRESSES N3=N2+NEQ1 N4=N3+NCOFS C C ID(NDS,2) MAXA(NEQ1) A(NCOFS) V(NEQ) CALL QUAD8(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),NDS,NEQ1,NCOFS * ,NEQ,ITYPE,NE,D) C C I F A V A I L A B L E C S U P E R P O S E N O D E F O R C E S C C ID(NDS,2) V(NEQ) IF(NLOAD.GT.0)CALL NDLOAD(SHARE(N1),SHARE(N4),NDS,NEQ,NLOAD) C C I F A V A I L A B L E I N C O R P O R A T E C P R E S C R I B E D N O D E D I S P L A C E M E N T S C C ID(NDS,2) MAXA(NEQ1) A(NCOFS) V(NEQ) IF(NDISP.GT.0)CALL NDDISP(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4), * NDS,NEQ1,NCOFS,NEQ,NDISP) 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 ID(NDS,2) MAXA(NEQ1) A(NCOFS) V(NEQ) IF(NDISP.GT.0)CALL NDFORC(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4), * NDS,NEQ1,NCOFS,NEQ,NDISP,MBAND) C C R E C O V E R D I S P L A C E M E N T S U(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 C ID(NDS,2) U(NDS,2) V(NEQ) CALL DISP(SHARE(N1),SHARE(N2),SHARE(N4),NDS,NEQ) C C C A L C U L A T E S T R E S S E S A N D S T R A I N S C C ADDRESSES N3=N2+2*NDS N4=N3+NE*64 N5=N4+NE*64 C C SYMMETRIC STRAIN TENSOR COMPONENTS C W.R.T. LAMINATE AXIS (H,Z) C ( EH , EZ , EHZ ) C DO 333 J=1,3 333 WRITE(6,334)(D(J,I),I=1,3) 334 FORMAT(3E10.3) C U(NDS,2),EH(NE,8,8),EZ(NE,8,8),EHZ(NE,8,8) CALL STR(SHARE(N2),SHARE(N3),SHARE(N4),SHARE(N5), * NE,NDS,D) C 888 STOP END C*********************************************************************** SUBROUTINE NEDATA(ID,H,Z,NDS,NE,SMH,SMZ,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 INTEGER NU12,NU13,NU23 C COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E1,E2,E3,NU12,NU13,NU23,G12,G13,G23,AL1,AL2,AL3,THETA 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 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 5TH 6TH 7TH 8TH 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,ND5,ND6,ND7,ND8,THIC 60 FORMAT(6X,8I6,F12.0) C C STORE NODES ASSOCIATED WITH EACH ELEMENT ON TAPE 1 C WRITE(1)ND1,ND2,ND3,ND4,ND5,ND6,ND7,ND8,THIC,TS1,TS2,TS3 C C NODE CO-ORDINATES ASSOCIATED WITH EACH ELEMENT RESCALED C H1=H(ND1)*SMH H2=H(ND2)*SMH H3=H(ND3)*SMH H4=H(ND4)*SMH H5=H(ND5)*SMH H6=H(ND6)*SMH H7=H(ND7)*SMH H8=H(ND8)*SMH Z1=Z(ND1)*SMZ Z2=Z(ND2)*SMZ Z3=Z(ND3)*SMZ Z4=Z(ND4)*SMZ Z5=Z(ND5)*SMZ Z6=Z(ND6)*SMZ Z7=Z(ND7)*SMZ Z8=Z(ND8)*SMZ C C STORE ELEMENT NODE POINT COORD.S ON TAPE 3 C WRITE(3)H1,H2,H3,H4,H5,H6,H7,H8,Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8 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 * ,ND5,ND6,ND7,ND8,THIC 70 FORMAT(9(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(16) C C CLEAR LM AND MAXA ARRAYS C DO 5 IEQN=1,NEQ1 5 MAXA(IEQN)=0 DO 10 I=1,16 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,ND5,ND6,ND7,ND8) 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,ND5,ND6,ND7,ND8,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(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 5TH NODE)=ID(ND5,1)=LM(9) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 5TH NODE)=ID(ND5,2)=LM(10) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 6TH NODE)=ID(ND6,1)=LM(11) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 6TH NODE)=ID(ND6,2)=LM(12) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 7TH NODE)=ID(ND7,1)=LM(13) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 7TH NODE)=ID(ND7,2)=LM(14) C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 8TH NODE)=ID(ND8,1)=LM(15) C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 8TH NODE)=ID(ND8,2)=LM(16) 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) LM(9)=ID(ND5,1) LM(10)=ID(ND5,2) LM(11)=ID(ND6,1) LM(12)=ID(ND6,2) LM(13)=ID(ND7,1) LM(14)=ID(ND7,2) LM(15)=ID(ND8,1) LM(16)=ID(ND8,2) C C CHOOSE NO. OF DEGREES OF FREEDOM PER ELEMENT, NDOFE C NDOFE=16 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.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 QUAD8(ID,MAXA,A,V,NDS,NEQ1,NCOFS,NEQ,ITYPE,NE,D) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M - QUAD8 * C * * C * CALCULATES THE STIFFNESS MATRICES FOR EIGHT NODED QUADRILATERIAL * C * 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 STIFFNESS * 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. 8) * C * THIC = THICKNESS OF ELEMENT * C * E = YOUNG'S MODULUS * C * NU = POISSON'S RATIO * C * * C * - - CALCULATES - - * C * AK(NDOFE,NDOFE) = ELEMENT STIFFNESS 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 * V(NEQ) = GLOBAL LOAD VECTOR * C * D(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 COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E1,E2,E3,G12,G13,G23,NU12,NU13,NU23,AL1,AL2,AL3,THETA DIMENSION ID(NDS,2),MAXA(NEQ1),A(NCOFS),V(NEQ) DIMENSION XG(8,8),WGT(8,8),D(3,3),C(2,8),B(3,16),DB(3) DIMENSION AK(16,16),F(16),LM(16) C C ZERO ALL ELEMENTS OF XG NAD WXG ARRAYS TO 15 SIGNIFICANT FIGURES C DO 5 I=1,8 DO 5 J=1,8 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) XG(1,5)=-0.906179845938664 XG(2,5)=-0.538469310105683 XG(4,5)=-XG(2,5) XG(5,5)=-XG(1,5) XG(1,6)=-0.932469514203152 XG(2,6)=-0.661209386466265 XG(3,6)=-0.238619186083197 XG(4,6)=-XG(3,6) XG(5,6)=-XG(2,6) XG(6,6)=-XG(1,6) XG(1,7)=-0.949107912342759 XG(2,7)=-0.741531185599394 XG(3,7)=-0.405845151377397 XG(5,7)=-XG(3,7) XG(6,7)=-XG(2,7) XG(7,7)=-XG(1,7) XG(1,8)=-0.960289856497536 XG(2,8)=-0.796666477413627 XG(3,8)=-0.525532409916329 XG(4,8)=-0.183434642495650 XG(5,8)=-XG(4,8) XG(6,8)=-XG(3,8) XG(7,8)=-XG(2,8) XG(8,8)=-XG(1,8) 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) WGT(1,5)=0.236926885056189 WGT(2,5)=0.478628670499366 WGT(3,5)=0.568888888888889 WGT(4,5)=WGT(2,5) WGT(5,5)=WGT(1,5) WGT(1,6)=0.171324492379170 WGT(2,6)=0.360761573048139 WGT(3,6)=0.467913934572691 WGT(4,6)=WGT(3,6) WGT(5,6)=WGT(2,6) WGT(6,6)=WGT(1,6) WGT(1,7)=0.129484966168870 WGT(2,7)=0.279705391489277 WGT(3,7)=0.381830050505119 WGT(4,7)=0.417959183673469 WGT(5,7)=WGT(3,7) WGT(6,7)=WGT(2,7) WGT(7,7)=WGT(1,7) WGT(1,8)=0.101228536290376 WGT(2,8)=0.222381034453374 WGT(3,8)=0.313706645877887 WGT(4,8)=0.362683783378362 WGT(5,8)=WGT(4,8) WGT(6,8)=WGT(3,8) WGT(7,8)=WGT(2,8) WGT(8,8)=WGT(1,8) 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 D(1,1)=CB11 D(1,2)=CB12 D(1,3)=CB16 D(2,1)=CB12 D(2,2)=CB22 D(2,3)=CB26 D(3,1)=CB16 D(3,2)=CB26 D(3,3)=CB66 IF(ITYPE.EQ.1)THIC=1. IF(ITYPE.EQ.1)GO TO 10 C C FOR PLANE STRESS ANALYSIS C D(1,1)=CB11-CB13**2/CB33 D(2,2)=CB22-CB23**2/CB33 D(3,3)=CB66-CB36**2/CB33 D(1,2)=CB12-CB13*CB23/CB33 D(2,1)=D(1,2) D(2,3)=CB26-CB23*CB36/CB33 D(3,2)=D(2,3) D(1,3)=CB16-CB13*CB36/CB33 D(3,1)=D(1,3) C C CLEAR GLOBAL ARRAYS C 10 DO 1 I=1,NCOFS 1 A(I)=0.0 DO 2 I=1,NEQ 2 V(I)=0.0 C DO 333 J=1,3 333 WRITE(6,334)(D(J,I),I=1,3) 334 FORMAT(3E10.3) 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,ND5,ND6,ND7,ND8,THIC,TS1,TS2,TS3 READ(3)C(1,1),C(1,2),C(1,3),C(1,4),C(1,5),C(1,6),C(1,7),C(1,8), * C(2,1),C(2,2),C(2,3),C(2,4),C(2,5),C(2,6),C(2,7),C(2,8) C C CREATE CONNECTIVITY ARRAY LM(NDOFE) FOR EACH ELEMENT C NDOFE=16 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) LM(9)=ID(ND5,1) LM(10)=ID(ND5,2) LM(11)=ID(ND6,1) LM(12)=ID(ND6,2) LM(13)=ID(ND7,1) LM(14)=ID(ND7,2) LM(15)=ID(ND8,1) LM(16)=ID(ND8,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 30 AK(I,J)=0.0 C NINT=4 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 8-NODED QUADRILATERAL ELEMENTS C CALL STDM(C,B,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,16) 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 DO 40 K=1,3 DB(K)=0.0 DO 40 L=1,3 C DB(K)=DB(K)+D(K,L)*B(L,J) 40 CONTINUE C DO 60 I=J,NDOFE STIFF=0.0 DO 50 L=1,3 STIFF=STIFF+B(L,I)*DB(L) 50 CONTINUE C AK(I,J)=AK(I,J)+STIFF*WT C 60 CONTINUE C F(J)=F(J)+(B(1,J)*(D(1,1)*TS1+D(1,2)*TS2+D(1,3)*TS3) * +B(2,J)*(D(1,2)*TS1+D(2,2)*TS2+D(2,3)*TS3) * +B(3,J)*(D(1,3)*TS1+D(2,3)*TS2+D(3,3)*TS3)) 70 CONTINUE 80 CONTINUE C C SYMMETRY C DO 90 J=1,NDOFE DO 90 I=1,NDOFE 90 AK(J,I)=AK(I,J) C C ASSEMBLE EACH ELEMENT INTO THE GLOBAL MATRICES C 4 CALL ASSEMB(NDOFE,AK,F,LM,MAXA,A,V,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 RETURN 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 RETURN END C*********************************************************************** SUBROUTINE STDM(XX,B,DET,R,S,NEL) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * TO EVALUATE THE STRAIN-DISPLACEMENT TRANSFORMATION MATRIX B AT * C * POINT (R,S) FOR A QUADRILATERAL 8-NODED ELEMENT * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C IMPLICIT REAL*8(A-H,O-Z) C REAL NU COMMON/RDK1/LW,LR,KEY DIMENSION XX(2,8),B(3,16),H(8),P(2,8),XJ(2,2),XJI(2,2) C C EVALUATE INTERPOLATION COEFF.S AND JACOBIAN TRANSFORMATION MATRICES C CALL JACOB(NEL,R,S,XX,H,P,XJ,XJI,DET) C C CLEAR B(I,J) ARRAY C DO 40 I=1,3 DO 40 J=1,16 B(I,J)=0.0 40 CONTINUE C C EVALUATE GLOBAL DERIVATIVE OPERATOR B FOR NONSINGULAR PART C K2=0 DO 60 K=1,8 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 AN ISOPARAMETRIC QUADRILATERIAL ELEMENT. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY DIMENSION XX(2,8),H(8),P(2,8),XJ(2,2),XJI(2,2) C RP=1.0+R SP=1.0+S RM=1.0-R SM=1.0-S R2M=1.0-R**2 S2M=1.0-S**2 C C INTERPOLATION FUNCTIONS C H(1)=0.25*(RP*SP-S2M*SP-S2M*RP) H(2)=0.25*(RM*SP-R2M*SP-S2M*RM) H(3)=0.25*(RM*SM-S2M*RM-R2M*SM) H(4)=0.25*(RP*SM-R2M*SM-S2M*RM) H(5)=0.5*R2M*SP H(6)=0.5*S2M*RM H(7)=0.5*R2M*SM H(8)=0.5*S2M*RP C C NATURAL COORDINATE DERIVATIVES OF THE INTERPOLATION FUNCTIONS C C 1. WITH RESPECT TO R C P(1,1)=0.25*(SP+2.*R*SP-S2M) P(1,2)=0.25*(-SP+2.*R*SP+S2M) P(1,3)=0.25*(-SM+S2M+2.*R*SM) P(1,4)=0.25*(SM+2.*R*SM-S2M) P(1,5)=-R*SP P(1,6)=-S2M/2. P(1,7)=-R*SM P(1,8)=S2M/2. C C 2. WITH RESPECT TO S C P(2,1)=0.25*(RP-R2M+2.*S*RP) P(2,2)=0.25*(RM-R2M+2.*S*RM) P(2,3)=0.25*(-RM+2.*S*RM+R2M) P(2,4)=0.25*(-RP+R2M+2.*S*RP) P(2,5)=R2M/2. P(2,6)=-S*RM P(2,7)=-R2M/2. P(2,8)=-S*RP 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,8 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,MAXA,A,V,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 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DIMENSION AK(16,16),F(16),LM(16) DIMENSION MAXA(NEQ1),A(NCOFS),V(NEQ) 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) 10 CONTINUE C V(L)=V(L)+F(I) 20 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE NDLOAD(ID,V,NDS,NEQ,NLOAD) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * ADDS NODAL LOADS TO THE EXISTING ASSEMBLED GLOBAL * C * RIGHT-HAND-SIDE LOAD VECTOR V(NEQ) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),V(NEQ) C C USING EQN NO.S STORED IN ID(NDS,2) ARRAY C RECOVER EQN NO.(IEQN) WHICH CORRESPONDS TO NODAL LOAD C DIRECTION (IDIRN) AND ADD NODE LOAD (FLOAD) TO THE C RIGHT-HAND-SIDE LOAD VECTOR (V(IEQN)). C C WRITE HEADER ON OUTPUT FILE C WRITE(LW,10)NLOAD 10 FORMAT(1H1,//,' NUMBER OF APPLIED LOADS=',I4,/) C C L O O P T H R U L O A D S C DO 60 ILOAD=1,NLOAD C C READ 1. NOD , NODE NO. WHERE I-TH LOAD IS APPLIED C 2. IDIRN, DIRECTION OF I-TH LOAD C 3. FLOAD, MAGNITUDE OF I-TH LOAD C READ(LR,20)NOD,IDIRN,FLOAD 20 FORMAT(2I6,F12.0) C C RECOVER EQN NO. C IEQN=ID(NOD,IDIRN) C WRITE(LW,30)ILOAD,NOD,IDIRN,IEQN,FLOAD,IEQN,V(IEQN) 30 FORMAT(/,5X,'ILOAD=',I4,' NODE NO.=',I4,' LOAD DIRECTION=', *I3,' EQN NO.=',I5,/,10X,'LOAD MAGNITUDE=',E10.3,/, *15X,'OLD RIGHT-HAND-SIDE LOAD VECTOR, V(',I4,')=',E11.4) C C ADD NODAL LOAD (FLOAD) TO EXISTING LOAD VECTOR C V(IEQN)=V(IEQN)+FLOAD C WRITE(LW,50)IEQN,V(IEQN) 50 FORMAT(15X,'NEW RIGHT-HAND-SIDE LOAD VECTOR,V(',I4,')=',E11.4,/) C 60 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE NDDISP(ID,MAXA,A,V,NDS,NEQ1,NCOFS,NEQ,NDISP) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M C * 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 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),MAXA(NEQ1),A(NCOFS),V(NEQ) C C USING EQN NO.S STORED IN ID(NDS,2) ARRAY C RECOVER EQN NO. (IEQN) CORRESPONDING TO NODE DISPLACEMENT (DISP). C THEN 1. ADD LARGE SPRING STIFFNESS (HUGEK) TO DIAGONAL C TERN ( A(MAXA(IEQN)) ) OF THE GLOBAL STIFFNESS ARRAY. C 2. ADD LARGE SPRING FORCE ( HUGEK*DISP ) TO THE C RIGHT-HAND-SIDE LOAD VECTOR ( V(IEQN) ). C HUGEK=1.0E+20 C C PRIOR TO ADDITION OF HUGE SPRING STIFFNESS C STORE GLOBAL STIFFNESS ARRAY ON TAPE 8 C WRITE(8)(A(I),I=1,NCOFS) C C WRITE HEADER ON OUTPUT FILE C WRITE(LW,10)NDISP 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 DO 60 IDISP=1,NDISP 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 I-TH DISPL. C READ(LR,20)NOD,IDIRN,DISP 20 FORMAT(2I6,F12.0) C C STORE PRESCRIBED DISPL. DATA ON TAPE9 C WRITE(9)NOD,IDIRN,DISP C C RECOVER EQN NO. (IEQN) AND DIAGONAL ADDRESSES MAXA(IEQN) C IEQN=ID(NOD,IDIRN) II=MAXA(IEQN) C WRITE(LW,30)IDISP,NOD,IDIRN,IEQN,DISP,II,A(II),IEQN,V(IEQN) 30 FORMAT(/,5X,' IDISP=',I4,' NODE NO.=',I4,' DISPL. DIRECTION=', *I3,/,6X,'EQN NO.=',I5,14X,'DISPL. MAGNITUDE=',E10.3,//, *5X,'OLD DIAGONAL A(',I4,')=',E11.4, *' OLD RIGHT-HAND-SIDE VECTOR V(',I4,')=',E11.4) C C ADD SPRING STIFFNESS AND SPRING FORCE C A(II)=A(II)+HUGEK V(IEQN)=V(IEQN)+HUGEK*DISP C WRITE(LW,50)II,A(II),IEQN,V(IEQN) 50 FORMAT(5X,'NEW DIAGONAL A(',I4,')=',E11.4, *' NEW RIGHT-HAND-SIDE VECTOR V(',I4,')=',E11.4,/) C 60 CONTINUE 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 NDFORC(ID,MAXA,A,V,NDS,NEQ1,NCOFS,NEQ,NDISP,MBAND) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * RECOVER FORCES AT NODES CORRESPONDING TO PRESCRIBED * C * NODAL DISPLACEMENTS. FORCES CALCULATED FROM COMPONENTS * C * OF GLOBAL ARRAY A(I) STORED ON TAPE 8 PRIOR TO SUPERPOSING * C * A LARGE SPRING STIFFNESS (HUGEK) ON DIAGONAL TERMS. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),MAXA(NEQ1),A(NCOFS),V(NEQ) C C CALCULATE BANDWITH (MB) NOT INCLUDING DIAGONAL TERM C CALCULATE (LIMIT) EQN NO. WHERE TRUNCATION OF TERMS TO C RIGHT OF THE DIAGONAL OCCURS FOR (IEQN.GT.LIMIT) C MB=MBAND-1 LIMIT=NEQ1-MBAND C C REPLACE GLOBAL STIFFNESS ARRAY TERMS A(NCOFS) USED IN COLSOL WITH C VALUES STORED ON TAPE 8 PRIOR TO THE ADDITION OF HUGEK STIFFNESS. C REWIND 8 READ(8)(A(I),I=1,NCOFS) C C USING EQN NO.S STORED IN ID(NDS,3) AND MAXA(NEQ+1) ARRAYS C RECOVER EQN NO. (IEQN) CORRESPONDING TO PRESCRIBED NODAL DISPL. C C TO CALCULATE THE FORCE IN THE DIRECTION OF THE PRESCRIBED C DISPLACEMENT, SUM UP THE PRODUCT OF KNOWN DISPLACEMENTS STORED C IN ARRAY V(I) WITH THE ROW OF STIFFNESSES A(I) TERMS STORED C ON TAPE 8, WHICH COMPRISES EQN NO. (IEQN). C C WRITE HEADER ON OUTPUT FILE C WRITE(LW,10)NDISP 10 FORMAT(1H1,//,' NUMBER OF FORCES CORRESPONDING TO PRESCRIBED NODA *L DISPLACEMENTS=',I4) WRITE(LW,999) 999 FORMAT(/,80(1H*),/) C REWIND 9 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 DO 70 IDISP=1,NDISP C C RECOVER PRESCRIBED DISPL. DATA STORED ON TAPE 9 C NOD , NODE NO. WHERE I-TH DISPL. IS PRESCRIBED C IDIRN, DIRECTION OF I-TH DISPL. C DISP , MAGNITUDE OF I-TH DISPL. C READ(9)NOD,IDIRN,DISP C IEQN=ID(NOD,IDIRN) 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 DISPL. 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 V(I) AND STIFFNESSES A(I) C FOR TERMS TO THE LEFT AND INCLUDING DIAGONAL C FORCE=0.0 C DO 30 K=1,N FORCE=FORCE+V(IN)*A(IM) C IF(KEY.EQ.7)WRITE(LW,20)IN,V(IN),IM,A(IM),FORCE 20 FORMAT(/,' DISPLACEMENT, V(',I4,')=',E11.4, *' STIFFNESS, A(',I4,')=',E11.4,' FORCE=',E11.4) C IM=IM-1 IN=IN+1 30 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 55 C DO 50 IB=1,MB IM=MAXA(IN)+IB IF(IM.GE.MAXA(IN+1))GO TO 40 FORCE=FORCE+V(IN)*A(IM) IF(KEY.EQ.7)WRITE(LW,20)IN,V(IN),IM,A(IM),FORCE 40 IN=IN+1 50 CONTINUE C 55 WRITE(LW,60)NOD,IDIRN,IEQN,DISP,FORCE 60 FORMAT(///,' NODE NO.=',I5,' DIRECTION=',I2,' EQN. NO.=',I5,/, *10X,'PRESCRIBED DISPLACEMENT=',E11.4,' RESULTING FORCE=',E11.4) WRITE(LW,999) 70 CONTINUE C C E N D P R E S C R I B E D D I S P L. S L O O P C RETURN END C*********************************************************************** SUBROUTINE DISP(ID,U,V,NDS,NEQ) 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 U(NDS,2). * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION ID(NDS,2),U(NDS,2),V(NEQ) C C RECOVER NODAL DISPLACEMENTS C 5 IEQN=1 DO 10 I=1,NDS DO 10 J=1,2 U(I,J)=0.0 IF(ID(I,J).EQ.0)GO TO 10 U(I,J)=V(IEQN) IEQN=IEQN+1 10 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE STR(U,EH,EZ,EHZ,NE,NDS,D) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATES 1. STR AINS FROM DISPLACEMENTS * C * 2. STR ESSESES FROM STRAINS * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON/RDK1/LW,LR,KEY C C DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON C DIMENSION U(NDS,2),EH(NE,8,8),EZ(NE,8,8),EHZ(NE,8,8) DIMENSION B(3,16),D(3,3),IND(8),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,ND5,ND6,ND7,ND8,THIC,TS1,TS2,TS3 REWIND 1 C C TAPE 2 NELS,LXS,LXY,((B(I,J),I=1,3),J=1,16) 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 DO 333 J=1,3 333 WRITE(6,334)(D(J,I),I=1,3) 334 FORMAT(3E10.3) C DO 70 NEL=1,NE C C READ ELEMENT DATA READ(1)(IND(I),I=1,8),DUM,TS1,TS2,TS3 C C CHOSE NO. OF QUASS POINTS C NINT=4 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,16) 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 TOTAL STRAINS FROM DISPLACEMENTS C 10 E1=B(1,1)*U(IND(1),1)+B(1,2)*U(IND(1),2) * +B(1,3)*U(IND(2),1)+B(1,4)*U(IND(2),2) * +B(1,5)*U(IND(3),1)+B(1,6)*U(IND(3),2) * +B(1,7)*U(IND(4),1)+B(1,8)*U(IND(4),2) * +B(1,9)*U(IND(5),1)+B(1,10)*U(IND(5),2) * +B(1,11)*U(IND(6),1)+B(1,12)*U(IND(6),2) * +B(1,13)*U(IND(7),1)+B(1,14)*U(IND(7),2) * +B(1,15)*U(IND(8),1)+B(1,16)*U(IND(8),2) C E2=B(2,1)*U(IND(1),1)+B(2,2)*U(IND(1),2) * +B(2,3)*U(IND(2),1)+B(2,4)*U(IND(2),2) * +B(2,5)*U(IND(3),1)+B(2,6)*U(IND(3),2) * +B(2,7)*U(IND(4),1)+B(2,8)*U(IND(4),2) * +B(2,9)*U(IND(5),1)+B(2,10)*U(IND(5),2) * +B(2,11)*U(IND(6),1)+B(2,12)*U(IND(6),2) * +B(2,13)*U(IND(7),1)+B(2,14)*U(IND(7),2) * +B(2,15)*U(IND(8),1)+B(2,16)*U(IND(8),2) C E3=B(3,1)*U(IND(1),1)+B(3,2)*U(IND(1),2) * +B(3,3)*U(IND(2),1)+B(3,4)*U(IND(2),2) * +B(3,5)*U(IND(3),1)+B(3,6)*U(IND(3),2) * +B(3,7)*U(IND(4),1)+B(3,8)*U(IND(4),2) * +B(3,9)*U(IND(5),1)+B(3,10)*U(IND(5),2) * +B(3,11)*U(IND(6),1)+B(3,12)*U(IND(6),2) * +B(3,13)*U(IND(7),1)+B(3,14)*U(IND(7),2) * +B(3,15)*U(IND(8),1)+B(3,16)*U(IND(8),2) C C RECOVER MECHANICAL (ELASTIC) STRAINS FROM TOTAL STRAINS C 40 E(1)=E1-TS1 E(2)=E2-TS2 E(3)=E3-TS3 C C CALCULATE STRESSES FROM MECHANICAL (ELASTIC) STRAINS C DO 50 J=1,3 S(J)=0.0 DO 50 I=1,3 50 S(J)=S(J)+D(J,I)*E(I) C C STORE STRAINS IN ARRAYS C EH(NEL,LX,LY)=E(1) EZ(NEL,LX,LY)=E(2) EHZ(NEL,LX,LY)=E(3) C C WRITE STRESSES ON OUTPUT FILE CONDITIONAL ON KEY C IF((KEY.EQ.6).OR.(KEY.EQ.7))WRITE(LW,60)NEL,LX,LY,S(1),S(2),S(3) 60 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 WRITE(LW,65) 65 FORMAT(80(1H*)) C 70 CONTINUE DO 335 J=1,3 335 WRITE(6,334)(D(J,I),I=1,3) C C E N D E L E M E N T L O O P C C WRITE STRAINS AND DISPLACEMENTS ON OUTPUT FILE C IF((KEY.EQ.5).OR.(KEY.EQ.7))GO TO 75 RETURN C C STRAINS C 75 WRITE(LW,80) 80 FORMAT(1H1,//,25X,'** STRAINS **',/, * ' ELE. LX LY',7X,'EH',13X,'EZ',13X,'EHZ') C NINT=4 DO 91 NEL=1,NE DO 81 LX=1,NINT DO 81 LY=1,NINT 81 WRITE(LW,82)NEL,LX,LY,EH(NEL,LX,LY),EZ(NEL,LX,LY),EHZ(NEL,LX,LY) 82 FORMAT(2X,I5,2I3,3(1X,E14.7)) WRITE(LW,65) 91 CONTINUE C C DISPLACEMENTS C WRITE(LW,83) 83 FORMAT(1H1,//,16X,29H NODAL DISPLACEMENTS,/,10X,46HNODE * NO. U-DISPLACEMENT V-DISPLACEMENT) C DO 84 I=1,NDS 84 WRITE(LW,85)I,U(I,1),U(I,2) 85 FORMAT(13X,I5,2(7X,E12.5)) RETURN END