PROGRAM LSFEA1 C LIMIT STATE FINITE ELEMENT ANALYSIS (TYPE-1) C SLOPE STABILITY, AND BEARING CAPACITY PROBLEMS C MODIFIED INITIAL STRESS METHOD C SHEAR-BAND = COULOMB FAILURE SURFACE C MONOTONOUSLY LOADING C PORE-WATER PRESSURE C IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB3/STR(1690,6),EPS(1690,3) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB16/DEP(1690,3,3),GG1(1690,3,2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB19/IWE(20),NRE,IRE(50),LRE(1690) COMMON/LB20/NTS,ITS(1500),JTS(1690) COMMON/LB21/PS1(1690),PS3(1690),QS1(1690),QS3(1690),LFA(1690) COMMON/LB22/FIR,RKH,SCD,ERC,FLO(60),DWE(20) C CALL INPUTS CALL PREPAR(IND) IF(IND.NE.1) GO TO 10 CALL INISTR C 10 CONTINUE CLOSE(6) CLOSE(8) CLOSE(9) CLOSE(10) CLOSE(11) CLOSE(12) CLOSE(13) STOP END C * * * * * SUBROUTINE INPUTS C *** DATA INPUT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3),IST(1690) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB19/IWE(20),NRE,IRE(50),LRE(1690) COMMON/LB20/NTS,ITS(1500),JTS(1690) COMMON/LB22/FIR,RKH,SCD,ERC,FLO(60),DWE(20) DIMENSION PNA(20),STII(50,3),TY(1690),EE1(30),PP1(30),HH1(30), * RO1(30),CC1(30),FA1(30),DL1(30),FAI(1690),DLT(1690) DIMENSION IDATA(12),IDATAQ(12),IBA(1690),NET(6),NM(1690), * MAT(1690),JSTP(30) C OPEN(5,FILE='DALSFE',STATUS='UNKNOWN') OPEN(6,FILE='PRLSFE',STATUS='UNKNOWN') OPEN(8,FILE='PRDISP',STATUS='UNKNOWN') OPEN(9,FILE='M-AVS.inp',STATUS='UNKNOWN') OPEN(10,FILE='L-AVS',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 C NN: MAXIMUM NUMBER OF DISPLACEMENT VARIABLES C NVM: MAXIMUM NUMBER OF INITIAL-STRESS VARIABLES NN=2850 NVM=2150 DO 2 I=1,12 2 IDATAQ(I)=0 DO 4 I=1,6 4 NET(I)=0 C PNA(I): PROJECT NAME READ(5,100) (PNA(I),I=1,15) WRITE(6,100) (PNA(I),I=1,15) C NNP: NUMBER OF NODES, NNE:NUMBER OF ELEMENTS C NMAT: NUMBER OF MATERIALS C NFL: NUMBER OF LOAD POINTS C NWE: NUMBER OF RETAINING-WALL ELEMENTS C NIS: NUMBER OF SETS OF ACTUAL INITIAL-STRESSES C ITER=0: NO ITERATION, 1:ITERATION C NTRY: NUMBER OF ITERATIONS FOR NO-TENSION ELEMENT C ICR=0: NOT CONSIDER CRITICAL STATE C 1: CONSIDER AFTER FIRST YIELDING (PERFECTLY PLASTIC) C ONLY FOR MOHR-COULOMB MATERIAL C 2: CONSIDER AFTER FIRST YIELDING (PERFECTLY PLASTIC) C BOTH FOR MOHR-COULOMB & COULOMB MATERIALS C ISTA=0: TAKE INPUT DATA AS UNYIELD STRESSES, WHEN FINDING YIELD STRESSES C 1: TAKE STRESSES DECIDED MEAN OR NORMAL STRESS 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 NSTP: NUMBER OF LOADING-STEPS C MSTP: ACTUAL CALCULATION NUMBER OF LOADING-STEPS C KSTP: NUMBER OF LOADING-STEPS FOR MICRO-AVS OUTPUT C JSTP: STEP NO. OUTPUT FOR MICRO-AVS C READ(5,101) NNP,NNE,NMAT,NFL,NWE,NIS,ITER,NTRY READ(5,101) ICR,ISTA,IPR,IPW READ(5,101) NSTP,MSTP,KSTP READ(5,101) (JSTP(N),N=1,KSTP) IF(NSTP.LE.0) NSTP=1 IF(MSTP.LE.0) MSTP=NSTP IF(NTRY.LE.0) NTRY=1 IDY=0 WRITE(6,200) NNP,NNE,NMAT,NFL,NWE,NIS,ITER,NTRY, * ICR,ISTA,IPR,IPW, * NSTP,MSTP,KSTP,(JSTP(N),N=1,KSTP) DO 6 N=1,NSTP 6 LSTP(N)=0 DO 8 I=1,KSTP N=JSTP(I) LSTP(N)=1 8 CONTINUE C 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 I=1,NNP READ(5,102) K,XX(I),YY(I),IX(I),IY(I),IQ(I) IF(IPR.EQ.0) GO TO 10 WRITE(6,102) K,XX(I),YY(I),IX(I),IY(I),IQ(I) 10 CONTINUE C *** ELEMENT DATA C K: ELEMENT NO., IJK(L,I): NODE NO. C K2(L)= 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(L): MATERAL NO. C IBA(L)=1: CONSIDER OWN WEIGHT, 0:NOT CONSIDER C IST(L): NO. OF SET OF ACTUAL INITIAL STRESSES C LRE(L): TYPE OF SHEAR BAND DIRECTION (1 OR 2) C *** IF DATA(I)=0, PRECEDING VALUE IS EMPLOYED C DO 12 L=1,NNE READ(5,103) K,(IJK(L,I),I=1,4),(IDATA(I),I=1,5) DO 14 I=1,12 IF(IDATA(I).EQ.0) IDATA(I)=IDATAQ(I) IF(IDATA(I).LT.0) IDATA(I)=0 14 IDATAQ(I)=IDATA(I) K2(L)= IDATA(1) MAT(L)=IDATA(2) IBA(L)=IDATA(3) IST(L)=IDATA(4) LRE(L)=IDATA(5) 12 CONTINUE C *** MATERIAL PARAMETERS C EE1(M): YOUNG'S MODULUS FOR K2=1-5, SHEAR RIGIDITY FOR K2=6 C PP1(M): POISSON'S RATIO FOR K2=4 & 5, C CROSS AREA FOR K2=1-3, YOUNG'S MODULUS FOR K2=6 C HH1(M): 1.0 FOR K2=4 & 5, MOMENT OF INERTIA FOR K2=3 C POISSON'S RATIO FOR K2=6 C RO1(M): DENSITY C CC1(M): COHESION, FA1(M):ANGLE OF SHEAR RESISTANCE C DL1(M): ANGLE OF DILATANCY C WRITE(6,213) DO 16 M=1,NMAT READ(5,107) M1,EE1(M),PP1(M),HH1(M),RO1(M),CC1(M), * FA1(M),DL1(M) WRITE(6,107) M1,EE1(M),PP1(M),HH1(M),RO1(M),CC1(M), * FA1(M),DL1(M) FA1(M)=FA1(M)*3.14159D0/180.D0 DL1(M)=DL1(M)*3.14159D0/180.D0 16 CONTINUE C DO 18 L=1,NNE M=MAT(L) EE(L)=EE1(M) PP(L)=PP1(M) HH(L)=HH1(M) ROU(L)=RO1(M) CC(L)=CC1(M) FI=FA1(M) FAI(L)=FI*180.D0/3.14159D0 SFI(L)=DSIN(FI) CFI(L)=DCOS(FI) DL=DL1(M) DLT(L)=DL*180.D0/3.14159D0 SDL(L)=DSIN(DL) CDL(L)=DCOS(DL) N=K2(L) NET(N)=NET(N)+1 IF(N.LE.3.OR.N.EQ.6) ROU(L)=0.D0 18 CONTINUE C *** ELEMENT TYPE WRITE(6,202) (NET(I),I=1,6) C *** EMBANKMENT ELEMENTS NBA=0 DO 20 L=1,NNE IB=IBA(L) IF(IB.EQ.0) GO TO 20 NBA=NBA+1 LBA(NBA)=L 20 CONTINUE IF(NBA.LE.0) GO TO 22 WRITE(6,203) NBA WRITE(6,101) (LBA(I),I=1,NBA) 22 CONTINUE C WRITE(6,205) DO 26 L=1,NNE IF(IPR.EQ.0) GO TO 26 WRITE(6,206) L,(IJK(L,I),I=1,4),K2(L),EE(L),PP(L),HH(L), * ROU(L),CC(L),FAI(L),DLT(L),IBA(L),IST(L),LRE(L) 26 CONTINUE C *** NODAL LOADS C NOD: NODE NO., IXY=1: X-DIRECTION, 2:Y-DIRECTION, 3:MOMENT, FLO(I):LOAD IF(NFL.EQ.0) GO TO 30 WRITE(6,207) DO 32 I=1,NFL READ(5,104) NOD,IXY,FLO(I) WRITE(6,104) NOD,IXY,FLO(I) 32 LOI(I)=3*NOD-(3-IXY) 30 CONTINUE C *** WALL ELEMENTS C IWE(I): WALL ELEMENT NO., DWE(I): LENGTH OF WALL ELEMENT (K2=6) IF(NWE.EQ.0) GO TO 34 READ(5,101) (IWE(I),I=1,NWE) WRITE(6,209) (IWE(I),I=1,NWE) READ(5,108) (DWE(I),I=1,NWE) WRITE(6,106) (DWE(I),I=1,NWE) 34 CONTINUE C *** ACTUALLY INITIAL STRESSES C K1: SET NO., STII(K,I) I=1,2,3: SIGUMA-X,Y,TAU-XY DO 40 L=1,NNE DO 40 I=1,3 40 PST(L,I)=0.D0 IF(NIS.LE.0) GO TO 50 WRITE(6,210) DO 44 K=1,NIS READ(5,107) K1,(STII(K,I),I=1,3) 44 WRITE(6,107) K1,(STII(K,I),I=1,3) DO 46 L=1,NNE IF(K2(L).LE.3) GO TO 46 K=IST(L) DO 48 I=1,3 48 PST(L,I)=STII(K,I) 46 CONTINUE READ(5,101) IU WRITE(6,212) IU IF(IU.LE.0) GO TO 50 OPEN(7,FILE='DASTRE',STATUS='UNKNOWN') DO 52 I=1,3 READ(7,108) (PST(L,I),L=1,NNE) 52 WRITE(6,108) (PST(L,I),L=1,NNE) 50 CONTINUE CLOSE(7) C *** SEISMIC COEFFICIENT RKH READ(5,108) RKH WRITE(6,204) RKH C *** CONVERGENCE CONSTANTS ETC. C SCD: SCALING FACTOR FOR DISPLACEMENT OUTPUT C ERC: CONV. CONST. FOR CONFINING PRESSURE C FSS: SAFETY FACTOR FS C READ(5,106) SCD,ERC,FSS IF(SCD.LT.0.01) SCD=1. IF(THB.LT.0.001) THB=0.001 IF(FSS.LE.0.) FSS=1.D0 WRITE(6,208) SCD,ERC,FSS C C *** DEVIDE BY SAFETY FACTOR FSS DO 56 L=1,NNE CC(L)=CC(L)/FSS FI1=FAI(L)*3.14159D0/180.D0 TF1=DSIN(FI1)/DCOS(FI1) TFI=TF1/FSS FI=DATAN(TFI) SFI(L)=DSIN(FI) CFI(L)=DCOS(FI) FAI(L)=FI*180.D0/3.14159D0 DL=DLT(L)*3.14159D0/180.D0 SDL(L)=DSIN(DL) CDL(L)=DCOS(DL) 56 CONTINUE WRITE(6,214) WRITE(6,108) (CC(L),L=1,NNE) WRITE(6,108) (FAI(L),L=1,NNE) C NDI: DISPLACEMENT OUTPUT NODE C FIR: RATIO TO TERZAGHI BEARING CAPACITY READ(5,102) NDI,FIR WRITE(6,227) NDI,FIR IS=0 R=0.D0 DI=0.D0 WRITE(8,102) IS,R,DI C *** INITIAL DISPLACEMENTS READ(5,101) IU WRITE(6,211) IU DO 60 I=1,NN 60 PTU(I)=0.D0 IF(IU.LE.0) GO TO 64 OPEN(7,FILE='DADISP',STATUS='UNKNOWN') READ(7,108) (PTU(I),I=1,NN) DO 62 I=1,NN 62 PTU(I)=-PTU(I) CLOSE(7) 64 CONTINUE C *** EQUIVALENT NODAL FORCE C NRE: NUMBER OF NODES FOR CALCULATING EQUIVALENT NODAL FORCE C IRE(I): NODE NO. READ(5,101) NRE IF(NRE.LE.0) GO TO 66 READ(5,101) (IRE(I),I=1,NRE) 66 CONTINUE C C *** PORE WATER PRESSURE WRITE(6,215) DO 68 L=1,NNE DO 68 I=1,2 68 PW(L,I)=0.D0 IF(IPW.LE.1) GO TO 78 DO 80 I=1,2 80 READ(5,108) (STB(L,I),L=1,NNE) DO 82 I=1,2 82 READ(5,108) (STC(L,I),L=1,NNE) DO 84 L=1,NNE DO 84 I=1,2 84 PW(L,I)=STB(L,I)-STC(L,I) DO 86 I=1,2 86 WRITE(6,108) (PW(L,I),L=1,NNE) 78 CONTINUE CLOSE(5) C C *** MICRO-AVS7 OUTPUT WRITE(9,216) WRITE(9,217) WRITE(9,218) WRITE(9,219) WRITE(9,101) NNP,NNE DO 90 I=1,NNP WRITE(9,220) I,XX(I),YY(I),ZZ(I) 90 CONTINUE DO 92 L=1,NNE TY(L)=' quad' IF(IJK(L,3).EQ.0.AND.IJK(L,4).EQ.0) TY(L)=' line' IF(IJK(L,3).NE.0.AND.IJK(L,4).EQ.0) TY(L)=' tri' WRITE(9,221) L,K2(L),TY(L),(IJK(L,I),I=1,4) 92 CONTINUE WRITE(9,222) WRITE(9,223) WRITE(9,224) WRITE(9,225) WRITE(9,226) C 100 FORMAT(15A4) 101 FORMAT(15I5) 102 FORMAT(I5,2F10.3,3I5) 103 FORMAT(10I5) 104 FORMAT(2I5,5F10.2/10X,5F10.2) 105 FORMAT(3I5,4E10.3,2I5) 106 FORMAT(8E10.3) 107 FORMAT(I5,7E10.3) 108 FORMAT(10E8.3) 200 FORMAT(/'NODES=',I5,' ELEMENTS=',I5,' MATERIALS=',I3 */' LOADS=',I3,' WALL-ELEMENTS=',I3,' INITIAL-STRESSES=',I3 */' ITERATION=',I2,' TRYALS=',I3,' CRITICAL=',I2 */' YIELD STRESSES=',I2,' PRINT=',I2,' PORE-WATER PRESSURE=',I2 */' LOAD-STEPS=',I5,' ACTUAL STEPS=',I5,' M-AVS STEPS='I5 */' M-AVS STEP NO.'/5X,10I5) 201 FORMAT(/'NODAL DATA'/'NO.',9X,'X',9X,'Y',' X-F Y-F R-F') 202 FORMAT(/'TRUSS=',I5,' BEAM=',I5,' TEXTILE=',I5 */' PLANE-STRAIN (ELASTIC)=',I5,' PLANE-STRAIN (NONLINEAR)=',I5 */' INTERFACE=',I5) 203 FORMAT(/'EMBANKMENT ELEMENTS=',I5) 204 FORMAT(/'SEISMIC COEFFCIENT'/(10F7.3)) 205 FORMAT(/'ELEMENT DATA' */'NO. NODES TYPE E P H ROH C FAI DEL EMB STRE SLIP') 206 FORMAT(I4,4I5,I2,3E9.2,F6.2,F7.2,2F6.2,3I2) 207 FORMAT(/'NODAL LOADS'/' NODE X-Y LOAD') 208 FORMAT(/'DISPLACEMENT SCALE=',F10.2 */' CRITICAL CHECK=',F10.4,' SAFETY-FACTOR FS=',F10.3) 209 FORMAT(/'WALL ELEMENTS'/(10I5)) 210 FORMAT(/'ACTUALLY INITIAL STRESSES'/' SET SX SY TAU') 211 FORMAT(/'INITIAL DISPLACEMENT 0=NO, 1=INPUT',I5) 212 FORMAT(/' INITIAL STRESSES 0=NO, 1=INPUT',I5) 213 FORMAT(/'MATERIAL PARAMETERS' */' NO. E NYU H ROH C FAI DEL') 214 FORMAT(/'STRENGTH PARAMETERS C & FAI DEVIDED BY FS') 215 FORMAT(/'PORE-WATER PRESSURE') 216 FORMAT('# UCD data') 217 FORMAT('1') 218 FORMAT('data_geom') 219 FORMAT('step1') 220 FORMAT(I5,3E12.4) 221 FORMAT(2I5,A5,4I5) 222 FORMAT('3 8') 223 FORMAT('3 1 1 1') 224 FORMAT('DISP_X,') 225 FORMAT('DISP_Y,') 226 FORMAT('DISP_Z,') 227 FORMAT(/'DISPLACEMENT OUTPUT NODE=',I5 */' FINAL PRESSURE RATIO=',F10.3) RETURN END C * * * * * SUBROUTINE PREPAR(IND) C *** SET DISPLACEMENT VECTOR AND STIFFNESS MATRIX IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB22/FIR,RKH,SCD,ERC,FLO(60),DWE(20) DIMENSION EML(1580) C DO 10 L=1,NNP 10 NDF(L)=2 DO 12 L=1,NNE IF(K2(L).NE.2) GO TO 12 I=IJK(L,1) J=IJK(L,2) NDF(I)=3 NDF(J)=3 12 CONTINUE DO 14 L=1,NNP LY(3*L-2)=IX(L) LY(3*L-1)=IY(L) IF(NDF(L).EQ.2) IQ(L)=1 14 LY(3*L)=IQ(L) NX=3*NNP LOC=0 DO 16 L=1,NX IF(LY(L).EQ.0) GO TO 18 LY(L)=NN GO TO 16 18 LOC=LOC+1 LX(LOC)=L LY(L)=LOC 16 CONTINUE NY=LOC WRITE(6,200) NY IF(NY.LT.NN) GO TO 20 WRITE(6,201) IND=10 RETURN C 20 DO 22 I=1,NY DO 22 J=1,NY 22 GKK(I,J)=0.D0 DO 30 L=1,NNE KOL2=K2(L) GO TO (31,32,31,34,34,34),KOL2 31 CALL TRUSS(L) GO TO 30 32 CALL BEAM(L) GO TO 30 34 CALL PLANES(L,KOL2) 30 CONTINUE C CALL DINV(GKK,NN,NY,0,DET,IND) IF(IND.EQ.1) GO TO 50 WRITE(6,202) IND C *** NODAL LOAD 50 CONTINUE DO 52 I=1,NY 52 TF(I)=0.D0 TF(NN)=0.D0 IF(NFL.EQ.0) GO TO 60 DO 54 I=1,NFL IT=LOI(I) IT=LY(IT) 54 TF(IT)=TF(IT)+FLO(I)/DFLOAT(NSTP) C *** EMBANKMENT LOAD 60 CONTINUE IF(NBA.LE.0) GO TO 80 DO 62 I=1,NNP 62 EML(I)=0.D0 DO 64 I=1,NBA LB=LBA(I) IF(K2(LB).LE.3.OR.K2(LB).EQ.6) GO TO 64 MM=IJK(LB,4) N1=4 C1=0.25 IF(MM.EQ.0) N1=3 IF(MM.EQ.0) C1=0.333 EMI=-AES(LB)*ROU(LB)*C1 DO 66 J=1,N1 II=IJK(LB,J) 66 EML(II)=EML(II)+EMI 64 CONTINUE DO 68 L=1,NNP IT=LY(3*L-2) JT=LY(3*L-1) TF(IT)=TF(IT)+EML(L)*RKH/DFLOAT(NSTP) TF(JT)=TF(JT)+EML(L)/DFLOAT(NSTP) 68 CONTINUE WRITE(6,203) (EML(I),I=1,NNP) C *** ORIGINAL LOAD VECTOR 80 CONTINUE DO 70 I=1,NY FF(I)=0.D0 DO 70 J=1,NY 70 FF(I)=FF(I)+GKK(I,J)*TF(J) FF(NN)=0.D0 C 200 FORMAT(/'NUMBER OF VARIABLES=',I4) 201 FORMAT(/'* VARIABLE MEMORY=OVER') 202 FORMAT(/'* INDEX=',I3) 203 FORMAT(/'EMBANK. LOAD'/(10F7.2)) RETURN END C * * * * * SUBROUTINE INISTR C *** APPLY MODIFIED INITIAL STRESS METHOD IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB3/STR(1690,6),EPS(1690,3) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB16/DEP(1690,3,3),GG1(1690,3,2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB19/IWE(20),NRE,IRE(50),LRE(1690) COMMON/LB20/NTS,ITS(1500),JTS(1690) COMMON/LB21/PS1(1690),PS3(1690),QS1(1690),QS3(1690),LFA(1690) COMMON/LB22/FIR,RKH,SCD,ERC,FLO(60),DWE(20) C NFAT=0 NTS=0 DO 2 L=1,NNE JFA(L)=0 MFA(L)=0 LFA(L)=0 JTS(L)=0 PS3(L)=0.D0 QS3(L)=0.D0 DO 4 I=1,3 SST(L,I)=0.D0 4 STAX(L,I)=0.D0 2 CONTINUE C DO 100 ITP=1,MSTP ISTP=ITP WRITE(6,200) ISTP DO 6 L=1,NNE DO 6 I=1,3 6 ST0(L,I)=0.D0 C DO 110 ITT=1,NTRY WRITE(6,201) ITT C CALL DISPLA CALL STRESS DO 10 I=1,NY 10 QTU(I)=PTU(I)+TU(I) DO 12 L=1,NNE DO 14 I=1,3 14 QST(L,I)=PST(L,I)+STR(L,I) IF(K2(L).LE.3) GO TO 12 DO 16 I=1,2 16 QST(L,I)=QST(L,I)-PW(L,I) 12 CONTINUE IF(ITT.EQ.NTRY) GO TO 302 CALL YIELD(QST) IF(ITER.EQ.0) GO TO 302 IF(ITT.EQ.1.AND.NFA.LE.0) GO TO 300 IF(ITT.GE.2.AND.NFA.LE.0) GO TO 302 WRITE(6,202) (IFA(I),I=1,NFA) CALL CALSTA DO 20 IE=1,NFA LI=IFA(IE) DO 22 I=1,3 22 STAX(LI,I)=STA(IE,I) SX=STAX(LI,1) SY=STAX(LI,2) TA=STAX(LI,3) CALL PRINCE(SX,SY,TA,S1,S3,TH) PI3(LI)=TH 20 CONTINUE CALL DEPMAT CALL YGRAD1 300 CONTINUE IF(NFAT.LE.0) GO TO 302 CALL CALSTC CALL NOTENS CALL YGRAD2 CALL YGRAD3 110 CONTINUE C 302 CONTINUE CALL NOTENS CALL YGRAD2 CALL YGRAD3 CALL DISPLA CALL STRESS DO 50 I=1,NY 50 QTU(I)=PTU(I)+TU(I) DO 52 L=1,NNE DO 54 I=1,3 54 QST(L,I)=PST(L,I)+STR(L,I) 52 CONTINUE DO 24 L=1,NNE DO 24 I=1,3 24 STB(L,I)=QST(L,I) IF(NFAT.LE.0) GO TO 304 CALL CALSTB(1) 304 CONTINUE CALL OUTPUT(1) OPEN(7,FILE='PSLSFE',STATUS='UNKNOWN') CALL OUTPUT(2) CLOSE(7) C DO 26 L=1,NNE 26 MFA(L)=JFA(L) DO 28 L=1,NNE DO 28 I=1,3 28 SST(L,I)=PST(L,I) DO 32 L=1,NNE DO 32 I=1,3 32 PST(L,I)=STB(L,I) DO 34 I=1,NY 34 PTU(I)=QTU(I) C IM=LSTP(ISTP) IF(IM.LE.0) GO TO 80 CALL MAVS1 CALL MAVS2 CALL MAVS3 80 CONTINUE 100 CONTINUE C 200 FORMAT(/'* LOADING STEP=',I5) 201 FORMAT('STAGE',I3) 202 FORMAT('NEW YIELD EL',10(/10I5)) RETURN END C * * * * * SUBROUTINE DISPLA C *** CALCULATE DISPLACEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB20/NTS,ITS(1500),JTS(1690) DIMENSION TR(2850),FE(8),BT(8,3) C DO 2 I=1,NY 2 TR(I)=0.D0 IF(NFAT.LE.0) GO TO 4 DO 10 IL=1,NFAT LI=KFA(IL) KOL2=K2(LI) NV=8 IF(IJK(LI,4).EQ.0) NV=6 DO 12 I=1,NV DO 12 J=1,3 BT(I,J)=0.D0 DO 12 K=1,3 12 BT(I,J)=BT(I,J)+BMX(LI,K,I)*TMX(LI,K,J) C DO 16 I=1,NV FE(I)=0.D0 DO 16 J=1,3 16 FE(I)=FE(I)+BT(I,J)*ST0(LI,J) A=AES(LI) DO 18 I=1,NV IT=LLL(LI,I) 18 TR(IT)=TR(IT)+FE(I)*A 10 CONTINUE 4 CONTINUE C DO 20 I=1,NY TU(I)=0.D0 DO 20 J=1,NY 20 TU(I)=TU(I)+GKK(I,J)*TR(J) DO 22 I=1,NY 22 TU(I)=TU(I)+FF(I) TU(NN)=0.D0 RETURN END C * * * * * SUBROUTINE STRESS C *** CALCULATE STRESSES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB3/STR(1690,6),EPS(1690,3) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) DIMENSION EU(8) C DO 10 L=1,NNE DO 20 I=1,6 20 STR(L,I)=0.D0 KOL2=K2(L) GO TO (11,12,11,14,14,14),KOL2 C *** TRUSS & TEXTILE 11 DO 22 I=1,4 IT=LLL(L,I) 22 EU(I)=TU(IT) E=EE(L) STR(L,1)=0.D0 DO 24 I=1,4 24 STR(L,1)=STR(L,1)+BMX(L,1,I)*EU(I)*E GO TO 10 C *** BEAM 12 DO 26 I=1,6 IT=LLL(L,I) 26 EU(I)=TU(IT) DO 28 I=1,6 STR(L,I)=0.D0 DO 28 J=1,6 28 STR(L,I)=STR(L,I)+DBM(L,I,J)*EU(J) GO TO 10 C *** PLANE-STRAIN AND INTERFACE 14 NV=8 IF(IJK(L,4).EQ.0) NV=6 DO 30 I=1,NV IT=LLL(L,I) 30 EU(I)=TU(IT) DO 32 I=1,3 EPS(L,I)=0.D0 DO 32 J=1,NV 32 EPS(L,I)=EPS(L,I)+BMX(L,I,J)*EU(J) DO 34 I=1,3 STR(L,I)=0.D0 DO 34 J=1,NV 34 STR(L,I)=STR(L,I)+DBM(L,I,J)*EU(J) 10 CONTINUE RETURN END C * * * * * SUBROUTINE YIELD(RST) C *** FIND YIELD ELEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB16/DEP(1690,3,3),GG1(1690,3,2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) DIMENSION RST(1690,3) C IL=0 IM=0 DO 10 L=1,NNE KOL2=K2(L) GO TO (10,10,10,14,14,16),KOL2 14 CONTINUE SX=RST(L,1) SY=RST(L,2) TA=RST(L,3) CALL MOHRCO(L,SX,SY,TA,F,S,P2,B0) GO TO 20 16 CONTINUE SG=RST(L,2) TA=RST(L,3) CALL COULOM(L,SG,TA,F,S,0) 20 CONTINUE PI2(L)=S IF(MFA(L).GE.1) GO TO 30 IF(F) 22,24,24 22 CONTINUE JFA(L)=0 GO TO 10 24 CONTINUE IF(JFA(L).GE.1) GO TO 30 IM=IM+1 IFA(IM)=L IFS(L)=ISTP 30 CONTINUE IL=IL+1 JFA(L)=IL KFA(IL)=L 10 CONTINUE NFA=IM NFAT=IL RETURN END C * * * * * SUBROUTINE NOTENS C *** FIND NO-TENSION ELEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB3/STR(1690,6),EPS(1690,3) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB20/NTS,ITS(1500),JTS(1690) DIMENSION STD(3),TQS(3),TST(3) C IT=0 DO 10 LI=1,NNE KOL2=K2(LI) GO TO (10,10,13,10,15,16),KOL2 13 CONTINUE SG=STR(LI,1) IF(SG.LE.0.) GO TO 10 GO TO 18 15 CONTINUE IF(JFA(LI).LE.0) GO TO 10 16 CONTINUE DO 20 I=1,3 STD(I)=0.D0 DO 20 J=1,3 20 STD(I)=STD(I)+TMI(LI,I,J)*(QST(LI,J)-STC(LI,J)) IF(STD(2).GE.0.) GO TO 10 18 IT=IT+1 ITS(IT)=LI JTS(LI)=IT 10 CONTINUE NTS=IT 80 RETURN END C * * * * * SUBROUTINE CALSTA C *** CALCULATE YIELD STRESSES SIGUMA-A IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) DIMENSION STD(3),RST(3),TST(3) C DO 10 IE=1,NFA LI=IFA(IE) DO 2 I=1,3 RST(I)=QST(LI,I) 2 TST(I)=PST(LI,I) DO 4 I=1,3 4 STD(I)=RST(I)-TST(I) C KOL2=K2(LI) GO TO (10,10,10,14,14,16),KOL2 14 CONTINUE SX=RST(1) SY=RST(2) TA=RST(3) CALL MOHRCO(LI,SX,SY,TA,F1,S,P2,B0) SX=TST(1) SY=TST(2) TA=TST(3) IF(ISTA.LE.0) GO TO 26 P1=(RST(1)+RST(2))*0.5 SX=P1 SY=P1 TA=0.D0 26 CONTINUE CALL MOHRCO(LI,SX,SY,TA,F0,S,P2,B0) IF(F0.LE.0.) GO TO 20 F1=F0 DO 22 I=1,3 RST(I)=PST(LI,I) 22 TST(I)=SST(LI,I) DO 24 I=1,3 24 STD(I)=RST(I)-TST(I) SX=TST(1) SY=TST(2) TA=TST(3) CALL MOHRCO(LI,SX,SY,TA,F0,S,P2,B0) 20 CONTINUE C1=-F0/(F1-F0) SX=TST(1)+STD(1)*C1 SY=TST(2)+STD(2)*C1 TA=TST(3)+STD(3)*C1 CALL MOHRCO(LI,SX,SY,TA,F2,S,P2,B0) B2=B0**(-0.5) SF=SFI(LI) A1=B2*P2-SF A2=B2*P2*(-1.)-SF A3=B2*4.*TA R=A1*STD(1)+A2*STD(2)+A3*STD(3) IF(R.EQ.0.) R=0.0001 C=C1-F2/R C IF(C.GT.1.) C=C1 GO TO 18 16 CONTINUE SG=RST(2) TA=RST(3) CALL COULOM(LI,SG,TA,F1,S,1) SG=TST(2) TA=TST(3) IF(ISTA.LE.0) GO TO 36 SG=RST(2) TA=0.D0 36 CONTINUE CALL COULOM(LI,SG,TA,F0,S,1) IF(F0.LE.0.) GO TO 30 F1=F0 DO 32 I=1,3 RST(I)=PST(LI,I) 32 TST(I)=SST(LI,I) DO 34 I=1,3 34 STD(I)=RST(I)-TST(I) SG=TST(2) TA=TST(3) CALL COULOM(LI,SG,TA,F0,S,1) 30 CONTINUE C=-F0/(F1-F0) 18 CONTINUE DO 50 I=1,3 50 STA(IE,I)=TST(I)+STD(I)*C C GO TO (10,10,10,54,54,56),KOL2 54 SX=STA(IE,1) SY=STA(IE,2) TA=STA(IE,3) CALL MOHRCO(LI,SX,SY,TA,F,S,P2,B0) GO TO 58 56 SG=STA(IE,2) TA=STA(IE,3) CALL COULOM(LI,SG,TA,F,S,0) 58 PI1(LI)=S 10 CONTINUE RETURN END C * * * * * SUBROUTINE CALSTB(IC) C *** CALCULATE CURRENT STRESSES SIGUMA-B IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB3/STR(1690,6),EPS(1690,3) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB16/DEP(1690,3,3),GG1(1690,3,2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB20/NTS,ITS(1500),JTS(1690) DIMENSION SBD(500,3),QSD(500,3),SAD(500,3),S0D(500,3) DIMENSION DP(3),STD(3),SCD(500,3),SRD(500,3) C DO 10 LI=1,NNE KOL2=K2(LI) GO TO (10,10,13,10,15,16),KOL2 13 CONTINUE IF(JTS(LI).LE.0) GO TO 10 STB(LI,1)=QST(LI,1)-ST0(LI,1) GO TO 10 15 CONTINUE IF(JFA(LI).LE.0) GO TO 10 16 CONTINUE IF(JFA(LI).GT.0) GO TO 18 IF(JTS(LI).LE.0) GO TO 10 18 CONTINUE DO 20 I=1,3 DP(I)=0.D0 DO 20 J=1,3 20 DP(I)=DP(I)+TMX(LI,I,J)*ST0(LI,J) DO 22 I=1,3 22 STB(LI,I)=QST(LI,I)-DP(I) 10 CONTINUE RETURN END C * * * * * SUBROUTINE CALSTC C *** CALCULATE PRECEDING STRESSES SIGUMA-C IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB3/STR(1690,6),EPS(1690,3) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) C DO 10 IL=1,NFAT LI=KFA(IL) DO 12 I=1,3 12 EPE(LI,I)=0.D0 C IS=IFS(LI)-ISTP IF(IS) 20,22,22 20 CONTINUE DO 14 I=1,3 14 STC(LI,I)=PST(LI,I) GO TO 10 22 CONTINUE DO 16 I=1,3 EPE(LI,I)=0.D0 DO 16 J=1,3 16 EPE(LI,I)=EPE(LI,I)+DIX(LI,I,J)*(STAX(LI,J)-PST(LI,J)) DO 18 I=1,3 18 STC(LI,I)=STAX(LI,I) 10 CONTINUE RETURN END C * * * * * SUBROUTINE DEPMAT C *** CALCULATE DEP-MATRIX IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB3/STR(1690,6),EPS(1690,3) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB16/DEP(1690,3,3),GG1(1690,3,2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB19/IWE(20),NRE,IRE(50),LRE(1690) DIMENSION T(3,3),D(3,3),TD(3,3),DPF(3,3) C DO 10 IE=1,NFA LI=IFA(IE) KOL2=K2(LI) C GO TO (10,10,10,10,15,16),KOL2 15 CONTINUE SI=1. IF(LRE(LI).GE.2) SI=-1. E=EE(LI) P=PP(LI) G=E/(2.*(1.+P)) C SIT=PI3(LI) TFI=SFI(LI)/CFI(LI) TDL=SDL(LI)/CDL(LI) ALF=3.14159*0.25+DATAN(TFI)*0.5 BET=(ALF+SIT)*(-1.) IF(LRE(LI).GE.2) BET=ALF-SIT PI4(LI)=BET IF(ICR.LE.0) GO TO 12 TFI=0.D0 TDL=0.D0 12 CONTINUE CB=DCOS(BET) SB=DSIN(BET) T(1,1)=CB*CB T(1,2)=SB*SB T(1,3)=-2.*SB*CB T(2,1)=T(1,2) T(2,2)=T(1,1) T(2,3)=2.*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(LI,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(LI,I,J)=T(I,J) DO 26 I=1,3 DO 26 J=1,3 26 T(I,J)=TMX(LI,I,J) C 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 CALL MULT(T,DPF,TD,3,3,3) CALL XULT(TD,T,DPF,3,3,3) GO TO 18 16 CONTINUE SI=-1. IF(STR(LI,3).LT.0.) SI=1. E=PP(LI) P=HH(LI) G=EE(LI) TFI=SFI(LI)/CFI(LI) TDL=SDL(LI)/CDL(LI) IF(ICR.LE.1) GO TO 30 TFI=0.D0 TDL=0.D0 30 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 44 I=1,3 DO 44 J=1,3 44 DEP(LI,I,J)=DPF(I,J) 10 CONTINUE RETURN END C * * * * * SUBROUTINE YGRAD1 C *** CALCULATE 'INITIAL STRESSES' (STEP-1) IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB16/DEP(1690,3,3),GG1(1690,3,2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB20/NTS,ITS(1500),JTS(1690) DIMENSION TR(2850),FE(8),GU(2850),BT(8,3) C DO 10 IE=1,NFA LI=IFA(IE) NV=8 IF(IJK(LI,4).EQ.0) NV=6 DO 12 I=1,NV DO 12 J=1,3 BT(I,J)=0.D0 DO 12 K=1,3 12 BT(I,J)=BT(I,J)+BMX(LI,K,I)*TMX(LI,K,J) A=AES(LI) C DO 20 II=2,3 DO 22 I=1,NY 22 TR(I)=0.D0 DO 24 I=1,NV 24 FE(I)=BT(I,II) DO 26 I=1,NV IT=LLL(LI,I) 26 TR(IT)=TR(IT)+FE(I)*A DO 28 I=1,NY GU(I)=0.D0 DO 28 J=1,NY 28 GU(I)=GU(I)+GKK(I,J)*TR(J) C DO 30 I=1,NY 30 GG1(LI,II,I)=GU(I) GG1(LI,II,NN)=0.D0 20 CONTINUE 10 CONTINUE RETURN END C * * * * * SUBROUTINE YGRAD2 C *** CALCULATE 'INITIAL STRESSES' (STEP-2) IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB16/DEP(1690,3,3),GG1(1690,3,2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB20/NTS,ITS(1500),JTS(1690) DIMENSION TR(2850),FE(8),GU(2850),BT(8,3) C IF(NTS.LE.0) GO TO 80 DO 10 IE=1,NTS LI=ITS(IE) IF(JFA(LI).GE.1) GO TO 10 KOL2=K2(LI) GO TO (10,10,13,10,10,16), KOL2 16 CONTINUE NV=8 DO 20 I=1,NV DO 20 J=1,3 BT(I,J)=0.D0 DO 20 K=1,3 20 BT(I,J)=BT(I,J)+BMX(LI,K,I)*TMX(LI,K,J) A=AES(LI) C II=2 DO 22 I=1,NY 22 TR(I)=0.D0 DO 24 I=1,NV 24 FE(I)=BT(I,II) DO 26 I=1,NV IT=LLL(LI,I) 26 TR(IT)=TR(IT)+FE(I)*A DO 28 I=1,NY GU(I)=0.D0 DO 28 J=1,NY 28 GU(I)=GU(I)+GKK(I,J)*TR(J) DO 30 I=1,NY 30 GG1(LI,II,I)=GU(I) GG1(LI,II,NN)=0.D0 GO TO 10 C 13 CONTINUE NV=4 DO 32 I=1,NV 32 BT(I,J)=BMX(LI,1,I) A=AES(LI) C II=1 DO 34 I=1,NY 34 TR(I)=0.D0 DO 36 I=1,NV 36 FE(I)=BT(I,II) DO 38 I=1,NV IT=LLL(LI,I) 38 TR(IT)=TR(IT)+FE(I)*A DO 40 I=1,NY GU(I)=0.D0 DO 40 J=1,NY 40 GU(I)=GU(I)+GKK(I,J)*TR(J) DO 42 I=1,NY 42 GG1(LI,II,I)=GU(I) GG1(LI,II,NN)=0.D0 EPE(LI,1)=0.D0 10 CONTINUE 80 CONTINUE RETURN END C * * * * * SUBROUTINE YGRAD3 C *** CALCULATE 'INITIAL STRESSES' (STEP-3) IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB16/DEP(1690,3,3),GG1(1690,3,2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) COMMON/LB20/NTS,ITS(1500),JTS(1690) DIMENSION EU(8),EP(3),FP(3),DM(3,3),RF(2150), * AA(2150,2150),SL0(2150) C DO 2 L=1,NNE DO 2 I=1,3 2 ST0(L,I)=0.D0 C NVS=NFAT+NTS IF(NVS.LE.NVM) GO TO 5 WRITE(6,200) NVS STOP 5 CONTINUE C DO 10 IL=1,NVS LI=KFA(IL) IF(IL.GT.NFAT) LI=ITS(IL-NFAT) KOL2=K2(LI) II=3 IF(IL.GT.NFAT) II=2 IF(IL.GT.NFAT.AND.KOL2.EQ.3) II=1 NVI=8 IF(IJK(LI,4).EQ.0) NVI=6 IF(KOL2.EQ.3) NVI=4 C DO 20 JL=1,NVS LJ=KFA(JL) IF(JL.GT.NFAT) LJ=ITS(JL-NFAT) KOL3=K2(LJ) IJ=3 IF(JL.GT.NFAT) IJ=2 IF(JL.GT.NFAT.AND.KOL3.EQ.3) IJ=1 DO 22 I=1,NVI IT=LLL(LI,I) 22 EU(I)=GG1(LJ,IJ,IT) NVST=3 IF(KOL2.EQ.3) NVST=1 DO 24 I=1,NVST EP(I)=0.D0 DO 24 J=1,NVI 24 EP(I)=EP(I)+BMX(LI,I,J)*EU(J) IF(IL-NFAT) 25,25,27 25 CONTINUE IF(JTS(LI).GE.1) GO TO 27 C *** WHEN SIGUMA-ST IS LESS THAN 0, CONSIDER TAU-ST AS 0. DO 26 I=1,3 DO 26 J=1,3 26 DM(I,J)=DMX(LI,I,J)-DEP(LI,I,J) GO TO 29 27 CONTINUE IF(KOL3.EQ.3) GO TO 40 DO 28 I=1,3 DO 28 J=1,3 28 DM(I,J)=DMX(LI,I,J) 29 CONTINUE DO 30 I=1,3 FP(I)=0.D0 DO 30 J=1,3 30 FP(I)=FP(I)+DM(I,J)*EP(J) DO 32 I=1,3 EP(I)=0.D0 DO 32 J=1,3 32 EP(I)=EP(I)+TMI(LI,I,J)*FP(J) GO TO 48 C 40 CONTINUE DM(1,1)=EE(LI) EP(1)=DM(1,1)*EP(1) C 48 C1=0.D0 IF(LJ.EQ.LI.AND.IJ.EQ.II) C1=1.D0 AA(IL,JL)=C1-EP(II) 20 CONTINUE 10 CONTINUE C DO 50 JL=1,NVS LJ=KFA(JL) IF(JL.GT.NFAT) LJ=ITS(JL-NFAT) KOL3=K2(LJ) NV=8 IF(IJK(LJ,4).EQ.0) NV=6 IF(KOL3.EQ.3) NV=4 DO 52 I=1,NV IT=LLL(LJ,I) 52 EU(I)=FF(IT) NVST=3 IF(KOL3.EQ.3) NVST=1 DO 54 I=1,NVST EP(I)=0.D0 DO 54 J=1,NV 54 EP(I)=EP(I)+BMX(LJ,I,J)*EU(J) IF(JL-NFAT) 55,55,57 55 CONTINUE IF(JTS(LJ).GE.1) GO TO 57 DO 56 I=1,3 DO 56 J=1,3 56 DM(I,J)=DMX(LJ,I,J)-DEP(LJ,I,J) GO TO 59 57 CONTINUE IF(KOL3.EQ.3) GO TO 64 DO 58 I=1,3 DO 58 J=1,3 58 DM(I,J)=DMX(LJ,I,J) 59 CONTINUE DO 60 I=1,3 FP(I)=0.D0 DO 60 J=1,3 60 FP(I)=FP(I)+DM(I,J)*EP(J) DO 62 I=1,3 EP(I)=0.D0 DO 62 J=1,3 62 EP(I)=EP(I)+TMI(LJ,I,J)*FP(J) GO TO 66 C 64 CONTINUE DM(1,1)=EE(LJ) EP(1)=DM(1,1)*EP(1) C 66 CONTINUE IJ=3 IF(JL.GT.NFAT) IJ=2 IF(KOL3.EQ.3) IJ=1 RF(JL)=EP(IJ) 50 CONTINUE C DO 70 JL=1,NVS LJ=KFA(JL) IF(JL.GT.NFAT) LJ=ITS(JL-NFAT) KOL3=K2(LJ) IF(JL-NFAT) 71,71,73 71 CONTINUE IF(JTS(LJ).GE.1) GO TO 73 DO 72 I=1,3 DO 72 J=1,3 72 DM(I,J)=DMX(LJ,I,J)-DEP(LJ,I,J) GO TO 75 73 CONTINUE IF(KOL3.EQ.3) GO TO 80 DO 74 I=1,3 DO 74 J=1,3 74 DM(I,J)=DMX(LJ,I,J) 75 CONTINUE DO 76 I=1,3 FP(I)=0.D0 DO 76 J=1,3 76 FP(I)=FP(I)+DM(I,J)*EPE(LJ,J) DO 78 I=1,3 EP(I)=0.D0 DO 78 J=1,3 78 EP(I)=EP(I)+TMI(LJ,I,J)*FP(J) GO TO 82 C 80 CONTINUE DM(1,1)=EE(LJ) EP(1)=DM(1,1)*EPE(LJ,1) C 82 IJ=3 IF(JL.GT.NFAT) IJ=2 IF(KOL3.EQ.3) IJ=1 RF(JL)=RF(JL)-EP(IJ) 70 CONTINUE C CALL GAUSSZ(AA,RF,NVM,NVS) C DO 84 JL=1,NVS LJ=KFA(JL) IF(JL.GT.NFAT) LJ=ITS(JL-NFAT) KOL3=K2(LJ) IJ=3 IF(JL.GT.NFAT) IJ=2 IF(KOL3.EQ.3) IJ=1 ST0(LJ,IJ)=RF(JL) 84 CONTINUE 200 FORMAT(/'* NUMBER OF INITIAL STRESSES IS OVER',I5) RETURN END C * * * * * SUBROUTINE MOHRCO(L,SX,SY,TA,F,S,P2,B0) C *** MOHR-COULOMB YIELD CONDITION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) P1=SX+SY P2=SX-SY A1=P1*SFI(L)+2.*CC(L)*CFI(L) B0=P2*P2+4.*TA*TA IF(B0.LT.0.0001) B0=0.0001 B1=DSQRT(B0) F=B1-A1 S=A1/B1 RETURN END C * * * * * SUBROUTINE COULOM(L,SG,TA,F,S,IC) C *** COULOMB YIELD CONDITION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) ST=CC(L)+SG*SFI(L)/CFI(L) IF(IC.LE.0) ST=VFUNC(ST) TA=ABS(TA) F=TA-ST S=0.D0 IF(TA.GT.0.001) S=ST/TA RETURN END C * * * * * SUBROUTINE TRUSS(L) C *** TRUSS ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) DIMENSION LL(4),BM(4),EK(4,4),T(4,4),TEK(4,4),TET(4,4) C II=IJK(L,1) JJ=IJK(L,2) C E=EE(L) A=PP(L) 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./EL BM(2)=0.D0 BM(3)=-1./EL BM(4)=0.D0 DO 2 I=1,4 DO 2 J=1,4 2 EK(I,J)=BM(I)*BM(J)*E*EL*A C DO 4 I=1,4 DO 4 J=1,4 4 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) C AES(L)=A*EL DO 6 J=1,4 BMX(L,1,J)=0.D0 DO 6 K=1,4 6 BMX(L,1,J)=BMX(L,1,J)+BM(K)*T(K,J) DBM(L,1,1)=E*EL LL(4)=3*JJ-1 LL(3)=LL(4)-1 LL(2)=3*II-1 LL(1)=LL(2)-1 DO 8 I=1,4 IT=LL(I) 8 LLL(L,I)=LY(IT) C NV=4 DO 10 I=1,NV IT=LLL(L,I) IF(IT.EQ.NN) GO TO 10 DO 12 J=1,NV JT=LLL(L,J) IF(JT.EQ.NN) GO TO 12 GKK(IT,JT)=GKK(IT,JT)+TET(I,J) 12 CONTINUE 10 CONTINUE RETURN END C * * * * * SUBROUTINE BEAM(L) C *** BEAM (RAHMEN) ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) DIMENSION LL(10),EK(6,6),T(6,6),TEK(6,6),TET(6,6),EKT(6,6) C II=IJK(L,1) JJ=IJK(L,2) C E=EE(L) A=PP(L) AI=HH(L) 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.*E*AI/EL G4=2.*G5 G3=3.*G5/EL G2=2.*G3/EL C DO 2 I=1,6 DO 2 J=1,6 T(I,J)=0.D0 EK(I,J)=0.D0 2 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 4 I=1,5 IP1=I+1 DO 6 J=IP1,6 6 EK(J,I)=EK(I,J) 4 CONTINUE DO 8 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. 8 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 10 I=1,6 DO 10 J=1,6 10 DBM(L,I,J)=EKT(I,J) C 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 12 I=1,6 IT=LL(I) 12 LLL(L,I)=LY(IT) C NV=6 DO 20 I=1,NV IT=LLL(L,I) IF(IT.EQ.NN) GO TO 20 DO 22 J=1,NV JT=LLL(L,J) IF(JT.EQ.NN) GO TO 22 GKK(IT,JT)=GKK(IT,JT)+TET(I,J) 22 CONTINUE 20 CONTINUE RETURN END C * * * * * SUBROUTINE PLANES(L,KOL2) C *** PLANE-STRAIN & INTERFACE ELEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB18/TMX(1690,3,3),TMI(1690,3,3) 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),TM(3,3) C II=IJK(L,1) JJ=IJK(L,2) KK=IJK(L,3) MM=IJK(L,4) C DO 2 I=1,3 DO 2 J=1,3 2 D(I,J)=0.D0 E=EE(L) P=PP(L) H=HH(L) GO TO (80,80,80,4,4,6),KOL2 4 GP=1.+P GM=1.-P GN=1.-2.*P G=E/(GP*GN) G1=GM*G G2=P*G G3=0.5*E/GP D(1,1)=G1 D(2,2)=G1 D(3,3)=G3 D(1,2)=G2 D(2,1)=G2 GO TO 8 6 GP=1.+H GM=1.-H GN=1.-2.*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) CB=DX/EL SB=DY/EL PI4(L)=DASIN(SB) C *** DEFINE TMX AND TMI FOR INTERFACE ELEMENT DO 10 I=1,3 DO 10 J=1,3 10 TM(I,J)=0.D0 DO 12 I=1,3 12 TM(I,I)=1.D0 DO 14 I=1,3 DO 14 J=1,3 14 TMX(L,I,J)=TM(I,J) DO 16 I=1,3 DO 16 J=1,3 16 TMI(L,I,J)=TM(I,J) C DO 18 I=1,8 DO 18 J=1,8 18 T(I,J)=0.D0 DO 20 I=1,7,2 T(I,I)=CB T(I,I+1)=SB T(I+1,I)=-SB T(I+1,I+1)=CB 20 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 22 I=1,8 COD(I)=0.D0 DO 22 J=1,8 22 COD(I)=COD(I)+T(I,J)*COR(J) DO 24 I=1,4 IE=2*I IO=IE-1 XXL(I)=COD(IO) YYL(I)=COD(IE) 24 CONTINUE C 8 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 IIT=IJK(L,I1) JJT=IJK(L,I2) KKT=IJK(L,I3) XI=XX(IIT) YI=YY(IIT) XJ=XX(JJT) YJ=YY(JJT) XK=XX(KKT) YK=YY(KKT) GO TO 48 46 XI=XXL(I1) YI=YYL(I1) XJ=XXL(I2) YJ=YYL(I2) XK=XXL(I3) YK=YYL(I3) 48 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.5 AE=AE+ABS(DA)*0.5 D2=0.125/DA IF(MM.EQ.0) D2=0.5/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./(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 IE=2*I IO=IE-1 B(1,IO)=B(3,IE) B(2,IE)=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 IE=IOE(I) DO 54 J=1,6 JE=IOE(J) ES(IE,JE)=ES(IE,JE)+EK(I,J) 54 CONTINUE 52 CONTINUE 40 CONTINUE C IF(MM.EQ.0) AE=ABS(DA) AES(L)=AE GO TO (80,80,80,64,64,66),KOL2 64 CONTINUE C1=0.25 IF(MM.EQ.0) C1=1. DO 70 I=1,NV DO 70 J=1,NV 70 ES(I,J)=ES(I,J)*AE*C1 DO 72 I=1,2 DO 72 J=1,NV 72 BMX(L,I,J)=0.D0 DO 74 J=1,NV 74 BMX(L,3,J)=BMT(J) N1=NV/2 DO 76 J=1,N1 JE=2*J JO=JE-1 BMX(L,1,JO)=BMX(L,3,JE) 76 BMX(L,2,JE)=BMX(L,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 ES(I,J)=ES(I,J)*AES(L)*0.25 DO 90 I=1,2 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(L,I,J)=0.D0 DO 96 K=1,8 96 BMX(L,I,J)=BMX(L,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(L,I)=LY(IT) DO 100 I=1,3 DO 100 J=1,NV DBM(L,I,J)=0.D0 DO 100 K=1,3 100 DBM(L,I,J)=DBM(L,I,J)+D(I,K)*BMX(L,K,J) DO 102 I=1,3 DO 102 J=1,3 102 DMX(L,I,J)=D(I,J) CALL DINV(D,3,3,0,DET,IND) DO 104 I=1,3 DO 104 J=1,3 104 DIX(L,I,J)=D(I,J) C DO 120 I=1,NV IT=LLL(L,I) IF(IT.EQ.NN) GO TO 120 DO 122 J=1,NV JT=LLL(L,J) IF(JT.EQ.NN) GO TO 122 GKK(IT,JT)=GKK(IT,JT)+ES(I,J) 122 CONTINUE 120 CONTINUE 80 RETURN END C * * * * * SUBROUTINE DINV(AA,N0,N1,N2,DET,IND) C *** INVERCE MATRIX IMPLICIT REAL*8(A-H,O-Z) DIMENSION AA(N0,N0),IPERM(2850),X(2850) 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.01-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 * * * * * SUBROUTINE GAUSSZ(A,X,N0,N) C *** GAUSS-ZAIDEL SWEEP OUT CALCULATION IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(N0,N0),X(N0),JUNJO(N0),SCALEF(N0),WORK(N0) DO 1 I=1,N 1 JUNJO(I)=I C DO 2 I=1,N P=DABS(A(1,1)) DO 3 J=2,N 3 P=DMAX1(P,ABS(A(I,J))) IF(P.EQ.0.) GO TO 9001 DO 4 J=1,N 4 A(I,J)=A(I,J)/P 2 X(I)=X(I)/P DO 5 J=1,N P=DABS(A(1,J)) DO 6 I=2,N 6 P=DMAX1(P,DABS(A(I,J))) IF(P.EQ.0.) GO TO 9002 DO 7 I=1,N 7 A(I,J)=A(I,J)/P 5 SCALEF(J)=P C DO 8 K=1,N-1 C P=DABS(A(K,K)) II=K JJ=K DO 9 I=K,N DO 10 J=K,N AA=DABS(A(I,J)) IF(AA.LE.P) GO TO 10 P=AA II=I JJ=J 10 CONTINUE 9 CONTINUE IF(P.LT.1.0E-8) GO TO 9003 DO 11 I=1,N W=A(I,K) A(I,K)=A(I,JJ) 11 A(I,JJ)=W DO 12 J=K,N W=A(K,J) A(K,J)=A(II,J) 12 A(II,J)=W W=X(K) X(K)=X(II) X(II)=W J=JUNJO(K) JUNJO(K)=JUNJO(JJ) JUNJO(JJ)=J C P=A(K,K) DO 13 J=K,N 13 A(K,J)=A(K,J)/P X(K)=X(K)/P DO 14 I=K+1,N Q=A(I,K) DO 15 J=K+1,N 15 A(I,J)=A(I,J)-Q*A(K,J) 14 X(I)=X(I)-Q*X(K) 8 CONTINUE X(N)=X(N)/A(N,N) C DO 16 L=2,N K=N-L+2 DO 16 I=1,K-1 16 X(I)=X(I)-X(K)*A(I,K) C DO 17 I=1,N 17 WORK(I)=X(I) DO 18 I=1,N J=JUNJO(I) 18 X(J)=WORK(I)/SCALEF(J) RETURN C 9001 WRITE(6,200) 200 FORMAT('I-TH ROW=0') GO TO 80 9002 WRITE(6,201) 201 FORMAT('I-TH COL=0') GO TO 80 9003 WRITE(6,202) 202 FORMAT('SMALL PIVOT') 80 STOP END C * * * * * FUNCTION VFUNC(C) C *** VFUNC(C)=0, WHEN C IS LESS THAN 0 IMPLICIT REAL*8(A-H,O-Z) VFUNC=0.D0 IF(C.GT.0.) VFUNC=C RETURN END C * * * * * SUBROUTINE PRINCE(SX,SY,TA,S1,S3,T) C *** CALCULATE 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.) 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) C *** MULTIPLY MATRIX A * B 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) C *** MULTIPLY MATRIX A(TRANSPOSE) * B 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) C *** MULTIPLY MATRIX A * B(TRANSPOSE) 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) C *** MULTIPLY MATRIX A * VECTOR U 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 OUTPUT(IC) C *** RESULT OUTPUT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB5/DBM(1690,6,8),BMX(1690,3,8),SDL(1690),CDL(1690) COMMON/LAB6/GKK(2850,2850),TF(2850),TU(2850),FF(2850) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LAB8/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LAB9/DMX(1690,3,3),DIX(1690,3,3) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB14/NSTP,MSTP,LSTP(300),ISTP,ISTA,NDI COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB19/IWE(20),NRE,IRE(50),LRE(1690) COMMON/LB20/NTS,ITS(1500),JTS(1690) COMMON/LB21/PS1(1690),PS3(1690),QS1(1690),QS3(1690),LFA(1690) COMMON/LB22/FIR,RKH,SCD,ERC,FLO(60),DWE(20) DIMENSION UX(1690),UY(1690),UM(1690),UZ(1690),WST(1690,7),EF(8), * TT(1690),YST(1500,6) DIMENSION IYIEL(1690),INOTE(1690),IP(1690) C IR=IC+5 QTU(NN)=0.D0 DO 2 L=1,NNE TT(L)=0.D0 IYIEL(L)=0 INOTE(L)=0 2 CONTINUE C *** DISPLACEMENTS DO 4 I=1,NNP JX=LY(3*I-2) JY=LY(3*I-1) JM=LY(3*I) SCU=SCD UX(I)=QTU(JX)*SCU UY(I)=QTU(JY)*SCU 4 UM(I)=QTU(JM)*SCU WRITE(IR,200) ISTP WRITE(IR,201) (UX(I),I=1,NNP) WRITE(IR,213) WRITE(IR,201) (UY(I),I=1,NNP) C IF(IC.LE.1) GO TO 80 R=FIR*FLOAT(ISTP)/FLOAT(NSTP) WRITE(8,220) ISTP,R,UY(NDI) WRITE(IR,202) C *** MICRO-AVS7 OUTPUT IM=LSTP(ISTP) IF(IM.LE.0) GO TO 7 WRITE(10,214) ISTP DO 6 I=1,NNP WRITE(10,221) I,UX(I),UY(I),UZ(I) 6 CONTINUE 7 CONTINUE C *** STRESSES DO 10 L=1,NNE KOL2=K2(L) A=1. IF(KOL2.EQ.1.OR.KOL2.EQ.3) A=PP(L) DO 12 I=1,3 12 WST(L,I)=STB(L,I)*A DO 14 I=4,7 14 WST(L,I)=0.D0 IF(KOL2.LE.3.OR.KOL2.GE.6) GO TO 10 SX=WST(L,1) SY=WST(L,2) TA=WST(L,3) CALL PRINCE(SX,SY,TA,S1,S3,TH) WST(L,4)=S1 WST(L,5)=S3 WST(L,6)=TH*180.D0/3.14159 WST(L,7)=(S1-S3)*0.5 10 CONTINUE DO 16 I=1,7 WRITE(IR,203) (WST(L,I),L=1,NNE) 16 WRITE(IR,213) WRITE(IR,204) (PI2(L),L=1,NNE) C *** YIELD ELEMENTS 80 IF(NFAT.EQ.0) GO TO 82 WRITE(IR,205) (KFA(K),K=1,NFAT) C WRITE(IR,205) (LFA(L),L=1,NNE) IF(IC.LE.1) GO TO 84 DO 30 K=1,NFAT L=KFA(K) IYIEL(L)=1 IP(K)=IFS(L) UX(K)=PI1(L) 30 CONTINUE C WRITE(IR,210) (IP(K),K=1,NFAT) C WRITE(IR,206) (UX(K),K=1,NFAT) WRITE(IR,207) DO 32 K=1,NFAT L=KFA(K) DO 32 I=1,3 32 YST(K,I)=ST0(L,I) DO 34 I=1,3 34 WRITE(IR,203) (YST(K,I),K=1,NFAT) WRITE(IR,208) DO 36 K=1,NFAT L=KFA(K) DO 36 I=1,3 36 YST(K,I)=STAX(L,I) DO 38 K=1,NFAT DO 40 I=4,6 40 YST(K,I)=0.D0 L=KFA(K) IF(K2(L).LE.3.OR.K2(L).GE.6) GO TO 38 SX=YST(K,1) SY=YST(K,2) TA=YST(K,3) CALL PRINCE(SX,SY,TA,S1,S3,TH) YST(K,4)=S1 YST(K,5)=S3 YST(K,6)=TH*180.D0/3.14159 38 CONTINUE DO 42 I=1,6 42 WRITE(IR,201) (YST(K,I),K=1,NFAT) C WRITE(IR,209) C DO 44 K=1,NFAT C DO 44 I=1,3 C 44 YST(K,I)=EPP(K,I)*100.D0 C DO 46 I=1,3 C 46 WRITE(IR,203) (YST(K,I),K=1,NFAT) C WRITE(IR,218) 84 CONTINUE C DO 50 K=1,NFAT C LI=KFA(K) C IP(K)=LRE(LI) C UX(K)=PI3(LI)*180.D0/3.14159 C 50 UY(K)=PI4(LI)*180.D0/3.14159 C WRITE(IR,203) (UX(K),K=1,NFAT) C WRITE(IR,203) (UY(K),K=1,NFAT) C WRITE(IR,219) (IP(K),K=1,NFAT) 82 CONTINUE C *** NO-TENSION ELEMENTS IF(NTS.LE.0) GO TO 85 WRITE(IR,212) (ITS(K),K=1,NTS) DO 51 I=1,NTS L=ITS(I) INOTE(L)=1 51 CONTINUE 85 CONTINUE C *** EQUIVALENT NODAL FORCES IF(NRE.LE.0) GO TO 86 WRITE(IR,216) DO 52 J=1,NRE L=IRE(J) NV=8 MM=IJK(L,4) IF(MM.LE.0) NV=6 NV1=NV/2 WRITE(IR,217) L,(IJK(L,I),I=1,NV1) DO 54 I=1,NV EF(I)=0.D0 DO 56 K=1,3 56 EF(I)=EF(I)+BMX(L,K,I)*QST(L,K)*AES(L) 54 CONTINUE WRITE(IR,215) (EF(I),I=1,NV) 52 CONTINUE 86 CONTINUE C *** RETAINING WALL ELEMENTS IF(NWE.EQ.0) GO TO 87 APN=0.D0 APS=0.D0 DO 60 K=1,NWE IW=IWE(K) SY=STB(IW,2)-PST(IW,2) WST(K,1)=VFUNC(SY) WST(K,2)=STB(IW,3)-PST(IW,3) APN=APN+WST(K,1)*DWE(K) APS=APS+WST(K,2)*DWE(K) 60 CONTINUE AP=APN**2+APS**2 AP=DSQRT(AP) WRITE(IR,211) AP,APN,APS WRITE(IR,232) (IWE(K),K=1,NWE) DO 62 I=1,2 62 WRITE(IR,203) (WST(K,I),K=1,NWE) 87 CONTINUE C *** MICRO-AVS7 OUTPUT IF(IC.LE.1) GO TO 90 IM=LSTP(ISTP) IF(IM.LE.0) GO TO 71 WRITE(10,222) WRITE(10,223) WRITE(10,224) WRITE(10,225) WRITE(10,226) WRITE(10,227) WRITE(10,228) WRITE(10,229) WRITE(10,230) DO 70 L=1,NNE 70 WRITE(10,231) L,(WST(L,I),I=1,5),TT(L),IYIEL(L),INOTE(L) 71 CONTINUE 90 CONTINUE C 200 FORMAT(/'LOADING STEP',I5/'DISP. X, Y') 201 FORMAT(10F8.4) 202 FORMAT('STRESS SX,SY,TA,S1,S3,TH') 203 FORMAT(10F8.3) 204 FORMAT('SAFETY-FACTOR'/(10F7.3)) 205 FORMAT('YIELD-ELEMENT'/(10I5)) 206 FORMAT('YIELD-FACTOR'/(10F7.3)) 207 FORMAT('INITIAL-STRESSES') 208 FORMAT('YIELD-STRESSES SX,SY,TA,S1,S3,TH') 209 FORMAT('PLASTIC STRAINS %') 210 FORMAT('FAILURE STEP'/(10I5)) 211 FORMAT('THRUST R,N,S',3F10.3) 212 FORMAT('NO-TENSION EL.'/(10I5)) 213 FORMAT(' ') 214 FORMAT(/'*** STEP',I5) 215 FORMAT(7E10.3) 216 FORMAT('EQUIVALENT NODAL FORCE'/' EL NODES') 217 FORMAT(16I4) 218 FORMAT('ANGLE OF S1, SLIP-SUR. DISP') 219 FORMAT('LEFT-1, RIGHT-2'/(10I5)) 220 FORMAT(I5,2F10.4) 221 FORMAT(I5,3E12.4) 222 FORMAT('8 1 1 1 1 1 1 1 1') 223 FORMAT('SX,') 224 FORMAT('SY,') 225 FORMAT('TA,') 226 FORMAT('S1,') 227 FORMAT('S3,') 228 FORMAT('TT,') 229 FORMAT('YIELD,') 230 FORMAT('NOTEN,') 231 FORMAT(I4,6E11.3,2I2) 232 FORMAT(10I5) 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/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB19/IWE(20),NRE,IRE(50),LRE(1690) COMMON/LB20/NTS,ITS(1500),JTS(1690) COMMON/LB22/FIR,RKH,SCD,ERC,FLO(60),DWE(20) DIMENSION SX(1690),SY(1690),XXX(1580),YYY(1580) DIMENSION J(1690) C DO 2 L=1,NNE J(L)=1 SX(L)=0.D0 SY(L)=0.D0 2 CONTINUE C DO 4 K=1,NFAT LI=KFA(K) SI=PI4(LI) SX(LI)=DCOS(SI) SY(LI)=DSIN(SI) 4 CONTINUE C NNE1=0 DO 6 L=1,NNE IF(IJK(L,3).EQ.0) GO TO 6 NNE1=NNE1+1 6 CONTINUE WRITE(11,301) NNP,NNE1 C DO 8 I=1,NNP 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,NNE 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,NNE 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/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LB10/LX(5000),LY(5000),NDF(5000),LLL(1690,8) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB21/PS1(1690),PS3(1690),QS1(1690),QS3(1690),LFA(1690) COMMON/LB22/FIR,RKH,SCD,ERC,FLO(60),DWE(20) DIMENSION UX(1580),UY(1580),UM(1580),WST(1690,6),PI9(1690) C QTU(NN)=0.D0 DO 2 I=1,NNP JX=LY(3*I-2) JY=LY(3*I-1) JM=LY(3*I) SCU=SCD UX(I)=QTU(JX)*SCU UY(I)=QTU(JY)*SCU 2 CONTINUE WRITE(12,310) 1 WRITE(12,311) WRITE(12,300) 1 C WRITE(12,301) NNP,NNE DO 10 I=1,NNP 10 WRITE(12,302) I,XX(I),YY(I),0 DO 12 L=1,NNE 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,NNP 14 WRITE(12,309) I,UX(I),UY(I),0 WRITE(12,313) WRITE(12,314) DO 20 L=1,NNE 20 PI9(L)=0.1D0 DO 22 L=1,NNE IF(JFA(L).GE.1) PI9(L)=2.D0 22 CONTINUE DO 24 L=1,NNE 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/NNP,NNE,NFL,NWE,NIS,ITER,NTRY,NY,NN,NFA,NFAT,ICR,NVM COMMON/LAB2/XX(1580),YY(1580),ZZ(1580),IX(1580),IY(1580),IQ(1580) COMMON/LAB4/K2(1690),IJK(1690,4),NBA,LBA(1690),LOI(60) COMMON/LAB7/AES(1690),EE(1690),PP(1690),HH(1690),PW(1690,2) COMMON/LB11/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB12/ST0(1690,3),STA(1500,3),STAX(1690,3),SST(1690,3) COMMON/LB13/STB(1690,3),STC(1690,3),EPE(1690,3),EPP(1500,3) COMMON/LB15/QST(1690,3),PST(1690,3),QTU(2850),PTU(2850) COMMON/LB17/IFA(1500),JFA(1690),KFA(1500),MFA(1690),IFS(1690) COMMON/LB21/PS1(1690),PS3(1690),QS1(1690),QS3(1690),LFA(1690) COMMON/LB22/FIR,RKH,SCD,ERC,FLO(60),DWE(20) DIMENSION WST(6),AX(1690),AY(1690),XXX(1580),YYY(1580) DIMENSION J(1690) C DO 2 L=1,NNE J(L)=0 AX(L)=0.D0 AY(L)=0.D0 2 CONTINUE C DO 10 L=1,NNE KOL2=K2(L) IF(KOL2.LE.3) GO TO 10 IF(JFA(L).GE.1) GO TO 12 DO 20 I=1,3 20 WST(I)=STB(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,NNE IF(IJK(L,3).EQ.0) GO TO 6 NNE1=NNE1+1 6 CONTINUE WRITE(13,301) NNP,NNE1 C DO 30 I=1,NNP 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,NNE 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,NNE 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