program sfem c PROGRAM SFEM(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: ( SFEM ) * C * S INGULAR F INITE E LEMENT M ODEL * 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 NU C COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E,NU,AL 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='sfem.dat',STATUS='OLD',ERR=888) OPEN(6,FILE='sfem.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)CL,SMH,SMZ 21 FORMAT(3F12.0,/) C WRITE(LW,22)NE,NDS,NLOAD,NDISP,CL,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' HORIZ. CRACK LENGTH ........................',F12.5,/, 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)E READ(LR,40)NU READ(LR,41)AL 40 FORMAT(8X,E10.3) 41 FORMAT(8X,E10.3,/) WRITE(LW,42)E 42 FORMAT(/,15X,'YOUNGS MODULUS E=',E10.3,/) WRITE(LW,43)NU 43 FORMAT(/,15X,'POISSONS RATIO NU=',E10.3,/) WRITE(LW,51)AL 51 FORMAT(/,15X,'THERMAL EXPANSION COEFFICIENT ALPHA=',E10.3,/) 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) READ(LR,61)JTYPE READ(LR,63)MTYPE 63 FORMAT(I6,/) 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 P R I N T H E A D E R F O R S I N G U L A R E L E M E N T S C IF(JTYPE.EQ.0)WRITE(LW,24) 24 FORMAT(//,' SINGULAR ELEMENTS ARE NOT PRESENT') IF(JTYPE.EQ.0)GO TO 25 WRITE(LW,23) 23 FORMAT(//,' SINGULAR ELEMENTS ARE PRESENT',/ * ' SINGULAR ELEMENTS REQUIRES EXTRA D.O.F.',//) IF(MTYPE.EQ.1)WRITE(LW,26) 26 FORMAT(' ONLY MODE-I S.I.F') IF(MTYPE.EQ.2)WRITE(LW,27) 27 FORMAT(' ONLY MODE-II S.I.F') IF(MTYPE.EQ.3)WRITE(LW,28) 28 FORMAT(' ONLY MODE-I AND MODE-II S.I.F') 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,CL,SMH,SMZ,JTYPE * ,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,NEQF1,NEQF2,JTYPE,MTYPE) 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 SINGULAR ELEMENTS REQUIRE TWO EXTRA EQN. NO.S, NEQF1 AND NEQF2 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,MTYPE,NEQF1,NEQF2 * ,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 ELEMNT(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),NDS,NEQ1,NCOFS * ,NEQ,ITYPE,MTYPE,NEQF1,NEQF2,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,NEQF1,NEQF2 * ,JTYPE,MTYPE,SIF1,SIF2) 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 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,SIF1,SIF2,D) C 888 STOP END C*********************************************************************** SUBROUTINE NEDATA(ID,H,Z,NDS,NE,CL,SMH,SMZ,JTYPE,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 COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E,NU,AL 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 TS=0.0 IF(LTYPE.EQ.2)TS=AL*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 IF(JTYPE.EQ.0)GO TO 40 IF(CL.EQ.0.0)GO TO 40 C C SHIFT GLOBAL COORDINATE ORIGIN TO CRACK TIP C H(I)=H(I)-CL C 40 CONTINUE C C WRITE SHIFTED NODE DATA ON OUTPUT FILE C IF(JTYPE.EQ.0)GO TO 45 IF(CL.EQ.0.0)GO TO 45 IF((KEY.EQ.2).OR.(KEY.EQ.7))GO TO 41 GO TO 45 41 WRITE(LW,42) 42 FORMAT(//,' ORIGIN OF GLOBAL COORD.S SHIFTED TO CRACK TIP',//, *5X,' S H I F T E D 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') DO 43 I=1,NDS 43 WRITE(LW,30)I,ID(I,1),ID(I,2),H(I),Z(I) 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 TYPE GRAD 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,KTYPE,MGRAD,THIC 60 FORMAT(6X,6I6,F12.0) C C STORE NODES ASSOCIATED WITH EACH ELEMENT ON TAPE 1 C WRITE(1)ND1,ND2,ND3,ND4,KTYPE,MGRAD,THIC,TS 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 Z1=Z(ND1)*SMZ Z2=Z(ND2)*SMZ Z3=Z(ND3)*SMZ Z4=Z(ND4)*SMZ 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,KTYPE * ,MGRAD,THIC 70 FORMAT(7(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,NEQF1,NEQF2,JTYPE,MTYPE) 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 NEQF1=0 NEQF2=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 C IF GRID DOES NOT INCLUDED SINGULAR ELEMENTS C THEN RETURN C IF(JTYPE.LE.0)RETURN C C WHEN SINGULAR ELEMENTS ARE PRESENT C ADD ADDITIONAL EQN.S TO THE SET OF SIMULTANEOUS EQUATIONS C C NEQ = TOTAL NO. OF EQN.S TO BE SOLVED C C NEQF1 = EQN. NO. ASSOCIATED WITH MODE I S.I.F. C NEQF1=0 IF((MTYPE.EQ.1).OR.(MTYPE.EQ.3))NEQ=NEQ+1 IF((MTYPE.EQ.1).OR.(MTYPE.EQ.3))NEQF1=NEQ C C NEQF2 = EQN. NO. ASSOCIATED WITH MODE II S.I.F. C NEQF2=0 IF(MTYPE.EQ.2)NEQ=NEQ+1 IF(MTYPE.EQ.2)NEQF2=NEQ IF(MTYPE.EQ.3)NEQ=NEQ+1 IF(MTYPE.EQ.3)NEQF2=NEQ C C ADDITIONAL EQN. ARE LOCATED AS THE LAST TWO EQN.S C OF THE ASSEMBLED (GLOBAL) STIFFNESS ARRAY SUCH THAT THE C COLUMN ORIENTED STORAGE OF THE NO. OF COEFF.S (NCOFS) IS MINIMIZED C C C WRITE ADDITIONAL EQUATION NUMBERS ON OUTPUT FILE C IF(MTYPE.EQ.1)WRITE(LW,50)NEQF1 50 FORMAT(' ADDITIONAL EQN:',I6,' ASSOCIATED WITH MODE I S.I.F.') IF(MTYPE.EQ.2)WRITE(LW,51)NEQF2 51 FORMAT(' ADDITIONAL EQN:',I6,' ASSOCIATED WITH MODE II S.I.F.') IF(MTYPE.EQ.3)WRITE(LW,52)NEQF1,NEQF2 52 FORMAT(' ADDITIONAL EQN.S:',I6,' ASSOCIATED WITH MODE I S.I.F.', * /,' ',I6,' ASSOCIATED WITH MODE II S.I.F.') C RETURN END C*********************************************************************** SUBROUTINE ADDRS(ID,MAXA,NE,NDS,NEQ1,MTYPE,NEQF1,NEQF2 * ,MBAND,NCOFS) 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(10) C C CLEAR LM AND MAXA ARRAYS C DO 5 IEQN=1,NEQ1 5 MAXA(IEQN)=0 DO 10 I=1,10 10 LM(I)=0 C C USING EQN. NO.S STORED IN ID(NDS,2) ARRAY TOGETHER WITH NODE C NO.S (ND1,ND2,ND3,ND4) STORED ON TAPE 1 CALCULATE CONNECTIVITY C 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,KTYPE,MGRAD,THIC,TS 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 FOR GRID WITH SINGULAR ELEMENTS INCLUDE EXTRA DOF.S C BY ASSINGING EQN NO. NEQF1 AND NEQF2 TO LM(9) AND/OR TO LM(10) C IF(MTYPE.EQ.1)LM(9) = NEQF1 IF(MTYPE.EQ.2)LM(9) = NEQF2 IF(MTYPE.EQ.3)LM(9) = NEQF1 IF(MTYPE.EQ.3)LM(10)= NEQF2 C C CHOOSE NO. OF DEGREES OF FREEDOM PER ELEMENT, NDOFE C C DEFINE NDOFE=8 WHEN GRID HAS NO SINGULAR ELEMENTS C NDOFE=8 C C DEFINE NDOFE=9 OR 10 WHEN GRID HAS SINGULAR ELEMENTS C IF((MTYPE.EQ.1).OR.(MTYPE.EQ.2))NDOFE=9 IF(MTYPE.EQ.3)NDOFE=10 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 ELEMNT(ID,MAXA,A,V,NDS,NEQ1,NCOFS,NEQ,ITYPE * ,MTYPE,NEQF1,NEQF2,NE,D) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATE ISOPARAMETRIC QUADRILATERAL ELEMENT STIFFNESS * C * MATRIX FOR PLANE STRESS AND PLANE STRAIN CONDITIONS. * C * ELEMENT TYPE: A (KTYPE=3) SINGULAR * C * ELEMENT TYPE: B (KTYPE=2) TRANSITIONAL * C * ELEMENT TYPE: C (KTYPE=1) NONSINGULAR * 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 * NEQF1 = NEQ+1: EQN. NO. ASSOC. WITH MODE I S.I.F. * C * NEQF2 = NEQ+2: EQN. NO. ASSOC. WITH NODE II S.I.F. * 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 NU COMMON/RDK1/LW,LR,KEY COMMON/RDK2/E,NU,AL DIMENSION ID(NDS,2),MAXA(NEQ1),A(NCOFS),V(NEQ) DIMENSION XG(8,8),WGT(8,8),D(3,3),C(2,4),B(3,10),DB(3) DIMENSION AK(10,10),F(10),LM(10) 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 Z=E/(1.+NU) G=Z*NU/(1.-2.*NU) H=Z+G C C FOR PLAIN STRAIN AYALYSIS C D(1,1)=H D(1,2)=G D(1,3)=0. D(2,1)=G D(2,2)=H D(2,3)=0. D(3,1)=0. D(3,2)=0. D(3,3)=Z/2. IF(ITYPE.EQ.1)THIC=1. IF(ITYPE.EQ.1)GO TO 10 C C FOR PLANE STRESS ANALYSIS C Z=E/(1-NU**2) G=(1.-NU)/2. C D(1,1)=Z D(2,2)=Z D(3,3)=Z*G D(1,2)=Z*NU D(2,1)=Z*NU D(2,3)=0. D(3,2)=0. D(3,1)=0. D(1,3)=0. 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 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,KTYPE,MGRAD,THIC,TS 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 IF ELEMENT TYPE IS SINGULAR ADD ADDITIONAL EQU. NO.S C IF(KTYPE.EQ.1)GO TO 20 IF((MTYPE.EQ.1).OR.(MTYPE.EQ.2))NDOFE=9 IF(MTYPE.EQ.3)NDOFE=10 IF(MTYPE.EQ.1)LM(9)=NEQF1 IF(MTYPE.EQ.2)LM(9)=NEQF2 IF(MTYPE.EQ.3)LM(9)=NEQF1 IF(MTYPE.EQ.3)LM(10)=NEQF2 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 IF(KTYPE.GE.2)NINT=7 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 TYPE(A,B,C) QUADRILATERAL ELEMENTS C CALL STDM(C,B,DET,RI,SI,NEL,ITYPE,E,NU,KTYPE,MGRAD) C C REASSIGN B(I,J) TERMS FOR CONSTRAINTS ON SINGULARITES WHEN MTYPE= 1 OR 2 C IF((MTYPE.EQ.0).OR.(MTYPE.EQ.3))GO TO 25 IF(MTYPE.EQ.2)GO TO 27 DO 26 I=1,3 B(I,10)=0.0 26 CONTINUE GO TO 25 27 DO 28 I=1,3 B(I,9)=B(I,10) B(I,10)=0.0 28 CONTINUE C C STORE B(I,J) ARRAY ON TAPE2 C 25 WRITE(2)NEL,LX,LY,((B(I,J),I=1,3),J=1,10) C C WRITE(LW,24)LX,LY C 24 FORMAT(' LX=',I2,' LY=',I2,/,' ********** B(I,J) **********',/) C DO 23 I=1,3 C 23 WRITE(LW,22)(B(I,J),J=1,10) C 22 FORMAT(10(1X,E10.3)) C C C WITH B(I,J) ARRAY, CALCULATE ELEMENTAL STIFFNESS MATRIX C WT=WGT(LX,NINT)*WGT(LY,NINT)*THIC*DET C C WRITE(LW,31)DET,THIC,LX,NINT,WGT(LX,NINT),LY,NINT,WGT(LY,NINT),WT C 31 FORMAT(' DET=',E11.4,' THIC=',E11.4,' WGT(',I2,',',I2,')=',E11.4 C *,' WGT(',I2,',',I2,')=',E11.4,/,' ******* WT=',E11.4) C DO 70 J=1,NDOFE DO 40 K=1,3 DB(K)=0.0 DO 40 L=1,3 C C IF(NEL.EQ.1)WRITE(LW,32)K,DB(K),K,L,D(K,L),L,J,B(L,J) C 32 FORMAT(' DB(',I2,')=',E11.4,' D(',I2,',',I2,')=',E11.4, C *' B(',I2,',',I2,')=',E11.4) 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 C AK(I,J)=AK(I,J)+STIFF*WT C C IF(NEL.EQ.1)WRITE(LW,34)I,J,AK(I,J),STIFF,WT C 34 FORMAT(' AK(',I2,',',I2,')=',E11.4,' STIFF=',E11.4,' WT=',E11.4) C 60 CONTINUE C F(J)=F(J)-(B(1,J)*(D(1,1)+D(1,2)) * +B(2,J)*(D(1,2)+D(2,2)) * +B(3,J)*(D(1,3)+D(2,3)))*TS 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 WRITE(LW,94)NEL,KTYPE,MGRAD C 94 FORMAT(' NEL=',I6,' KTYPE=',I2,' MGRAD=',I2,/,' AK(I,J)',/) C DO 91 I=1,NDOFE C IF(NDOFE.EQ.8)WRITE(LW,93)(AK(I,J),J=1,NDOFE) C IF(NDOFE.EQ.10)WRITE(LW,92)(AK(I,J),J=1,NDOFE) C 91 CONTINUE C 92 FORMAT(10(1X,E10.3)) C 93 FORMAT(8(1X,E10.3)) 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,ITYPE,E,NU,KTYPE,MGRAD) 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: * C * TYPE C QUADRILATERAL NONSINGULAR ELEMENT (KTYPE=1) * C * TYPE B QUADRILATERAL TRANSITION ELEMENT (KTYPE=2) * C * TYPE A QUADRILATERAL SINGULAR ELEMENT (KTYPE=3) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C IMPLICIT REAL*8(A-H,O-Z) C REAL NU COMMON/RDK1/LW,LR,KEY DIMENSION XX(2,4),B(3,10),H(4),P(2,8),XJ(2,2),XJI(2,2) DIMENSION XY(2),DQDX(2,2),DQDY(2,2),DQDR(2,2),DQDS(2,2) DIMENSION RT(2),RTB(2),Q(2,2),QB(2,2,4) 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,10 B(I,J)=0.0 40 CONTINUE C C EVALUATE GLOBAL DERIVATIVE OPERATOR B FOR NONSINGULAR PART 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 C END CALCULATIONS FOR TYPE C (KTYPE=1) ELEMENTS C IF(KTYPE.EQ.1)RETURN C C E V A L U A T E S I N G U L A R T E R M S O F B(I,J) C C F O R T Y P E A S I N G U L A R E L E M E N T S C C TRANSFORM QUASS POINTS (RI,SI) TO GLOBAL POLAR COORDS. (R,THETA) C C CALCULATE GLOBAL COORDS. (X,Y) AT QUASS POINTS (RI,SI) C DO 80 I=1,2 DUM=0.0 DO 70 J=1,4 70 DUM=DUM+H(J)*XX(I,J) 80 XY(I)=DUM C C TRANSFORM (X,Y) GLOBAL COORDS.OF THE GAUSS POINTS (R,S) TO (R,THETA) COORDS. C DUM=XY(1)**2+XY(2)**2 IF(DUM.EQ.0.0)WRITE(6,300)XY(1),XY(2) 300 FORMAT(' ZERO RADIUS: X=',E11.4,' Y=',E11.4) IF(DUM.EQ.0.0)STOP RT(1)=SQRT(DUM) RT(2)=ATAN2(XY(2),XY(1)) C C WRITE(LW,310)R,S,XY(1),XY(2),RT(1),RT(2) C 310 FORMAT(' R=',E11.4,' S=',E11.4,' X=',E11.4,' Y=',E11.4,/, C * ' RT(1)=',E11.4,' RT(2)=',E11.4) C C EVALUATE Q(I,J) AND THEIR X,Y DERIVATIVES AT (R,THETA) QUASS POINTS C PI=ACOS(-1.) PI2=PI*2. IF(ITYPE.EQ.1)RK=3.-4.*NU IF(ITYPE.EQ.2)RK=(3.-NU)/(1.+NU) G=E/(2.*(1.+NU)) DEN=2.*G*SQRT(2.*PI*RT(1)) CON=SQRT(RT(1)/PI2)/G CT=COS(RT(2)) ST=SIN(RT(2)) CTD2=COS(RT(2)/2.) STD2=SIN(RT(2)/2.) C RUM=(RK-1)/2.+STD2**2 DQDX(1,1)=(CT*CTD2*RUM-ST*STD2*(CTD2**2-RUM/2.)*2.)/DEN DQDY(1,1)=(ST*CTD2*RUM+CT*STD2*(CTD2**2-RUM/2.)*2.)/DEN Q(1,1)=CON*CTD2*RUM C C WRITE(LW,320)DQDX(1,1),DQDY(1,1),Q(1,1) C 320 FORMAT(' DQDX(1,1)=',E11.4,' DQDY(1,1)=',E11.4,' Q(1,1)=',E11.4) C RUM=(RK+1)/2.+CTD2**2 DQDX(2,1)=(CT*STD2*RUM-ST*CTD2*(-STD2**2+RUM/2.)*2.)/DEN DQDY(2,1)=(ST*STD2*RUM+CT*CTD2*(-STD2**2+RUM/2.)*2.)/DEN Q(2,1)=CON*STD2*RUM C C WRITE(LW,330)DQDX(2,1),DQDY(2,1),Q(2,1) C 330 FORMAT(' DQDX(2,1)=',E11.4,' DQDY(2,1)=',E11.4,' Q(2,1)=',E11.4) C RUM=(RK+1)/2.-CTD2**2 DQDX(1,2)=(CT*STD2*RUM-ST*CTD2*(STD2**2+RUM/2.)*2.)/DEN DQDY(1,2)=(ST*STD2*RUM+CT*CTD2*(STD2**2+RUM/2.)*2.)/DEN Q(1,2)=CON*STD2*RUM C C WRITE(LW,340)DQDX(1,2),DQDY(1,2),Q(1,2) C 340 FORMAT(' DQDX(1,2)=',E11.4,' DQDY(1,2)=',E11.4,' Q(1,2)=',E11.4) C RUM=(1-RK)/2.+STD2**2 DQDX(2,2)=(CT*CTD2*RUM-ST*STD2*(CTD2**2-RUM/2.)*2.)/DEN DQDY(2,2)=(ST*CTD2*RUM+CT*STD2*(CTD2**2-RUM/2.)*2.)/DEN Q(2,2)=CON*CTD2*RUM C C WRITE(LW,350)DQDX(2,2),DQDY(2,2),Q(2,2) C 350 FORMAT(' DQDX(2,2)=',E11.4,' DQDY(2,2)=',E11.4,' Q(2,2)=',E11.4) C C EVALUATE DERIVATIVES W.R.T. LOCAL COORDS. (R,S) AT QUASS POINTS (RI,SI) C DO 90 I=1,2 DO 90 J=1,2 DQDR(I,J)=XJ(1,1)*DQDX(I,J)+XJ(1,2)*DQDY(I,J) DQDS(I,J)=XJ(2,1)*DQDX(I,J)+XJ(2,2)*DQDY(I,J) C C WRITE(LW,360)I,J,DQDR(I,J),I,J,DQDS(I,J) C 360 FORMAT(' DQDR(',I2,',',I2,')=',E11.4,' DQDS(',I2,',',I2,')=' C *,E11.4) C 90 CONTINUE C C EVALUATE SINGULAR FUNCTIONS AT NODE POINTS C DO 100 K=1,4 XY(1)=XX(1,K) XY(2)=XX(2,K) RTB(1)=SQRT(XY(1)**2+XY(2)**2) RTB(2)=0.0 IF(RTB(1).NE.0.0)RTB(2)=ATAN2(XY(2),XY(1)) CONB=SQRT(RTB(1)/PI2)/G CTD2B=COS(RTB(2)/2.) STD2B=SIN(RTB(2)/2.) QB(1,1,K)=CONB*CTD2B*((RK-1)/2.+STD2B**2) QB(2,1,K)=CONB*STD2B*((RK+1)/2.+CTD2B**2) QB(1,2,K)=CONB*STD2B*((RK+1)/2.-CTD2B**2) QB(2,2,K)=CONB*CTD2B*((1-RK)/2.+STD2B**2) C C WRITE(LW,370)K,XY(1),XY(2),RTB(1),RTB(2),QB(1,1,K),QB(2,1,K) C * ,QB(1,2,K),QB(2,2,K) C 370 FORMAT(' NODE NO.=',I1,' X=',E11.4,' Y=',E11.4,/, C * ' RT(1)=',E11.4,' RT(2)=',E11.4,/, C *' QB(1,1)=',E11.4,' QB(2,1)=',E11.4,' QB(1,2)=',E11.4, C *' QB(2,2)=',E11.4) C 100 CONTINUE C C CALCULATE COEFFS. OF KI AND KII. DERIV.S ARE IN LOCAL COORD.S (R,S) C DO 120 I=1,2 DO 120 J=1,2 SUMR=DQDR(I,J) SUMS=DQDS(I,J) DO 110 K=1,4 SUMR=SUMR-P(1,K)*QB(I,J,K) SUMS=SUMS-P(2,K)*QB(I,J,K) 110 CONTINUE IF((I.EQ.1).AND.(J.EQ.1))P(1,5)=SUMR IF((I.EQ.1).AND.(J.EQ.1))P(2,5)=SUMS IF((I.EQ.2).AND.(J.EQ.1))P(1,6)=SUMR IF((I.EQ.2).AND.(J.EQ.1))P(2,6)=SUMS IF((I.EQ.1).AND.(J.EQ.2))P(1,7)=SUMR IF((I.EQ.1).AND.(J.EQ.2))P(2,7)=SUMS IF((I.EQ.2).AND.(J.EQ.2))P(1,8)=SUMR IF((I.EQ.2).AND.(J.EQ.2))P(2,8)=SUMS 120 CONTINUE C C WRITE(LW,380)P(1,5),P(1,6),P(1,7),P(1,8) C WRITE(LW,390)P(2,5),P(2,6),P(2,7),P(2,8) C 380 FORMAT(' P15',E11.4,'P16=',E11.4,'P17=',E11.4,'P18=',E11.4) C 390 FORMAT(' P25',E11.4,'P26=',E11.4,'P27=',E11.4,'P28=',E11.4) C C TEST FOR ELEMENT TYPE(KTYPE) C IF(KTYPE.EQ.3)GO TO 140 C C C U P A D A T E S I N G U L A R T E R M S O F B(I,J) C C F O R T Y P E B T R A N S I T I O N E L E M E N T S C C TEST GRADIENT OF TYPE B ELEMENT C IF(MGRAD.GT.0)GO TO 146 C WRITE(LW,145)NEL,KTYPE,MGRAD 145 FORMAT(' GRADIENT NOT DEFINED FOR TRANSITION ELEMENT=',I6,/, * ' KTYPE=',I2,' MGRAD=',I2) STOP C C CALCULATE TRANSITION COEFFICIENTS FOR TYPE B ELEMENTS C 146 R1=RM(MGRAD,R,S) R2=RMR(MGRAD,R,S) R3=RMS(MGRAD,R,S) C C WRITE(LW,391)R1,R2,R3 C 391 FORMAT(' RM=',E11.4,' RMR=',E11.4,' RMS=',E11.4) C C C CALCULATE COEFF.S OF KI AND KII FOR TRANSITIONAL ELEMENT TYPE C W1=Q(1,1) W2=Q(1,2) Z1=Q(2,1) Z2=Q(2,2) DO 130 K=1,4 W1=W1-H(K)*QB(1,1,K) W2=W2-H(K)*QB(1,2,K) Z1=Z1-H(K)*QB(2,1,K) Z2=Z2-H(K)*QB(2,2,K) 130 CONTINUE C C WRITE(LW,400)W1,W2,Z1,Z2 C 400 FORMAT(' W1=',E11.4,' W2=',E11.4,' Z1=',E11.4,' Z2=',E11.4) C C UPDATE P(1-2,5-8) COEFF.S FOR KI AND KII TERMS C P(1,5)=R1*P(1,5)+R2*W1 P(2,5)=R1*P(2,5)+R3*W1 P(1,6)=R1*P(1,6)+R2*Z1 P(2,6)=R1*P(2,6)+R3*Z1 P(1,7)=R1*P(1,7)+R2*W2 P(2,7)=R1*P(2,7)+R3*W2 P(1,8)=R1*P(1,8)+R2*Z2 P(2,8)=R1*P(2,8)+R3*Z2 C C WRITE(LW,380)P(1,5),P(1,6),P(1,7),P(1,8) C WRITE(LW,390)P(2,5),P(2,6),P(2,7),P(2,8) C C CALCULATE ADDITIONAL B(I,J) TERMS CORRESPONDING TO KI AND KII S.I.F. C 140 DO 150 I=1,2 B(1,9)= B(1,9)+ XJI(1,I)*P(I,5) B(1,10)=B(1,10)+XJI(1,I)*P(I,6) B(2,9)= B(2,9)+ XJI(2,I)*P(I,7) B(2,10)=B(2,10)+XJI(2,I)*P(I,8) B(3,9)= B(3,9)+ XJI(1,I)*P(I,7)+XJI(2,I)*P(I,5) B(3,10)=B(3,10)+XJI(1,I)*P(I,8)+XJI(2,I)*P(I,6) 150 CONTINUE C RETURN END C*********************************************************************** FUNCTION RM(MGRAD,R,S) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATE TRANSITION COEFF. RM FOR GRADIENT NO. (MGRAD) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C GO TO(1,2,3,4,5,6,7,8),MGRAD C 1 RM=(1.-R)/2. GO TO 9 C 2 RM=(1.+R)/2. GO TO 9 C 3 RM=(1.-S)/2. GO TO 9 C 4 RM=(1.+S)/2. GO TO 9 C 5 RM=(1.-R)*(1.-S)/4. GO TO 9 C 6 RM=(1.-R)*(1.+S)/4. GO TO 9 C 7 RM=(1.+R)*(1.+S)/4. GO TO 9 C 8 RM=(1.+R)*(1.-S)/4. C 9 RETURN END C*********************************************************************** FUNCTION RMR(MGRAD,R,S) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATE DERIVATIVE OF THE TRANSITION COEFF. (RM) W.R.T. R * C * CORRESPONDING TO THE GRADIENT NO. (MGRAD) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C GO TO(1,2,3,4,5,6,7,8),MGRAD C 1 RMR=-0.5 GO TO 9 C 2 RMR=+0.5 GO TO 9 C 3 RMR=0.0 GO TO 9 C 4 RMR=0.0 GO TO 9 C 5 RMR=-0.25*(1.-S) GO TO 9 C 6 RMR=-0.25*(1.+S) GO TO 9 C 7 RMR=+0.25*(1.+S) GO TO 9 C 8 RMR=+0.25*(1.-S) C 9 RETURN END C*********************************************************************** FUNCTION RMS(MGRAD,R,S) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * CALCULATE DERIVATIVE OF THE TRANSITION COEFF. (RM) W.R.T. S * C * CORRESPONDING TO THE GRADIENT NO. (MGRAD) * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C GO TO(1,2,3,4,5,6,7,8),MGRAD C 1 RMS=0.0 GO TO 9 C 2 RMS=0.0 GO TO 9 C 3 RMS=-0.5 GO TO 9 C 4 RMS=+0.5 GO TO 9 C 5 RMS=-0.25*(1.-R) GO TO 9 C 6 RMS=+0.25*(1.-R) GO TO 9 C 7 RMS=+0.25*(1.+R) GO TO 9 C 8 RMS=-0.25*(1.+R) C 9 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,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 JOCOBIAN 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(10,10),F(10),LM(10) 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 PRIESCRIBEC 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 PRESCRIBEC 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,NEQF1,NEQF2,JTYPE,MTYPE,SIF1,SIF2) C CALLING PROGRAM - MAIN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * P R O G R A M * C * 1. RECOVERS NODAL DISPL.S STORED IN V(NEQ) BY COLSOL. * C * AND ASSIGNS NODE DISPL.S TO ARRAY U(NDS,2). * C * 2. RECOVERS STRESS INTENSITY FACTORS FROM LAST TWO EQUAT.S * 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 MODE-I AND MODE-II STRESS INTENSITY FACTORS C IF SINGULAR ELEMENTS ARE PRESENT C IF(JTYPE.EQ.0)GO TO 5 IF(MTYPE.EQ.3)SIF1=V(NEQF1) IF(MTYPE.EQ.3)SIF2=V(NEQF2) IF(MTYPE.EQ.2)SIF1=0.0 IF(MTYPE.EQ.2)SIF2=V(NEQF2) IF(MTYPE.EQ.1)SIF1=V(NEQF1) IF(MTYPE.EQ.1)SIF2=0.0 WRITE(LW,1)SIF1,SIF2 1 FORMAT(1H1,' STRESS INTENSITY FACTORS',/, * ' MODE-I = ',E12.5,/, * ' MODE-II = ',E12.5,///) 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,SIF1,SIF2,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,10),D(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,KTYPE,MGRAD,THIC,TS REWIND 1 C C TAPE 2 NELS,LXS,LXY,((B(I,J),I=1,3),J=1,10) REWIND 2 C C WRITE HEADER FOR STRESSES C 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),KTYPE,IDUM,DUM2,TS C C CHOSE NO. OF QUASS POINTS C NINT=4 IF(KTYPE.GE.2)NINT=7 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,10) 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) 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) 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) C C IF SINGULAR ELEMENTS ARE PRESENT ADD CONTRIBUTION OF SINGULAR TERMS C IF(KTYPE.EQ.1)GO TO 40 E1=E1+B(1,9)*SIF1+B(1,10)*SIF2 E2=E2+B(2,9)*SIF1+B(2,10)*SIF2 E3=E3+B(3,9)*SIF1+B(3,10)*SIF2 C C RECOVER MECHANICAL (ELASTIC) STRAINS FROM TOTAL STRAINS C 40 E(1)=E1-TS E(2)=E2-TS E(3)=E3-TS 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 C 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 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 REWIND 1 DO 91 NEL=1,NE READ(1)(IND(I),I=1,4),KTYPE,IDUM,DUM2,DUM3 NINT=4 IF(KTYPE.GE.2)NINT=7 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