PROGRAM DYNA2 C *** TIME INTEGRATION IMPLICIT ALGORITHM C * PLANE STRAIN (CONSTANT STRAIN 3 AND 4 NODED ELEMENTS) C * ORDER OF STRAINS: X, Y, XY, Z C * NON-ASSOCIATED FLOW RULE C * SHEAR-BANDING STRESS-STRAIN RELATIONSHIP C * COMPRESSIVE STRESS IS POSITIVE C * GROUND-WATER LEVEL (STATIC PORE-WATER PRESSURE) C CAN BE CONSIDERED. C * WHEN USING INTERFACE-ELEMENTS IN ACTUAL SIZE MODEL, C DIVERGENCE MAY OCCUR DUE TO VIBRATION BY GRAVITY FORCE, C BECAUSE THE ELEMENTS BEAR NO TENSILE STRESS. C * IN ACTUAL SIZE MODEL, CONSIDER GRAVITY FORCE. C TIME LENGTH OF A FEW SECONDS IS REQUIRED FOR STABILIZING C THE EFEECT OF GRAVITY FORCE C IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA11/IREQD(10),IREQS(10),IWA(1260) COMMON/LA12/ACCEH(6600),ACCEV(6600) COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA14/XMASS(1850,1850),STIFF(1850,1850),STIFS(1850,1850) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA17/VELOI(1850),VELOL(1850),VELOT(1850) COMMON/LA18/ACCEI(1850),ACCEL(1850) COMMON/LA19/DAMPI(1850,1850),RESID(1850),STRIN(1260,3) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) COMMON/LA22/NTENS,KTENS,ITENS(1260),JTENS(1260) COMMON/LA23/TMX(1260,3,3),TMI(1260,3,3),EFFST(1260) C OPEN(5,FILE='DADYNA-2',STATUS='UNKNOWN') OPEN(6,FILE='PRDYNA-2',STATUS='UNKNOWN') OPEN(8,FILE='PRDISP',STATUS='UNKNOWN') OPEN(9,FILE='PRVELO',STATUS='UNKNOWN') OPEN(10,FILE='PRACCE',STATUS='UNKNOWN') OPEN(11,FILE='M-AVS1.INP',STATUS='UNKNOWN') OPEN(12,FILE='M-AVS2.INP',STATUS='UNKNOWN') OPEN(13,FILE='M-AVS3.INP',STATUS='UNKNOWN') C CALL INPUTD CALL LINKIN CALL INTIME CALL PREVOS CALL LOADPL CALL LUMASS CALL GSTIFF C DO 10 ISTEP=1,NSTEP DO 20 IITER=1,MITER CALL IMPEXP(IITER,ISTEP) CALL RESEPL(IITER,ISTEP) CALL ITRATE(IITER) IF(NCHEK.EQ.1) GO TO 15 20 CONTINUE 15 CALL OUTDYN(IITER,ISTEP) IOUT=NOUTD(ISTEP) IF(IOUT.LE.0) GO TO 10 CALL MAVS1 CALL MAVS2 CALL MAVS3 10 CONTINUE C CLOSE(5) CLOSE(6) CLOSE(8) CLOSE(9) CLOSE(10) CLOSE(11) CLOSE(12) CLOSE(13) STOP END C * * * * * * * * * * * * SUBROUTINE INPUTD C *** INPUT ROUTINE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA11/IREQD(10),IREQS(10),IWA(1260) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) DIMENSION TITLE(20) DIMENSION IDA(10),IDAQ(10),NET(10) C C NN: MAXIMUM NUMBER OF DISPLACEMENT VARIABLES NN=1850 DO 2 I=1,10 2 IDAQ(I)=0 DO 4 I=1,6 4 NET(I)=0 C TITLE(I): PROJECT NAME READ(5,100) (TITLE(I),I=1,15) WRITE(6,100) (TITLE(I),I=1,15) C NPOIN: NUMBER OF NODES, NELEM:NUMBER OF ELEMENTS C NMATS: NUMBER OF MATERIALS C NIST : NUMBER OF INITIAL STRESS GROUPS C NPREV=1: CONSIDER INITIAL STRESSES IN EACH ELEMENT C 0: NOT CONSIDER C ICRI=1: ELASTIC-PERFECTLY PLASTIC C 0: ELASTIC-PLASTIC C IPR=0: NOT PRINT NODAL & ELEMENT DATA, 1:PRINT C IPW=0: NOT CONSIDER WATER PRESSURE C 1: CALCULATE PRE-EFFECTIVE STRESSES C PERFORM FE-ANALYSIS USING UNDER-WATER UNIT-WEIGHT C 2: CALCULATE & CONSIDER WATER PRESSURE C READ(5,101) NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,IPR,IPW WRITE(6,200) NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,IPR,IPW C *** NODAL DATA C K:NODE NO., XX(I):X-COORDINATE, YY(I):Y-COORDINATE C IX(I)=0: X-DIRECTION FREE, 1=FIXED C IY(I)=0: Y-DIRECTION FREE, 1=FIXED C IQ(I)=0: ROTATION FREE, 1=FIXED C WRITE(6,201) DO 10 IP=1,NPOIN READ(5,102) K,XX(IP),YY(IP),IX(IP),IY(IP),IQ(IP) IF(IPR.EQ.0) GO TO 10 WRITE(6,102) K,XX(IP),YY(IP),IX(IP),IY(IP),IQ(IP) 10 CONTINUE C *** ELEMENT DATA C K: ELEMENT NO., IJK(L,I): NODE NO. C K2(IE)=ELEMENT TYPE: TRUSS=1, BEAM=2, TEXTILE=3, C PLANE-STRAIN (ELASTIC=4, ELASTIC-PLASTIC=5), C INTERFACE=6 (I-J:SHEAR DIRECTION, ANTI-CLOCKWISE) C MAT(IE): NO. OF MATERIAL C IST(IE): NO. OF SET OF ACTUAL INITIAL STRESSES C LRE(IE): DIRECTION OF SHEAR-BAND C =1: LEFT SIDE, 2: RIGHT SIDE AT MAXIMUM PRINCIPAL STRESS C *** IF DATA(I)=0, PRECEDING VALUE IS EMPLOYED C DO 12 IE=1,NELEM READ(5,103) K,(IJK(IE,I),I=1,4),(IDA(I),I=1,5) DO 14 I=1,5 IF(IDA(I).EQ.0) IDA(I)=IDAQ(I) IF(IDA(I).LT.0) IDA(I)=0 IDAQ(I)=IDA(I) 14 CONTINUE K2(IE)= IDA(1) MAT(IE)=IDA(2) IST(IE)=IDA(3) LRE(IE)=IDA(4) IWA(IE)=IDA(5) N=K2(IE) NET(N)=NET(N)+1 12 CONTINUE WRITE(6,202) (NET(I),I=1,6) WRITE(6,203) DO 20 IE=1,NELEM IF(IPR.EQ.0) GO TO 20 WRITE(6,204) IE,(IJK(IE,I),I=1,4), * K2(IE),MAT(IE),IST(IE),LRE(IE),IWA(IE) 20 CONTINUE C *** MATERIAL PROPERTIES C TRUSS 1:YOUNG'S MODULUS, 2:CROSS AREA, 3:0, 4:RHO, 5:0, 6:0, 7:0 C 8:ALFA, 9:BETA C BEAM 1;YOUNG'S MODULUS, 2:CROSS AREA, 3:MOMENT OF INERTIA, 4:RHO C 5:0, 6:0, 7:0, 8:ALFA, 9:BETA C PLANE-STRAIN 1:YOUNG'S MODULUS, 2:POISSON'S RATIO, 3:1.0, 4:DENSITY C 5:COHESION C, 6:ANGLE OF SHEAR RESISTANCE, 7:ANGLE OF DILATANCY C 8:ALFA, 9:BETA (DAMPING PARAMETERS) C 10:RATIO OF REDUCING COMPRESSIVE STRENGTH TO TENSILE STRENGTH C INTERFACE 1:SHEAR RIGIDITY, 2:YOUNG'S MODULUS, 3:POISSON'S RATIO, C 4:DENSITY, 5:COHESION C, 6:ANGLE OF SHEAR RESISTANCE C 7:ANGLE OF DILATANCY,8:ALFA, 9:BETA C CALCULATE DAMPING PARAMETERS, ACCORDING TO SPECIFIC PERIOD. C WRITE(6,205) DO 30 MT=1,NMATS READ(5,105) MT1,(PRP(MT,I),I=1,10) WRITE(6,206) MT1,(PRP(MT,I),I=1,10) COHES=PRP(MT,5) COMST=2.D0*COHES TENST=COMST*PRP(MT,10) C HEREAFTER PRP(MT,10) MEANS TENSILE STRENGTH PRP(MT,10)=TENST 30 CONTINUE C C *** PORE WATER PRESSURE WRITE(6,207) DO 40 IE=1,NELEM DO 40 I=1,3 40 PW(IE,I)=0.D0 IF(IPW.LE.1) GO TO 50 C USE STRSG(IE,I) & STRAG(IE,I) TEMPORALILY DO 42 I=1,2 42 READ(5,106) (STRSG(IE,I),IE=1,NELEM) DO 44 I=1,2 44 READ(5,106) (STRAG(IE,I),IE=1,NELEM) DO 46 IE=1,NELEM IF(IWA(IE).LE.0) GO TO 46 DO 47 I=1,2 47 PW(IE,I)=STRSG(IE,I)-STRAG(IE,I) 46 CONTINUE DO 48 I=1,3 48 WRITE(6,106) (PW(IE,I),IE=1,NELEM) 50 CONTINUE C 100 FORMAT(15A4) 101 FORMAT(15I5) 102 FORMAT(I5,2F10.3,3I5) 103 FORMAT(16I5) 104 FORMAT(8E10.3) 105 FORMAT(I2,10E10.3) 106 FORMAT(10F8.3) 200 FORMAT(/'NODES=',I5,' ELEMENTS=',I5,' MATERIALS=',I3 */' INITIAL STRESS GROUPS=',I3 */' GIVE INITIAL STRESSES IN EACH ELEMENT (YES=1)=',I2 */' PERFECTLY PLASTIC=',I2 */' PRINT=',I2,' PORE-WATER PRESSURE=',I2) 201 FORMAT(/'NODAL DATA'/' NO.',9X,'X',9X,'Y',' X-F Y-F R-F') 202 FORMAT(/'TRUSS=',I3,' BEAM=',I3,' TEXTILE=',I3,' P-STRA=',I3 */' P-STRA(NL)=',I3,' INTER=',I3) 203 FORMAT(/'ELEMENT DATA'/' NO. NODES TYPE MAT IST L-R') 204 FORMAT(16I5) 205 FORMAT(/'MATERIAL PARAMETERS' * /'NO. E NYU THICK RHO C FAI DELTA ALFA BETA TENS'/) 206 FORMAT(I2,10E10.3) 207 FORMAT(/'PORE-WATER PRESSURE') RETURN END C * * * * * * * * * * SUBROUTINE LINKIN C *** LINKS WITH PROFILE SOLVER IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) DIMENSION LX(3500),NDF(1280) C DO 10 IP=1,NPOIN 10 NDF(IP)=2 DO 12 IE=1,NELEM IF(K2(IE).NE.2) GO TO 12 I=IJK(IE,1) J=IJK(IE,2) NDF(I)=3 NDF(J)=3 12 CONTINUE DO 20 IP=1,NPOIN LY(3*IP-2)=IX(IP) LY(3*IP-1)=IY(IP) IF(NDF(IP).EQ.2) IQ(IP)=1 20 LY(3*IP)=IQ(IP) NX=3*NPOIN LOC=0 DO 22 L=1,NX IF(LY(L).EQ.0) GO TO 24 LY(L)=NN GO TO 22 24 LOC=LOC+1 LX(LOC)=L LY(L)=LOC 22 CONTINUE NSIZE=LOC WRITE(6,200) NSIZE 200 FORMAT(/'NUMBER OF VARIABLES',I5) RETURN END C * * * * * * * * * * * * SUBROUTINE INTIME C *** INITIAL VALUES AND TIME INTEGRATION DATA IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA11/IREQD(10),IREQS(10),IWA(1260) COMMON/LA12/ACCEH(6600),ACCEV(6600) DIMENSION LOUTD(100) C C *** READ TIME STEPPING AND SELECTIVE OUTPUT PARAMETERS C NSTEP: NUMBER OF TOTAL TIME-STEPS C NOUTP: NUMBER OF SELECTED TIME-STEPS FOR OUTPUT C NREQD: NUMBER OF SELECTED NODES FOR OUTPUT C NREQS: NUMBER OF SELECTED ELEMENTS FOR OUTPUT C NACCE: NUMBER OF STEPS FOR ACCERALATION INPUT C IFUNC=0: INPUT ACCERALATION TIME RECORD C =1: HEAVISIDE FUNCTION C =2: FUNCTION=AZERO+BZERO*SIN(ARGUM) C ARGUM=OMEGA*JSTEP*DTIME C IFIXD: DIRECTION OF ACCERALATION INPUT C =0: HORIZONTAL & VERTICAL, 1:VERTICAL, 2:HORIZONTAL C MITER: MAXIMUM NUMBER OF ITERATIONS (20 IS RECOMENNDED) C IPRED: TIME-STEPS FOR GRAVITY FORCE ONLY C READ(5,100) NSTEP,NOUTP,NREQD,NREQS,NACCE, * IFUNC,IFIXD,MITER,IPRED WRITE(6,200) NSTEP,NOUTP,NREQD,NREQS,NACCE, * IFUNC,IFIXD,MITER,IPRED C DTIME: TIME-STEP INCREMENT C DTEND: END-TIME C DTREC: TIME-STEP AT ACCERALATION RECORD C DELTA, GAMMA: CONSTANTS IN NEWMARK METHOD C DELTA=0.25, GAMMA=0.5 C AZERO, BZERO, OMEGA C FUNCTS=AZERO+BZERO*SIN(OMEGA*STEP*DTIME) C TOLER: CONVERGENCE CONSTANT (TOLER=0.0001) C READ(5,101) DTIME,DTEND,DTREC,DELTA,GAAMA, * AZERO,BZERO,OMEGA,TOLER WRITE(6,201) DTIME,DTEND,DTREC,DELTA,GAAMA, * AZERO,BZERO,OMEGA,TOLER C *** SELECTED TIME-STEPS FOR OUTPUT READ(5,100) (LOUTD(N),N=1,NOUTP) WRITE(6,206) (LOUTD(N),N=1,NOUTP) DO 2 N=1,NSTEP 2 NOUTD(N)=0 DO 4 L=1,NOUTP N=LOUTD(L) 4 NOUTD(N)=1 C *** SELECTED NODES AND ELEMENTS FOR OUTPUT WRITE(6,202) READ(5,100) (IREQD(I),I=1,NREQD) WRITE(6,203) (IREQD(I),I=1,NREQD) READ(5,100) (IREQS(I),I=1,NREQS) WRITE(6,204) (IREQS(I),I=1,NREQS) IF(IFUNC.NE.0) GO TO 50 C *** READ ACCELEROGRAM DATA C IFIXD=0: H & V AC, 1: V, 2: H C AFACT=DTREC/DTIME IF(IFIXD-1) 40,42,44 40 WRITE(6,209) DTREC READ(5,103) (ACCEH(N),N=1,NACCE) WRITE(6,103) (ACCEH(N),N=1,NACCE) WRITE(6,210) DTREC READ(5,103) (ACCEV(N),N=1,NACCE) WRITE(6,103) (ACCEV(N),N=1,NACCE) GO TO 50 42 WRITE(6,210) DTREC READ(5,103) (ACCEV(N),N=1,NACCE) WRITE(6,103) (ACCEV(N),N=1,NACCE) GO TO 50 44 WRITE(6,209) DTREC READ(5,103) (ACCEH(N),N=1,NACCE) WRITE(6,103) (ACCEH(N),N=1,NACCE) 50 CONTINUE 80 CONTINUE C 100 FORMAT(20I5) 101 FORMAT(10F10.4) 102 FORMAT(I5,2F10.5) 103 FORMAT(10F8.4) 104 FORMAT(2I5,F10.5) 200 FORMAT(/'TIME STEPPING PARAMETERS'/ * /2X,'NSTEP=',I5,12X,'NOUTP=',I5,12X,'NREQD=',I5 * /2X,'NREQS=',I5,12X,'NACCE=',I5,12X,'IFUNC=',I5 * /2X,'IFIXD=',I5,12X,'MITER=',I5,12X,'IPRED=',I5) 201 FORMAT(/2X,'DTIME=',G12.4,5X,'DTEND=',G12.4,5X,'DTREC=',G12.4 * /2X,'DELTA=',G12.4,5X,'GAAMA=',G12.4,5X,'AZERO=',G12.4 * /2X,'BZERO=',G12.4,5X,'OMEGA=',G12.4,5X,'TOLER=',G12.4) 202 FORMAT(/'SELECTIVE OUTPUT' ) 203 FORMAT(/2X,'NODE',10I5) 204 FORMAT(/2X,'ELEM',10I5) 206 FORMAT(/'SELECTIVE TIME-STEPS',30(/10I5)) 209 FORMAT(/'HORIZONTAL ACCELERATION ORDINATES AT',F9.4,2X,'SEC'/) 210 FORMAT(/'VERTICAL ACCELERATION ORDINATES AT',F9.4,2X,'SEC'/) RETURN END C * * * * * * * * * * * * SUBROUTINE PREVOS C *** INITIAL STRESSES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LA11/IREQD(10),IREQS(10),IWA(1260) COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA17/VELOI(1850),VELOL(1850),VELOT(1850) COMMON/LA19/DAMPI(1850,1850),RESID(1850),STRIN(1260,3) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) DIMENSION STREI(10,6) C DO 5 IT=1,NSIZE DISPI(IT)=0.D0 DISPT(IT)=0.D0 DISPL(IT)=0.D0 VELOI(IT)=0.D0 5 CONTINUE DO 10 IE=1,NELEM DO 12 I=1,6 STRAG(IE,I)=0.D0 STRSG(IE,I)=0.D0 12 CONTINUE 10 CONTINUE DO 14 IE=1,NELEM IYIEL(IE)=0 IYIEP(IE)=0 IYIEH(IE)=0 IUNLO(IE)=0 14 CONTINUE C *** INPUT INITIAL STRESS GROUPS WRITE(6,200) IF(NIST.LE.0) GO TO 26 DO 20 IS=1,NIST READ(5,100) IS1,(STREI(IS,I),I=1,6) WRITE(6,100) IS1,(STREI(IS,I),I=1,6) 20 CONTINUE DO 22 IE=1,NELEM IS=IST(IE) DO 24 I=1,6 24 STRIN(IE,I)=STREI(IS,I) 22 CONTINUE 26 CONTINUE C *** INPUT INITIAL STRESSES IN EACH ELEMENT IF(NPREV.EQ.0) GO TO 80 OPEN(7,FILE='DASTRE',STATUS='UNKNOWN') DO 30 I=1,3 READ(7,101) (STRIN(IE,I),IE=1,NELEM) WRITE(6,101) (STRIN(IE,I),IE=1,NELEM) 30 CONTINUE CLOSE(7) 80 CONTINUE 100 FORMAT(I2,6F8.3) 101 FORMAT(10F8.3) 200 FORMAT(/'INITIAL STRESSES') RETURN END C * * * * * * * * * * * * SUBROUTINE LOADPL C *** STANDARD LOAD ROUTINE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) DIMENSION ELOAD(8) C DO 2 IV=1,NSIZE 2 FORCE(IV)=0.D0 C *** LOAD DATA C NFL: NUMBER OF NODAL POINT LOADS C IGRAV=0: NOT CONSIDER GRAVITY LOADS, 1:CONSIDER C READ (5,100) NFL,IGRAV WRITE(6,200) NFL,IGRAV C *** READ NODAL POINT LOADS C NOD: NODE NO. C IXY: DIRECTION OF LOAD (1=HORIZONTAL, 2=VERTICAL, 3=MOMENT) C FLO: LOAD C IF(NFL.EQ.0) GO TO 12 WRITE(6,201) DO 10 IL=1,NFL READ(5,101) NOD,IXY,FLO WRITE(6,101) NOD,IXY,FLO IT=(NOD-1)*3+IXY IT=LY(IT) IF(IT.EQ.NN) GO TO 10 FORCE(IT)=FORCE(IT)+FLO 10 CONTINUE 12 CONTINUE C *** GIVE GRAVITY LOADS IF(IGRAV.EQ.0) GO TO 22 GRAVZ=-9.81D0 DO 20 IE=1,NELEM KOL2=K2(IE) GO TO (31,32,31,34,34,34),KOL2 31 CALL TRUSS(IE) NNODE=2 NEVAB=6 GO TO 30 32 CALL BEAM(IE) NNODE=2 NEVAB=6 GO TO 30 34 CALL PLANES(IE,KOL2) NNODE=4 NEVAB=8 MM=IJK(IE,4) IF(MM.LE.0) NNODE=3 IF(MM.LE.0) NEVAB=6 30 CONTINUE C MT=MAT(IE) RHO= PRP(MT,4) IF(RHO.LE.0.D0) RHO=0.D0 DMASS=AES(IE)*RHO*GRAVZ/DFLOAT(NNODE) NEVA2=NEVAB/2 IF(KOL2.EQ.2) NEVA2=2 DO 36 I=1,NEVA2 I2=I*2 IF(KOL2.EQ.2) I2=3*I-1 IT=LLL(IE,I2) IF(IT.EQ.NN) GO TO 36 FORCE(IT)=FORCE(IT)+DMASS 36 CONTINUE 20 CONTINUE 22 CONTINUE WRITE(6,202) WRITE(6,203) (FORCE(IV),IV=1,NSIZE) C 100 FORMAT(10I5) 101 FORMAT(2I5,F10.2) 200 FORMAT(/'NODAL LOADS=',I5,' GRAVITY LOADS=',I5/) 201 FORMAT(/'NODAL LOAD'/' NODE X-Y LOAD') 202 FORMAT(/'FORCE') 203 FORMAT(10E8.2) 204 FORMAT(/'FORCE PORE-WATER PRESSURE') 80 RETURN END C * * * * * * * * * * * * SUBROUTINE LUMASS C *** CALCULATES LUMPED MASS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LA11/IREQD(10),IREQS(10),IWA(1260) COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA14/XMASS(1850,1850),STIFF(1850,1850),STIFS(1850,1850) COMMON/LA19/DAMPI(1850,1850),RESID(1850),STRIN(1260,3) DIMENSION EMASS(8,8) C DO 5 IV=1,NSIZE DO 5 JV=1,NSIZE XMASS(IV,JV)=0.D0 5 DAMPI(IV,JV)=0.D0 C DO 100 IE=1,NELEM KOL2=K2(IE) GO TO (31,32,31,34,34,34),KOL2 31 CALL TRUSS(IE) NNODE=2 NEVAB=4 GO TO 30 32 CALL BEAM(IE) NNODE=2 NEVAB=6 GO TO 30 34 CALL PLANES(IE,KOL2) NNODE=4 NEVAB=8 MM=IJK(IE,4) IF(MM.LE.0) NNODE=3 IF(MM.LE.0) NEVAB=6 30 CONTINUE C MT=MAT(IE) RHO= PRP(MT,4) IF(RHO.LE.0.D0) RHO=0.01D0 ALFA=PRP(MT,8) DO 12 I=1,NEVAB DO 12 J=1,NEVAB 12 EMASS(I,J)=0.D0 DMASS=AES(IE)*RHO/FLOAT(NNODE) DO 14 I=1,NEVAB 14 EMASS(I,I)=DMASS C DO 20 I=1,NEVAB IT=LLL(IE,I) IF(IT.EQ.NN) GO TO 20 DO 22 J=1,NEVAB JT=LLL(IE,J) IF(JT.EQ.NN) GO TO 22 XMASS(IT,JT)=XMASS(IT,JT)+EMASS(I,J) DAMPI(IT,JT)=DAMPI(IT,JT)+EMASS(I,J)*ALFA 22 CONTINUE 20 CONTINUE 100 CONTINUE RETURN END C * * * * * * * * * * * * SUBROUTINE GSTIFF C *** EVALUATES STIFFNESS MATRIX IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA14/XMASS(1850,1850),STIFF(1850,1850),STIFS(1850,1850) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA19/DAMPI(1850,1850),RESID(1850),STRIN(1260,3) C DO 4 IV=1,NSIZE DO 4 JV=1,NSIZE 4 STIFF(IV,JV)=0.D0 C DO 100 IE=1,NELEM MT=MAT(IE) BETA=PRP(MT,9) KOL2=K2(IE) GO TO (31,32,31,34,34,34),KOL2 31 NEVAB=4 GO TO 30 32 NEVAB=6 GO TO 30 34 NEVAB=8 IF(IJK(IE,4).EQ.0) NEVAB=6 30 CONTINUE DO 20 I=1,NEVAB IT=LLL(IE,I) IF(IT.EQ.NN) GO TO 20 DO 22 J=1,NEVAB JT=LLL(IE,J) IF(JT.EQ.NN) GO TO 22 STIFF(IT,JT)=STIFF(IT,JT)+EKK(IE,I,J) DAMPI(IT,JT)=DAMPI(IT,JT)+EKK(IE,I,J)*BETA 22 CONTINUE 20 CONTINUE 100 CONTINUE C *** CALCULATES K-STAR MATRICES CONSC=1.D0/(DTIME*DTIME*DELTA) CONSD=GAAMA/(DTIME*DELTA) DO 36 IV=1,NSIZE DO 36 JV=1,NSIZE 36 STIFS(IV,JV)=XMASS(IV,JV)*CONSC+DAMPI(IV,JV)*CONSD+STIFF(IV,JV) C CALL MATINV(STIFS,NN,NSIZE,0,DET,IND) C RETURN END C * * * * * * * * * * * * SUBROUTINE IMPEXP(IITER,ISTEP) C *** GENERATES PARTIAL EFFECTIVE LOAD VECTOR IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA12/ACCEH(6600),ACCEV(6600) COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA14/XMASS(1850,1850),STIFF(1850,1850),STIFS(1850,1850) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA17/VELOI(1850),VELOL(1850),VELOT(1850) COMMON/LA18/ACCEI(1850),ACCEL(1850) COMMON/LA19/DAMPI(1850,1850),RESID(1850),STRIN(1260,3) DIMENSION TXMAS(1850,1850),TACCE(1850),ACCEJ(1850),ACCEK(1850) C IF(ISTEP.GT.1.OR.IITER.GT.1) GO TO 100 CONSA=DTIME*DTIME*(0.5D0-DELTA) CONSB=DTIME*(1.D0-GAAMA) CONSC=1.D0/(DTIME*DTIME*DELTA) CONSD=GAAMA/(DTIME*DELTA) C DO 10 IP=1,NPOIN DO 12 ID=1,3 IT=(IP-1)*3+ID IT=LY(IT) IF(IT.EQ.NN) GO TO 12 ACCEI(IT)=1.D0 ACCEL(IT)=0.D0 IF(ID.EQ.1.OR.ID.EQ.3) GO TO 12 ACCEI(IT)=0.D0 ACCEL(IT)=1.D0 12 CONTINUE 10 CONTINUE C *** CALCULATES VECTORS FOR HORIZONTAL AND VERTICAL EXCITATION CALL MULTPY(ACCEK,XMASS,ACCEL) CALL MULTPY(ACCEJ,XMASS,ACCEI) CALL MULTPY(DISPL,STIFF,DISPI) C *** CALCULATES INITIAL ACCELERATION CALL MULTPY(VELOL,DAMPI,VELOI) DO 14 IV=1,NSIZE 14 TACCE(IV)=FORCE(IV)-DISPL(IV)-VELOL(IV) DO 16 IV=1,NSIZE DO 16 JV=1,NSIZE 16 TXMAS(IV,JV)=XMASS(IV,JV) CALL MATINV(TXMAS,NN,NSIZE,0,DET,IND) CALL MULTPY(ACCEI,TXMAS,TACCE) WRITE(6,200) WRITE(6,201) (ACCEI(IV),IV=1,NSIZE) C 100 CONTINUE IF(IITER.GT.1) GO TO 110 C *** CALCULATES PREDICTED DISPLACEMENT AND VELOCITY VECTOR DO 20 IV=1,NSIZE DISPI(IV)=DISPI(IV)+DTIME*VELOI(IV)+CONSA*ACCEI(IV) VELOI(IV)=VELOI(IV)+CONSB*ACCEI(IV) DISPT(IV)=DISPI(IV) VELOT(IV)=VELOI(IV) ACCEI(IV)=CONSC*(DISPT(IV)-DISPI(IV)) 20 CONTINUE C *** CALCULATES LOAD VECTORS FACTS=FUNCTS(ISTEP) FACTH=FUNCTA(ACCEH,AFACT,ISTEP) FACTV=FUNCTA(ACCEV,AFACT,ISTEP) C 110 CONTINUE IF(ISTEP.EQ.1) GO TO 112 C *** CALCULATES PARTIAL EFFECTIVE LOAD VECTOR CALL MULTPY(VELOL,DAMPI,VELOT) 112 DO 40 IV=1,NSIZE IF(IFUNC.NE.0) GO TO 42 IF(IFIXD.EQ.2) * DISPL(IV)=-VELOL(IV)-FACTH*ACCEJ(IV)+FORCE(IV) IF(IFIXD.EQ.1) * DISPL(IV)=-VELOL(IV)-FACTV*ACCEK(IV)+FORCE(IV) IF(IFIXD.EQ.0) * DISPL(IV)=-VELOL(IV)-FACTH*ACCEJ(IV)-FACTV*ACCEK(IV) * +FORCE(IV) IF(IFUNC.EQ.0) GO TO 40 42 CONTINUE DISPL(IV)=-VELOL(IV)+FORCE(IV)*FACTS 40 CONTINUE 200 FORMAT(/'INITIAL ACCELERATION '/) 201 FORMAT(10E8.2) RETURN END C * * * * * * * * * * * * SUBROUTINE RESEPL(IITER,ISTEP) C *** EVALUATES RESIDUAL FORCES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA19/DAMPI(1850,1850),RESID(1850),STRIN(1260,3) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) COMMON/LA22/NTENS,KTENS,ITENS(1260),JTENS(1260) COMMON/LA23/TMX(1260,3,3),TMI(1260,3,3),EFFST(1260) DIMENSION ELOAD(8),TU(2220),EU(6),DSTRE(4), * STRAN(6),STRES(6),SIGMA(6),SGTOT(6),DESIG(6) C DO 2 IV=1,NSIZE RESID(IV)=0.D0 TU(IV)=DISPT(IV) 2 CONTINUE TU(NN)=0.D0 IYIEL(IE)=0 PI1(IE)=0.D0 PI2(IE)=0.D0 C DO 10 IE=1,NELEM KOL2=K2(IE) MT=MAT(IE) YOUNG=PRP(MT,1) POISS=PRP(MT,2) THICK=PRP(MT,3) COHES=PRP(MT,5) FRICT=PRP(MT,6)*0.017453292D0 DAILT=PRP(MT,7)*0.017453292D0 UNIAX=COHES*DCOS(FRICT) PREYS=UNIAX C DO 8 I=1,6 8 STRES(I)=0.D0 GO TO (11,12,11,14,14,14),KOL2 11 NEVAB=4 NSTRE=1 GO TO 20 12 NEVAB=6 NSTRE=6 GO TO 20 14 NEVAB=8 NSTRE=3 IF(IJK(IE,4).EQ.0) NEVAB=6 20 CONTINUE C CALL LINGNL(IE,KOL2,TU,YOUNG,POISS,STRAN,STRES) DO 22 I=1,NSTRE 22 STRAG(IE,I)=STRAG(IE,I)+STRAN(I) IF(ISTEP.GT.1.OR.IITER.GT.1) GO TO 23 DO 24 I=1,NSTRE 24 STRES(I)=STRES(I)+STRIN(IE,I) 23 CONTINUE DO 25 I=1,NSTRE DESIG(I)=STRES(I) SIGMA(I)=STRSG(IE,I)+STRES(I) 25 CONTINUE IF(KOL2.LE.4) GO TO 30 C IF(ISTEP.LE.IPRED) GO TO 27 DO 28 I=1,2 28 SIGMA(I)=SIGMA(I)-PW(IE,I) 27 CONTINUE C CALL YIELDJ(MT,KOL2,SIGMA,YIELD) ESPRE=EFFST(IE)-PREYS C ESPRE.GE.0: THE ELEMENT YIELDED AT PRECEDING STEP IF(ESPRE.GE.0.D0) GO TO 32 ESCUR=YIELD-PREYS C ESCUR.GT.0: THE ELEMENT HAS YIELDED AT PRESENT STEP IF(ESCUR.LE.0.D0) IYIEL(IE)=0 IF(ESCUR.LE.0.D0) GO TO 38 RFACT=ESCUR/(YIELD-EFFST(IE)) GO TO 34 32 CONTINUE ESCUR=YIELD-EFFST(IE) C ESCUR.LE.0: THE ELEMENT DOES NOT YIELD AT PRESENT STEP IF(ESCUR.LE.0.D0) IYIEL(IE)=0 IF(ESCUR.LE.0.D0) GO TO 38 RFACT=1.D0 34 CONTINUE IYIEL(IE)=1 C MSTEP=ESCUR*8.D0/UNIAX+1.D0 IF(MSTEP.GT.10) MSTEP=10 ASTEP=MSTEP REDUC=1.D0-RFACT C C CALL CALSTA(IE,SIGMA,ISTEP) DO 36 I=1,NSTRE SGTOT(I)=STRSG(IE,I)+REDUC*STRES(I) IF(ISTEP.GT.IPRED) SGTOT(I)=SGTOT(I)-PW(IE,I) STRAN(I)=RFACT*STRAN(I)/ASTEP 36 CONTINUE C DO 40 JSTEP=1,MSTEP CALL DEPMAT(IE,MT,KOL2,SGTOT) DO 44 I=1,NSTRE DSTRE(I)=0.D0 DO 44 J=1,NSTRE 44 DSTRE(I)=DSTRE(I)+DEP(IE,I,J)*STRAN(J) DO 45 I=1,NSTRE 45 SGTOT(I)=SGTOT(I)+DSTRE(I) 40 CONTINUE CALL YIELDJ(MT,KOL2,SGTOT,YIELD) CURYS=UNIAX BRING=1.D0 IF(YIELD.GT.CURYS) BRING=CURYS/YIELD DO 46 I=1,NSTRE 46 STRSG(IE,I)=BRING*SGTOT(I) EFFST(IE)=BRING*YIELD IF(ISTEP.LE.IPRED) GO TO 47 DO 48 I=1,2 48 STRSG(IE,I)=STRSG(IE,I)+PW(IE,I) 47 CONTINUE C *** ALTERNATIVE LOCATION OF STRESS REDUCTION LOOP TERMINATION CARD GO TO 50 38 CONTINUE DO 54 I=1,NSTRE 54 STRSG(IE,I)=SIGMA(I) IF(ISTEP.LE.IPRED) GO TO 53 DO 55 I=1,2 55 STRSG(IE,I)=STRSG(IE,I)+PW(IE,I) 53 CONTINUE GO TO 52 C 30 CONTINUE DO 56 I=1,NSTRE 56 STRSG(IE,I)=STRSG(IE,I)+DESIG(I) 52 CONTINUE EFFST(IE)=YIELD 50 CONTINUE C *** NO-TENSION CALL NOTENS(IE,KOL2) C C *** CALCULATE THE EQUIVALENT NODAL FORCES DVOLU=AES(IE) NSTR1=NSTRE IF(KOL2.EQ.2) GO TO 60 IF(KOL2.GE.4) NSTR1=3 DO 62 I=1,NEVAB ELOAD(I)=0.D0 DO 62 J=1,NSTR1 62 ELOAD(I)=ELOAD(I)+BMX(IE,J,I)*STRSG(IE,J)*DVOLU GO TO 64 60 CONTINUE DO 66 I=1,NEVAB IT=LLL(IE,I) 66 EU(I)=TU(IT) DO 68 I=1,NEVAB ELOAD(I)=0.D0 DO 68 J=1,NEVAB 68 ELOAD(I)=ELOAD(I)+EKK(IE,I,J)*EU(J) 64 CONTINUE DO 70 I=1,NEVAB IT=LLL(IE,I) IF(IT.EQ.NN) GO TO 70 RESID(IT)=RESID(IT)+ELOAD(I) 70 CONTINUE 10 CONTINUE RETURN END C * * * * * * * * * * SUBROUTINE ITRATE(IITER) C *** CALCULATES INCREMENT IN DISPLACEMENT AND APPLIES CONVERGENCE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT COMMON/LA14/XMASS(1850,1850),STIFF(1850,1850),STIFS(1850,1850) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA17/VELOI(1850),VELOL(1850),VELOT(1850) COMMON/LA18/ACCEI(1850),ACCEL(1850) COMMON/LA19/DAMPI(1850,1850),RESID(1850),STRIN(1260,3) DIMENSION TVECT(1850),DISPD(1850) C CONSD=DTIME*(1.D0-GAAMA) CONSF=1.D0/(DTIME*DTIME*DELTA) NCHEK=0 CALL MULTPY(ACCEL,XMASS,ACCEI) C *** CALCULATES TOTAL EFFECTIVE LOAD VECTOR DO 10 IV=1,NSIZE 10 TVECT(IV)=DISPL(IV)-ACCEL(IV)-RESID(IV) C *** CALCULATES DELTA DISPLACEMENT CALL MULTPY(DISPD,STIFS,TVECT) C *** APPLIES CONVERGENCE SUMPP=0.D0 SUMPQ=0.D0 DO 12 IV=1,NSIZE DISPP=DISPD(IV) DISPQ=DISPT(IV)+DISPP DISPT(IV)=DISPQ SUMPP=SUMPP+DISPP*DISPP SUMPQ=SUMPQ+DISPQ*DISPQ 12 CONTINUE DO 14 IV=1,NSIZE ACCEI(IV)=CONSF*(DISPT(IV)-DISPI(IV)) 14 VELOT(IV)=VELOI(IV)+CONSD*ACCEI(IV) SUMPP=DSQRT(SUMPP/SUMPQ) WRITE(6,200) IITER,SUMPP,SUMPQ IF(SUMPP.GT.TOLER) GO TO 30 NCHEK=1 GO TO 20 30 IF(IITER.LT.MITER) GO TO 32 20 DO 22 IV=1,NSIZE VELOI(IV)=VELOT(IV) 22 DISPI(IV)=DISPT(IV) 32 CONTINUE 200 FORMAT(I5,2E10.3) RETURN END C * * * * * * * * * * * * SUBROUTINE LINGNL(IE,KOL2,TU,YOUNG,POISS,STRAN,STRES) C *** ELASTIC STRAINS AND STRESSES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) DIMENSION EU(8),TU(1850),STRAN(6),STRES(6) C DO 5 I=1,6 STRAN(I)=0.D0 5 STRES(I)=0.D0 GO TO (11,12,11,14,14,14),KOL2 C *** TRUSS & TEXTILE 11 CONTINUE DO 20 I=1,4 IT=LLL(IE,I) 20 EU(I)=TU(IT) STRAN(1)=0.D0 DO 22 I=1,4 22 STRAN(1)=STRAN(1)+BMX(IE,1,I)*EU(I) STRAN(1)=STRAN(1)-STRAG(IE,1) STRES(1)=STRAN(1)*YOUNG GO TO 10 C *** BEAM 12 CONTINUE DO 30 I=1,6 IT=LLL(IE,I) 30 EU(I)=TU(IT) DO 32 I=1,6 STRES(I)=0.D0 DO 32 J=1,6 32 STRES(I)=STRES(I)+DBM(IE,I,J)*EU(J) DO 34 I=1,6 STRAN(I)=0.D0 34 STRES(I)=STRES(I)-STRSG(IE,I) C *** STRES(I) MEAN PARALLEL & NORMAL NODAL FORCES C AND BENDING MOMENT AT EACH NODE GO TO 10 C *** SOLID AND INTERFACE 14 CONTINUE NEVAB=8 IF(IJK(IE,4).LE.0) NEVAB=6 DO 40 I=1,NEVAB IT=LLL(IE,I) 40 EU(I)=TU(IT) DO 42 I=1,3 STRAN(I)=0.D0 DO 42 J=1,NEVAB 42 STRAN(I)=STRAN(I)+BMX(IE,I,J)*EU(J) DO 44 I=1,3 44 STRAN(I)=STRAN(I)-STRAG(IE,I) DO 46 I=1,3 STRES(I)=0.D0 DO 46 J=1,3 46 STRES(I)=STRES(I)+DMX(IE,I,J)*STRAN(J) 10 RETURN END C * * * * * * * * * * * * SUBROUTINE YIELDJ(MT,KOL2,STEMP,YIELD) C *** STRESS INVARIANTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) DIMENSION STEMP(4) C PHIRA=PRP(MT,6)*0.017453292 SNPHI=DSIN(PHIRA) COPHI=DCOS(PHIRA) GO TO (10,10,10,10,15,16),KOL2 15 CONTINUE C *** MOHR-COULOMB P1=(STEMP(1)+STEMP(2))*0.5D0 P2=(STEMP(1)-STEMP(2))*0.5D0 P3=P2*P2+STEMP(3)**2 P4=0.D0 IF(P3.GT.0.D0) P4=DSQRT(P3) YIELD=P4-P1*SNPHI GO TO 10 16 CONTINUE C *** COULOMB TAU=STEMP(3) TAU=DABS(TAU) SIG=STEMP(2) IF(SIG.LT.0.D0) SIG=0.D0 YIELD=TAU*COPHI-SIG*SNPHI 10 CONTINUE RETURN END C * * * * * * * * * * * * SUBROUTINE CALSTA(IE,SIGMA,ISTEP) C *** CALCULATES YIELD STRESSES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) DIMENSION SIGMA(3),STD(3),STP(3) C MT=MAT(IE) FAI=PRP(MT,6) FAI=FAI*3.14159D0/180.D0 DO 2 I=1,3 STP(I)=STRSG(IE,I) IF(ISTEP.LE.IPRED) GO TO 2 STP(I)=STRSG(IE,I)-PW(IE,I) 2 CONTINUE DO 4 I=1,3 4 STD(I)=SIGMA(I)-STP(I) KOL2=K2(IE) GO TO (10,10,10,10,15,16),KOL2 15 CONTINUE SX=SIGMA(1) SY=SIGMA(2) TA=SIGMA(3) CALL MOHRCO(IE,SX,SY,TA,F1,S,P2,B0) SX=STP(1) SY=STP(2) TA=STP(3) CALL MOHRCO(IE,SX,SY,TA,F0,S,P2,B0) C1=-F0/(F1-F0) SX=STP(1)+STD(1)*C1 SY=STP(2)+STD(2)*C1 TA=STP(3)+STD(3)*C1 CALL MOHRCO(IE,SX,SY,TA,F2,S,P2,B0) B2=B0**(-0.5) SF=DSIN(FAI) A1=B2*P2-SF A2=B2*P2*(-1.D0)-SF A3=B2*4.D0*TA R=A1*STD(1)+A2*STD(2)+A3*STD(3) IF(R.LE.0.D0) R=0.0001D0 C=C1-F2/R IF(C.GT.1.D0) C=C1 GO TO 18 16 SG=SIGMA(2) TA=SIGMA(3) CALL COULOM(IE,SG,TA,F1,S,1) SG=STP(2) TA=STP(3) CALL COULOM(IE,SG,TA,F0,S,1) C=-F0/(F1-F0) 18 CONTINUE DO 22 I=1,3 22 STAX(IE,I)=STP(I)+STD(I)*C 10 CONTINUE RETURN END C * * * * * * * * * * * * SUBROUTINE DEPMAT(IE,MT,KOL2,SIGMA) C *** CALCULATES DEP MATRIX IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) COMMON/LA23/TMX(1260,3,3),TMI(1260,3,3),EFFST(1260) DIMENSION T(3,3),TD(3,3),DPF(3,3),SIGMA(6) C DO 2 I=1,3 2 STAX(IE,I)=SIGMA(I) GO TO (10,10,10,10,15,16),KOL2 15 CONTINUE SX=SIGMA(1) SY=SIGMA(2) TA=SIGMA(3) CALL MOHRCO(IE,SX,SY,TA,F,S,P2,B0) CALL PRINCE(SX,SY,TA,S1,S3,AG) PI3(IE)=AG C E=PRP(MT,1) P=PRP(MT,2) G=E/(2.D0*(1.D0+P)) FAI=PRP(MT,6)*3.14159D0/180.D0 DLT=PRP(MT,7)*3.14159D0/180.D0 TFI=DSIN(FAI)/DCOS(FAI) TDL=DSIN(DLT)/DCOS(DLT) C SI=1.D0 IF(LRE(IE).GE.2) SI=-1.D0 SIT=PI3(IE) ALF=3.14159D0*0.25D0+DATAN(TFI)*0.5D0 BET=(ALF+SIT)*(-1.D0) IF(LRE(IE).GE.2) BET=ALF-SIT PI4(IE)=BET IF(ICRI.LE.0) GO TO 20 C *** PERFECTLY PLASTIC TFI=0.D0 TDL=0.D0 20 CONTINUE CB=DCOS(BET) SB=DSIN(BET) T(1,1)=CB*CB T(1,2)=SB*SB T(1,3)=-2.D0*SB*CB T(2,1)=SB*SB T(2,2)=CB*CB T(2,3)=2.D0*SB*CB T(3,1)=SB*CB T(3,2)=-SB*CB T(3,3)=CB*CB-SB*SB DO 22 I=1,3 DO 22 J=1,3 22 TMX(IE,I,J)=T(I,J) CALL DINV(T,3,3,0,DET,IND) DO 24 I=1,3 DO 24 J=1,3 24 TMI(IE,I,J)=T(I,J) DO 26 I=1,3 DO 26 J=1,3 26 T(I,J)=TMX(IE,I,J) C C1=E*(1.D0-P)/((1.D0+P)*(1.D0-2.D0*P)) C2=E*P/((1.D0+P)*(1.D0-2.D0*P)) B1=1.D0/(C1*TFI*TDL+G) DPF(1,1)=C1 DPF(1,2)=C2 DPF(1,3)=0.D0 DPF(2,1)=C2 DPF(2,2)=C1 DPF(2,3)=0.D0 DPF(3,1)=-SI*C2*TFI DPF(3,2)=-SI*C1*TFI DPF(3,3)=0.D0 CALL MULT(T,DPF,TD,3,3,3) CALL XULT(TD,T,DPF,3,3,3) GO TO 18 16 CONTINUE G=PRP(MT,1) E=PRP(MT,2) P=PRP(MT,3) FAI=PRP(MT,6)*3.14159D0/180.D0 DLT=PRP(MT,7)*3.14159D0/180.D0 TFI=DSIN(FAI)/DCOS(FAI) TDL=DSIN(DLT)/DCOS(DLT) SI=-1.D0 IF(SIGMA(3).LT.0.) SI=1.D0 IF(ICRI.LE.0) GO TO 28 C *** PERFECTLY PLASTIC TFI=0.D0 TDL=0.D0 28 CONTINUE C1=E*(1.-P)/((1.+P)*(1.-2.*P)) C2=E*P/((1.+P)*(1.-2.*P)) B1=1./(C1*TFI*TDL+G) DPF(1,1)=C1 DPF(1,2)=C2 DPF(1,3)=0.D0 DPF(2,1)=DPF(1,2) DPF(2,2)=C1 DPF(2,3)=0.D0 DPF(3,1)=-SI*C2*TFI DPF(3,2)=-SI*C1*TFI DPF(3,3)=0.D0 18 CONTINUE DO 30 I=1,3 DO 30 J=1,3 30 DEP(IE,I,J)=DPF(I,J) 10 CONTINUE RETURN END C * * * * * * * * * * * * SUBROUTINE NOTENS(IE,KOL2) C *** FINDS NO-TENSION ELEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA22/NTENS,KTENS,ITENS(1260),JTENS(1260) DIMENSION SIG(6) C DO 2 I=1,3 2 SIG(I)=STRSG(IE,I) JTENS(IE)=0 MT=MAT(IE) TENST=-PRP(MT,10) GO TO (10,10,11,10,15,16),KOL2 11 CONTINUE SG=SIG(1) IF(SG.LE.0.D0) GO TO 18 SIG(1)=0.D0 C KTENS=KTENS+1 C ITENS(KTENS)=IE JTENS(IE)=1 GO TO 18 15 CONTINUE SX=SIG(1) SY=SIG(2) TA=SIG(3) CALL PRINCE(SX,SY,TA,S1,S3,TH) IF(S3.GE.TENST) GO TO 18 C KTENS=KTENS+1 C ITENS(KTENS)=IE JTENS(IE)=1 C=1.D0 IF(TA.LT.0.D0) C=-1.D0 C IF(S1.GT.TENST) GO TO 20 SIG(1)=TENST SIG(2)=TENST SIG(3)=0.D0 GO TO 18 20 CONTINUE IF(SY.GE.SX) GO TO 30 C SX.GT.SY IF(SX) 50,50,40 40 CONTINUE SY=TENST+(S1-SX) SIG(2)=SY SIG(1)=SX GO TO 35 30 CONTINUE C SY.GE.SX IF(SY) 52,52,42 42 CONTINUE SX=TENST+(S1-SY) SIG(1)=SX SIG(2)=SY 35 CONTINUE A=((S1-TENST)*0.5D0)**2-((SX-SY)*0.5D0)**2 IF(A.LT.0.000001D0) A=0.000001D0 SIG(3)=C*DSQRT(A) GO TO 18 50 C1=1.D0 GO TO 55 52 C1=-1.D0 55 CONTINUE R=(S1-S3)*0.5D0 B2=R*R-TA*TA IF(B2.GT.0.D0) GO TO 57 T=0.5D0*3.14159D0 GO TO 58 57 CONTINUE B=DSQRT(B2) AT=TA/B T=DATAN(AT) 58 CONTINUE S=(S1-TENST)*0.5D0*DCOS(T) SIG(1)=(S1+TENST)*0.5D0+C1*S SIG(2)=(S1+TENST)*0.5D0-C1*S SIG(3)=(S1-TENST)*0.5D0*DSIN(T)*C GO TO 18 16 CONTINUE SG=SIG(2) IF(SG.GE.0.D0) GO TO 18 KTENS=KTENS+1 ITENS(KTENS)=IE JTENS(IE)=KTENS SIG(2)=0.001D0 TA=SIG(3) C=1.D0 IF(TA.LT.0.D0) C=-1.D0 MT=MAT(IE) SIG(3)=C*PRP(MT,5) 18 CONTINUE 10 CONTINUE DO 4 I=1,3 4 STRSG(IE,I)=SIG(I) RETURN END C * * * * * * * * * * * * SUBROUTINE MOHRCO(IE,SX,SY,TA,F,S,P2,B0) C *** MOHR-COULOMB YIELD CONDITION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) MT=MAT(IE) CC=PRP(MT,5) FAI=PRP(MT,6) FAI=FAI*3.14159D0/180.D0 SFI=DSIN(FAI) CFI=DCOS(FAI) P1=SX+SY P2=SX-SY A1=P1*SFI+2.D0*CC*CFI B0=P2*P2+4.D0*TA*TA IF(B0.LT.0.0001D0) B0=0.0001D0 B1=DSQRT(B0) F=B1-A1 S=A1/B1 RETURN END C * * * * * * * * * * * * SUBROUTINE COULOM(IE,SG,TA,F,S,IC) C *** COULOMB YIELD CONDITION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) MT=MAT(IE) CC=PRP(MT,5) FAI=PRP(MT,6) FAI=FAI*3.14159D0/180.D0 SFI=DSIN(FAI) CFI=DCOS(FAI) ST=CC+SG*SFI/CFI IF(IC.LE.0) ST=VFUNC(ST) TA=DABS(TA) F=TA-ST S=0.D0 IF(TA.GT.0.001D0) S=ST/TA RETURN END C * * * * * * * * * * * * SUBROUTINE TRUSS(IE) C *** TRUSS ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) DIMENSION LL(4),BM(4),EK(4,4),T(4,4),TEK(4,4),TET(4,4) C II=IJK(IE,1) JJ=IJK(IE,2) C MT=MAT(IE) E=PRP(MT,1) A=PRP(MT,2) XI=XX(II) YI=YY(II) XJ=XX(JJ) YJ=YY(JJ) DX=XJ-XI DY=YJ-YI EL=DSQRT(DX*DX+DY*DY) C=DX/EL S=DY/EL BM(1)=1.D0/EL BM(2)=0.D0 BM(3)=-1.D0/EL BM(4)=0.D0 DO 10 I=1,4 DO 10 J=1,4 10 EK(I,J)=BM(I)*BM(J)*E*EL*A C DO 20 I=1,4 DO 20 J=1,4 20 T(I,J)=0.D0 T(1,1)=C T(1,2)=S T(2,1)=-S T(2,2)=C T(3,3)=C T(3,4)=S T(4,3)=-S T(4,4)=C CALL WULT(T,EK,TEK,4,4,4) CALL MULT(TEK,T,TET,4,4,4) DO 22 I=1,4 DO 22 J=1,4 22 EKK(IE,I,J)=TET(I,J) AES(IE)=A*EL DO 24 J=1,4 BMX(IE,1,J)=0.D0 DO 24 K=1,4 24 BMX(IE,1,J)=BMX(IE,1,J)+BM(K)*T(K,J) LL(4)=3*JJ-1 LL(3)=LL(4)-1 LL(2)=3*II-1 LL(1)=LL(2)-1 DO 30 I=1,4 IT=LL(I) 30 LLL(IE,I)=LY(IT) RETURN END C * * * * * * * * * * * * SUBROUTINE BEAM(IE) C *** BEAM ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) DIMENSION LL(10),EK(6,6),T(6,6),TEK(6,6),TET(6,6),EKT(6,6) C II=IJK(IE,1) JJ=IJK(IE,2) C MT=MAT(IE) E= PRP(MT,1) A= PRP(MT,2) AI=PRP(MT,3) XI=XX(II) YI=YY(II) XJ=XX(JJ) YJ=YY(JJ) DX=XJ-XI DY=YJ-YI EL=DSQRT(DX*DX+DY*DY) C=DX/EL S=DY/EL G=E*A/EL G5=2.D0*E*AI/EL G4=2.D0*G5 G3=3.D0*G5/EL G2=2.D0*G3/EL C DO 20 I=1,6 DO 20 J=1,6 T(I,J)=0.D0 EK(I,J)=0.D0 20 CONTINUE EK(1,1)=G EK(2,2)=G2 EK(3,3)=G4 EK(4,4)=G EK(5,5)=G2 EK(6,6)=G4 EK(1,4)=-G EK(2,3)=G3 EK(2,5)=-G2 EK(2,6)=G3 EK(3,5)=-G3 EK(3,6)=G5 EK(5,6)=-G3 DO 22 I=1,5 IP1=I+1 DO 24 J=IP1,6 24 EK(J,I)=EK(I,J) 22 CONTINUE DO 26 K=1,4,3 T(K,K)=C T(K,K+1)=S T(K+1,K)=-S T(K+1,K+1)=C T(K+2,K+2)=1.D0 26 CONTINUE CALL WULT(T,EK,TEK,6,6,6) CALL MULT(TEK,T,TET,6,6,6) CALL MULT(EK,T,EKT,6,6,6) DO 28 I=1,6 DO 28 J=1,6 DBM(IE,I,J)=EKT(I,J) EKK(IE,I,J)=TET(I,J) 28 CONTINUE AES(IE)=A*EL LL(6)=3*JJ LL(5)=LL(6)-1 LL(4)=LL(5)-1 LL(3)=3*II LL(2)=LL(3)-1 LL(1)=LL(2)-1 DO 30 I=1,6 IT=LL(I) 30 LLL(IE,I)=LY(IT) RETURN END C * * * * * * * * * * * * SUBROUTINE PLANES(IE,KOL2) C *** PLANE STRAIN & INTERFACE ELEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) DIMENSION LL(10),D(3,3),ES(8,8),IOE(6),B(3,6),DB(3,6), * EK(6,6),BM(3,8),BMT(8),T(8,8),TEK(8,8),XXL(4),YYL(4), * COR(8),COD(8) C II=IJK(IE,1) JJ=IJK(IE,2) KK=IJK(IE,3) MM=IJK(IE,4) MT=MAT(IE) C DO 2 I=1,3 DO 2 J=1,3 2 D(I,J)=0.D0 C E=PRP(MT,1) P=PRP(MT,2) H=PRP(MT,3) GO TO (80,80,80,14,14,16),KOL2 14 CONTINUE GP=1.D0+P GM=1.D0-P GN=1.D0-2.D0*P G=E/(GP*GN) G1=GM*G G2=P*G G3=0.5D0*E/GP D(1,1)=G1 D(2,2)=G1 D(3,3)=G3 D(1,2)=G2 D(2,1)=G2 GO TO 18 16 CONTINUE GP=1.D0+H GM=1.D0-H GN=1.D0-2.D0*H G=P/(GP*GN) G1=GM*G G2=H*G D(1,1)=G1 D(2,2)=G1 D(3,3)=E D(1,2)=G2 D(2,1)=G2 DX=XX(JJ)-XX(II) DY=YY(JJ)-YY(II) EL=DSQRT(DX*DX+DY*DY) C=DX/EL S=DY/EL DO 20 I=1,8 DO 20 J=1,8 20 T(I,J)=0.D0 DO 22 I=1,7,2 T(I,I)=C T(I,I+1)=S T(I+1,I)=-S T(I+1,I+1)=C 22 CONTINUE COR(1)=XX(II) COR(2)=YY(II) COR(3)=XX(JJ) COR(4)=YY(JJ) COR(5)=XX(KK) COR(6)=YY(KK) COR(7)=XX(MM) COR(8)=YY(MM) DO 24 I=1,8 COD(I)=0.D0 DO 24 J=1,8 24 COD(I)=COD(I)+T(I,J)*COR(J) DO 26 I=1,4 IEV=2*I IOD=IEV-1 XXL(I)=COD(IOD) YYL(I)=COD(IEV) 26 CONTINUE C 18 CONTINUE DO 30 I=1,3 DO 30 J=1,6 30 B(I,J)=0.D0 DO 32 I=1,8 DO 32 J=1,8 32 ES(I,J)=0.D0 DO 34 J=1,8 34 BMT(J)=0.D0 AE=0.D0 I1=0 I2=1 I3=2 N4=4 NV=8 IF(MM.EQ.0) N4=1 IF(MM.EQ.0) NV=6 DO 40 LI=1,N4 I1=I1+1 I2=I2+1 IF(I2.EQ.5) I2=1 I3=I3+1 IF(I3.EQ.5) I3=1 I1E=2*I1 I1O=I1E-1 I2E=2*I2 I2O=I2E-1 I3E=2*I3 I3O=I3E-1 IOE(1)=I1O IOE(2)=I1E IOE(3)=I2O IOE(4)=I2E IOE(5)=I3O IOE(6)=I3E GO TO (80,80,80,44,44,46),KOL2 44 CONTINUE IIT=IJK(IE,I1) JJT=IJK(IE,I2) KKT=IJK(IE,I3) XI=XX(IIT) YI=YY(IIT) XJ=XX(JJT) YJ=YY(JJT) XK=XX(KKT) YK=YY(KKT) GO TO 48 46 CONTINUE XI=XXL(I1) YI=YYL(I1) XJ=XXL(I2) YJ=YYL(I2) XK=XXL(I3) YK=YYL(I3) 48 CONTINUE AI=XJ*YK-XK*YJ AJ=XK*YI-XI*YK AK=XI*YJ-XJ*YI BI=YJ-YK BJ=YK-YI BK=YI-YJ CI=XK-XJ CJ=XI-XK CK=XJ-XI DA=(AI+AJ+AK)*0.5D0 AE=AE+DABS(DA)*0.5D0 D2=0.125D0/DA IF(MM.EQ.0) D2=0.5D0/DA BMT(I1E)=BMT(I1E)-BI*D2 BMT(I1O)=BMT(I1O)-CI*D2 BMT(I2E)=BMT(I2E)-BJ*D2 BMT(I2O)=BMT(I2O)-CJ*D2 BMT(I3E)=BMT(I3E)-BK*D2 BMT(I3O)=BMT(I3O)-CK*D2 DEL2=1.D0/(AI+AJ+AK) B(3,1)=-CI*DEL2 B(3,2)=-BI*DEL2 B(3,3)=-CJ*DEL2 B(3,4)=-BJ*DEL2 B(3,5)=-CK*DEL2 B(3,6)=-BK*DEL2 DO 50 I=1,3 IL=2*I IO=IL-1 B(1,IO)=B(3,IL) B(2,IL)=B(3,IO) 50 CONTINUE CALL MULT(D,B,DB,3,3,6) CALL WULT(B,DB,EK,6,3,6) DO 52 I=1,6 IL=IOE(I) DO 54 J=1,6 JE=IOE(J) ES(IL,JE)=ES(IL,JE)+EK(I,J) 54 CONTINUE 52 CONTINUE 40 CONTINUE C IF(MM.EQ.0) AE=DABS(DA) AES(IE)=AE GO TO (80,80,80,64,64,66),KOL2 64 CONTINUE C1=0.25D0 IF(MM.EQ.0) C1=1.D0 DO 70 I=1,NV DO 70 J=1,NV 70 EKK(IE,I,J)=ES(I,J)*AE*C1 DO 72 I=1,3 DO 72 J=1,NV 72 BMX(IE,I,J)=0.D0 DO 74 J=1,NV 74 BMX(IE,3,J)=BMT(J) N1=NV/2 DO 76 J=1,N1 JE=2*J JO=JE-1 BMX(IE,1,JO)=BMX(IE,3,JE) 76 BMX(IE,2,JE)=BMX(IE,3,JO) GO TO 68 66 CONTINUE CALL WULT(T,ES,TEK,8,8,8) CALL MULT(TEK,T,ES,8,8,8) DO 78 I=1,8 DO 78 J=1,8 78 EKK(IE,I,J)=ES(I,J)*AES(IE)*0.25D0 DO 90 I=1,3 DO 90 J=1,8 90 BM(I,J)=0.D0 DO 92 J=1,8 92 BM(3,J)=BMT(J) DO 94 J=1,4 JE=2*J JO=JE-1 BM(1,JO)=BM(3,JE) BM(2,JE)=BM(3,JO) 94 CONTINUE DO 96 I=1,3 DO 96 J=1,8 BMX(IE,I,J)=0.D0 DO 96 K=1,8 96 BMX(IE,I,J)=BMX(IE,I,J)+BM(I,K)*T(K,J) 68 CONTINUE LL(8)=3*MM-1 LL(7)=LL(8)-1 LL(6)=3*KK-1 LL(5)=LL(6)-1 LL(4)=3*JJ-1 LL(3)=LL(4)-1 LL(2)=3*II-1 LL(1)=LL(2)-1 DO 98 I=1,NV IT=LL(I) 98 LLL(IE,I)=LY(IT) DO 100 I=1,3 DO 100 J=1,NV DBM(IE,I,J)=0.D0 DO 100 K=1,3 100 DBM(IE,I,J)=DBM(IE,I,J)+D(I,K)*BMX(IE,K,J) DO 102 I=1,3 DO 102 J=1,3 102 DMX(IE,I,J)=D(I,J) CALL DINV(D,3,3,0,DET,IND) DO 104 I=1,3 DO 104 J=1,3 104 DMI(IE,I,J)=D(I,J) 80 CONTINUE RETURN END C * * * * * * * * * * * * FUNCTION FUNCTS(JSTEP) C *** HEAVISIDE AND HARMONIC TIME FUNCTION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA10/AZERO,BZERO,OMEGA,TOLER,AFACT C FUNCTS=0.D0 IF(IFUNC.EQ.0.OR.IFUNC.GE.3) GO TO 80 IF(JSTEP.EQ.0.OR.FLOAT(JSTEP)*DTIME.GT.DTEND) GO TO 80 IF(IFUNC.EQ.1) FUNCTS=1.D0 IF(IFUNC.EQ.2) ARGUM=OMEGA*JSTEP*DTIME IF(IFUNC.EQ.2) FUNCTS=AZERO+BZERO*DSIN(ARGUM) 80 RETURN END C * * * * * * * * * * * * FUNCTION FUNCTA(ACC,AFACT,IST) C *** ACCELEROGRAM INTERPOLATION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA DIMENSION ACC(2000) C FUNCTA=0.D0 IF(IFUNC.NE.0) RETURN TE=FLOAT(IST)*DTIME IF(IST.EQ.0.OR.TE.GT.DTEND) RETURN X=(FLOAT(IST)-1.D0)/AFACT+1.D0 M=X N=M+1 X=X-FLOAT(M) FUNCTA=ACC(M)*(1.D0-X)+X*ACC(N) RETURN END C * * * * * * * * * * * * SUBROUTINE MULTPY(FINAL,AMATX,START) C *** TO EVALUATE PRODUCT OF B TIMES RR AND STORE RESULT IN TT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN DIMENSION FINAL(1850),AMATX(1850,1850),START(1850) C DO 10 IV=1,NSIZE FINAL(IV)=0.D0 DO 10 JV=1,NSIZE 10 FINAL(IV)=FINAL(IV)+AMATX(IV,JV)*START(JV) RETURN END C * * * * * SUBROUTINE DINV(AA,N0,N1,N2,DET,IND) IMPLICIT REAL*8(A-H,O-Z) DIMENSION AA(N0,N0),IPERM(2600),X(2600) N=N1 M=N+N2 NMX=N0+1 IF(0.GE.N.OR.N2.LT.0.OR.N.GE.NMX.OR.M.GE.NMX) GO TO 80 IND=1 DO 1 I=1,N 1 IPERM(I)=I EPS=0.D0 DO 2 K=1,N RMAX=0.D0 DO 3 J=K,N V=DABS(AA(K,J)) IF(RMAX-V) 4,3,3 4 RMAX=V L=J 3 CONTINUE IF(EPS-RMAX) 5,6,6 6 IF(EPS*0.01D0-RMAX) 7,8,8 8 DET=0.D0 IND=3 DO 9 I=1,N DO 9 J=1,N 9 AA(I,J)=1.0038D0 WRITE(6,200) GO TO 11 7 IND=2 5 PIVOT=AA(K,L) PIVI=1.D0/PIVOT IF(L-K) 12,13,12 12 IW=IPERM(K) IPERM(K)=IPERM(L) IPERM(L)=IW DO 14 I=1,N W=AA(I,K) AA(I,K)=AA(I,L) AA(I,L)=W 14 CONTINUE 13 CONTINUE DO 15 J=1,M X(J)=AA(K,J)*PIVI AA(K,J)=X(J) 15 CONTINUE DO 16 I=1,N IF(I-K) 17,16,17 17 W=AA(I,K) IF(W) 18,16,18 18 DO 19 J=1,M IF(J-K) 20,19,20 20 AA(I,J)=-W*X(J)+AA(I,J) 19 CONTINUE AA(I,K)=-W*PIVI 16 CONTINUE AA(K,K)=PIVI EPS=DMAX1(RMAX*1.D-33,EPS) 2 CONTINUE DO 21 I=1,N 22 K=IPERM(I) IF(K-I) 23,21,23 23 IW=IPERM(K) IPERM(K)=IPERM(I) IPERM(I)=IW DO 24 J=1,M W=AA(I,J) AA(I,J)=AA(K,J) AA(K,J)=W 24 CONTINUE GO TO 22 21 CONTINUE 11 RETURN 80 CONTINUE WRITE(6,201) N,N2 IND=4 GO TO 11 201 FORMAT(/'N1=',I5,' N2=',I5,' MEMORY-OVER') 200 FORMAT(/'THE GIVEN MATRIX IS SINGULAR') END C * * * * * * * * * * * * SUBROUTINE MATINV(AA,N0,N1,N2,DET,IND) IMPLICIT REAL*8(A-H,O-Z) DIMENSION AA(N0,N0),IPERM(1850),X(1850) N=N1 M=N+N2 NMX=N0+1 IF(0.GE.N.OR.N2.LT.0.OR.N.GE.NMX.OR.M.GE.NMX) GO TO 80 IND=1 DO 1 I=1,N 1 IPERM(I)=I EPS=0.D0 DO 2 K=1,N RMAX=0.D0 DO 3 J=K,N V=DABS(AA(K,J)) IF(RMAX-V) 4,3,3 4 RMAX=V L=J 3 CONTINUE IF(EPS-RMAX) 5,6,6 6 IF(EPS*0.01D0-RMAX) 7,8,8 8 DET=0.D0 IND=3 DO 9 I=1,N DO 9 J=1,N 9 AA(I,J)=1.0038D0 WRITE(6,200) GO TO 11 7 IND=2 5 PIVOT=AA(K,L) PIVI=1.D0/PIVOT IF(L-K) 12,13,12 12 IW=IPERM(K) IPERM(K)=IPERM(L) IPERM(L)=IW DO 14 I=1,N W=AA(I,K) AA(I,K)=AA(I,L) AA(I,L)=W 14 CONTINUE 13 CONTINUE DO 15 J=1,M X(J)=AA(K,J)*PIVI AA(K,J)=X(J) 15 CONTINUE DO 16 I=1,N IF(I-K) 17,16,17 17 W=AA(I,K) IF(W) 18,16,18 18 DO 19 J=1,M IF(J-K) 20,19,20 20 AA(I,J)=-W*X(J)+AA(I,J) 19 CONTINUE AA(I,K)=-W*PIVI 16 CONTINUE AA(K,K)=PIVI EPS=DMAX1(RMAX*1.E-33,EPS) 2 CONTINUE DO 21 I=1,N 22 K=IPERM(I) IF(K-I) 23,21,23 23 IW=IPERM(K) IPERM(K)=IPERM(I) IPERM(I)=IW DO 24 J=1,M W=AA(I,J) AA(I,J)=AA(K,J) AA(K,J)=W 24 CONTINUE GO TO 22 21 CONTINUE 11 RETURN 80 CONTINUE WRITE(6,201) N,N2 IND=4 GO TO 11 201 FORMAT(/'N1=',I5,' N2=',I5,' MEMORY-OVER') 200 FORMAT(/'THE GIVEN MATRIX IS SINGULAR') END C * * * * * * * * * * * * FUNCTION VFUNC(C) IMPLICIT REAL*8(A-H,O-Z) VFUNC=0.D0 IF(C.GT.0.D0) VFUNC=C RETURN END C * * * * * * * * * * * * SUBROUTINE PRINCE(SX,SY,TA,S1,S3,T) C *** CALCULATES PRINCIPAL STRESSES IMPLICIT REAL*8(A-H,O-Z) C=0.5D0*(SX+SY) A1=(SY-SX)*0.5D0 A2=A1*A1+TA*TA A=0.D0 IF(A2.GT.0.D0) A=DSQRT(A2) S1=C+A S3=C-A IF(SY.EQ.S3) GO TO 1 T1=TA/(SY-S3) T=DATAN(T1) GO TO 2 1 T=3.14159D0*0.5D0 2 RETURN END C * * * * * * * * * * * * SUBROUTINE MULT(A,B,C,L,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(L,M),B(M,N),C(L,N) DO 1 I=1,L DO 1 J=1,N C(I,J)=0.D0 DO 1 K=1,M 1 C(I,J)=C(I,J)+A(I,K)*B(K,J) RETURN END C * * * * * * * * * * * * SUBROUTINE WULT(A,B,C,L,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M,L),B(M,N),C(L,N) DO 1 I=1,L DO 1 J=1,N C(I,J)=0.D0 DO 1 K=1,M 1 C(I,J)=C(I,J)+A(K,I)*B(K,J) RETURN END C * * * * * * * * * * * * SUBROUTINE XULT(A,B,C,L,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(L,M),B(N,M),C(L,N) DO 1 I=1,L DO 1 J=1,N C(I,J)=0.D0 DO 1 K=1,M 1 C(I,J)=C(I,J)+A(I,K)*B(J,K) RETURN END C * * * * * * * * * * * * SUBROUTINE VULT(A,U,V,M,N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M,N),U(N),V(M) DO 1 I=1,M V(I)=0.D0 DO 1 J=1,N 1 V(I)=V(I)+A(I,J)*U(J) RETURN END C * * * * * * * * * * * * SUBROUTINE OUTDYN(IITER,ISTEP) C *** OUTPUT ROUTINE IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB5/DMX(1260,3,3),DMI(1260,3,3),BMX(1260,3,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LAB8/IFIXD,MITER,IPRED,NCHEK COMMON/LAB9/DTIME,DTEND,DTREC,DELTA,GAAMA COMMON/LA11/IREQD(10),IREQS(10),IWA(1260) COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA17/VELOI(1850),VELOL(1850),VELOT(1850) COMMON/LA18/ACCEI(1850),ACCEL(1850) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) COMMON/LA22/NTENS,KTENS,ITENS(1260),JTENS(1260) COMMON/LA23/TMX(1260,3,3),TMI(1260,3,3),EFFST(1260) DIMENSION STRPR(1260,7),UX(1280),UY(1280),UXH(20), * VX(1280),AX(1280),VXH(20),AXH(20),TT(1260) DIMENSION IYITE(2150) C KSTEP=ISTEP IOUT=NOUTD(KSTEP) C XTIME=FLOAT(KSTEP)*DTIME WRITE(6,200) KSTEP,XTIME C *** REARRANGE DISPLACEMENT VECTOR DO 10 IP=1,NPOIN I1=(IP-1)*3+1 IT=LY(I1) UX(IP)=DISPI(IT) VX(IP)=VELOI(IT) AX(IP)=ACCEI(IT) I2=(IP-1)*3+2 IT=LY(I2) UY(IP)=DISPI(IT) 10 CONTINUE C IF(IOUT.LE.0) GO TO 11 WRITE(6,201) DO 12 IP=1,NPOIN 12 WRITE(6,202) IP,UX(IP),UY(IP) 11 CONTINUE DO 52 ID=1,NREQD IP=IREQD(ID) WRITE(6,202) IP,UX(IP),UY(IP) UXH(ID)=UX(IP) VXH(ID)=VX(IP) AXH(ID)=AX(IP) 52 CONTINUE WRITE(8,205) XTIME,(UXH(ID),ID=1,NREQD) WRITE(9,205) XTIME,(VXH(ID),ID=1,NREQD) WRITE(10,205) XTIME,(AXH(ID),ID=1,NREQD) C *** WRITES STRESSES DO 14 IE=1,NELEM KOL2=K2(IE) GO TO (21,21,21,21,24,21),KOL2 21 CONTINUE DO 28 I=1,6 28 STRPR(IE,I)=STRSG(IE,I) STRPR(IE,7)=0.D0 A=1.D0 IF(KOL2.EQ.1.OR.KOL2.EQ.3) A=PRP(MT,2) TT(IE)=STRSG(IE,1)*A GO TO 20 24 CONTINUE SX=STRSG(IE,1) SY=STRSG(IE,2) TA=STRSG(IE,3) CALL PRINCE(SX,SY,TA,S1,S3,AG) STRPR(IE,1)=SX STRPR(IE,2)=SY STRPR(IE,3)=TA STRPR(IE,4)=S1 STRPR(IE,5)=S3 STRPR(IE,6)=PI3(IE) STRPR(IE,7)=PI4(IE) 20 CONTINUE IF(IOUT.LE.0) GO TO 14 WRITE(6,204) IE,IYIEL(IE),(STRPR(IE,I),I=1,7) 14 CONTINUE DO 30 ID=1,NREQS IE=IREQS(ID) WRITE(6,204) IE,IYIEL(IE),(STRPR(IE,I),I=1,7) 30 CONTINUE DO 32 IE=1,NELEM IYITE(IE)=0 IF(JTENS(IE).GE.1) IYITE(IE)=2 IF(IYIEL(IE).GE.1) IYITE(IE)=1 32 CONTINUE WRITE(6,207) (IYITE(IE),IE=1,NELEM) C *** OUTPUT FOR INITIAL STRESSES OPEN(7,FILE='DASTRE0',STATUS='UNKNOWN') DO 34 I=1,3 34 WRITE(7,206) (STRPR(IE,I),IE=1,NELEM) CLOSE(7) C 200 FORMAT(/'TIME STEP',I5,3X,'TIME ',E12.4) 201 FORMAT(/'NODE',6X,'X-DISP',6X,'Y-DISP') 202 FORMAT(I5,2E12.3) 203 FORMAT(/'STRESSES'/' EL YIELD UNLOAD YIELD-HISTORY', * 1X,'X-ST',5X,'Y-ST', * 4X,'XY-ST',5X,'Z-ST',6X,'MAX',6X,'MIN',1X,'S1-ANG') 204 FORMAT(I5,I3,2X,9E10.3) 205 FORMAT(5F12.8) 206 FORMAT(F6.3,8F8.2) 207 FORMAT(24I2) RETURN END C * * * * * * * * * * * * SUBROUTINE MAVS1 C *** MICRO-AVS10 OUTPUT (PART 1: DIRECTION OF SHEAR BAND) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*5 AC(900) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB6/XX(1280),YY(1280) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) COMMON/LA22/NTENS,KTENS,ITENS(1260),JTENS(1260) DIMENSION SX(1260),SY(1260),XXX(1280),YYY(1280) DIMENSION J(1260) C DO 2 L=1,NELEM J(L)=1 SX(L)=0.D0 SY(L)=0.D0 IF(IYIEL(L).LE.0.AND.JTENS(L).LE.0) GO TO 2 IF(JTENS(L).GE.1) SI=3.14159D0*0.5D0 IF(IYIEL(L).GE.1) SI=PI4(L) SX(L)=DCOS(SI) SY(L)=DSIN(SI) 2 CONTINUE C NELEM1=0 DO 6 L=1,NELEM IF(IJK(L,3).EQ.0) GO TO 6 NELEM1=NELEM1+1 6 CONTINUE WRITE(11,301) NPOIN,NELEM1 C DO 8 I=1,NPOIN XXX(I)=XX(I)*100.D0 YYY(I)=YY(I)*100.D0 WRITE(11,302) I,XXX(I),YYY(I),0 8 CONTINUE C DO 10 L=1,NELEM AC(L)='quad' IF(IJK(L,4).EQ.0) AC(L)='tri' IF(IJK(L,3).EQ.0) GO TO 10 WRITE(11,303) L,J(L),AC(L),(IJK(L,I),I=1,4) 10 CONTINUE WRITE(11,304) WRITE(11,305) WRITE(11,306) DO 12 L=1,NELEM 12 WRITE(11,307) L,SX(L),SY(L) C 301 FORMAT(2I5,4X,'0',4X,'2',4X,'0') 302 FORMAT(I5,2F12.3,F9.3) 303 FORMAT(2I5,A8,I4,20I5) 304 FORMAT(3X,'2',3X,'1',4X,'1') 305 FORMAT('strX , (kN/m2)') 306 FORMAT('strY , (kN/m2)') 307 FORMAT(I5,3E13.4) RETURN END C * * * * * * * * * * * * SUBROUTINE MAVS2 C *** MICRO-AVS10 OUTPUT (PART 2: DISPLACEMENTS & YIELD ELEMENTS) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*5 AC(900) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB2/LY(3500),LLL(1260,8),IX(1280),IY(1280),IQ(1280) COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB6/XX(1280),YY(1280) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) COMMON/LA22/NTENS,KTENS,ITENS(1260),JTENS(1260) DIMENSION UX(1280),UY(1280),UM(1280),WST(1260,6),PI9(1260) C DO 2 IP=1,NPOIN I1=(IP-1)*3+1 IT=LY(I1) UX(IP)=DISPI(IT) I2=(IP-1)*3+2 IT=LY(I2) UY(IP)=DISPI(IT) 2 CONTINUE WRITE(12,310) 1 WRITE(12,311) WRITE(12,300) 1 C WRITE(12,301) NPOIN,NELEM DO 10 I=1,NPOIN 10 WRITE(12,302) I,XX(I),YY(I),0 DO 12 L=1,NELEM AC(L)='quad' IF(IJK(L,4).EQ.0) AC(L)='tri' IF(IJK(L,3).EQ.0) AC(L)='line' WRITE(12,303) L,K2(L),AC(L),(IJK(L,I),I=1,4) 12 CONTINUE WRITE(12,304) WRITE(12,305) WRITE(12,306) WRITE(12,307) WRITE(12,308) DO 14 I=1,NPOIN 14 WRITE(12,309) I,UX(I),UY(I),0 WRITE(12,313) WRITE(12,314) DO 20 L=1,NELEM PI9(L)=0.1D0 IF(JTENS(L).GE.1) PI9(L)=1.D0 IF(IYIEL(L).GE.1) PI9(L)=2.D0 20 CONTINUE DO 24 L=1,NELEM 24 WRITE(12,315) L,PI9(L) C 300 FORMAT('step',I1) 301 FORMAT(2I5) 302 FORMAT(I5,3F10.4) 303 FORMAT(2I5,A6,20I5) 304 FORMAT(4X,'3',4X,'1') 305 FORMAT(4X,'3',4X,'1',4X,'1',4X,'1') 306 FORMAT(1X,'disp_x, m') 307 FORMAT(1X,'disp_y, m') 308 FORMAT(1X,'disp_z, m') 309 FORMAT(I5,3E14.4) 310 FORMAT(I5) 311 FORMAT('data_geom') 313 FORMAT(4X,'1',4X,'1') 314 FORMAT(1X,'no') 315 FORMAT(I5,E12.3) RETURN END C * * * * * * * * * * * * SUBROUTINE MAVS3 C *** MICRO-AVS10 OUTPUT (PART 3: PRINCIPAL STRESSES) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*5 AC(900) COMMON/LAB1/NPOIN,NELEM,NMATS,NIST,NPREV,ICRI,NSIZE,NN COMMON/LAB3/IJK(1260,4),K2(1260),MAT(1260),IST(1260),LRE(1260) COMMON/LAB4/PRP(20,12),AES(1260),EKK(1260,8,8),DBM(1260,6,8) COMMON/LAB6/XX(1280),YY(1280) COMMON/LAB7/NSTEP,NOUTP,NOUTD(6600),NREQD,NREQS,NACCE,IFUNC COMMON/LA13/FORCE(1850),DEP(1260,3,3),STAX(1260,3) COMMON/LA15/STRAG(1260,6),STRSG(1260,6),PW(1260,3) COMMON/LA16/DISPI(1850),DISPL(1850),DISPT(1850) COMMON/LA20/IYIEL(1260),IYIEP(1260),IUNLO(1260),IYIEH(1260) COMMON/LA21/PI1(1260),PI2(1260),PI3(1260),PI4(1260) COMMON/LA22/NTENS,KTENS,ITENS(1260),JTENS(1260) DIMENSION WST(6),AX(1260),AY(1260),XXX(1280),YYY(1280) DIMENSION J(1690) C DO 2 L=1,NELEM J(L)=0 AX(L)=0.D0 AY(L)=0.D0 2 CONTINUE C DO 10 L=1,NELEM KOL2=K2(L) IF(KOL2.LE.3) GO TO 10 IF(IYIEL(L).GE.1) GO TO 12 DO 20 I=1,3 20 WST(I)=STRSG(L,I) GO TO 14 12 CONTINUE DO 22 I=1,3 22 WST(I)=STAX(L,I) 14 CONTINUE SX=WST(1) SY=WST(2) TA=WST(3) CALL PRINCE(SX,SY,TA,S1,S3,TH) AX(L)=S1*DSIN(TH) AY(L)=S1*DCOS(TH) 10 CONTINUE C NNE1=0 DO 6 L=1,NELEM IF(IJK(L,3).EQ.0) GO TO 6 NNE1=NNE1+1 6 CONTINUE WRITE(13,301) NPOIN,NNE1 C DO 30 I=1,NPOIN XXX(I)=XX(I)*100.D0 YYY(I)=YY(I)*100.D0 WRITE(13,302) I,XXX(I),YYY(I),0 30 CONTINUE DO 32 L=1,NELEM AC(L)='quad' IF(IJK(L,4).EQ.0) AC(L)='tri' IF(IJK(L,3).EQ.0) GO TO 32 WRITE(13,303) L,J(L),AC(L),(IJK(L,I),I=1,4) 32 CONTINUE WRITE(13,304) WRITE(13,305) WRITE(13,306) DO 34 L=1,NELEM 34 WRITE(13,307) L,AX(L),AY(L) C 301 FORMAT(2I5,4X,'0',4X,'2',4X,'0') 302 FORMAT(I5,2F12.3,F9.3) 303 FORMAT(2I5,A8,I4,20I5) 304 FORMAT(3X,'2',3X,'1',4X,'1') 305 FORMAT('strX , (kN/m2)') 306 FORMAT('strY , (kN/m2)') 307 FORMAT(I5,3E13.4) RETURN END