PROGRAM LSFEA2 C LIMIT STATE FINITE ELEMENT ANALYSIS (TYPE-2) C SLOPE STABILITY, AND EARTH PRESSURE PROBLEMS C ORIGINAL INITIAL STRESS METHOD C PORE-WATER PRESSURE C IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB2/ISB,IDL,JBA,ICR,NRE,ISTP,MSTP,NDI,FIR COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB13/ST0(1690,3),S0P(1690,3),STA(960,3),STAX(1690,3) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB18/DWE(20),ERM,ERG,ERC,SCD COMMON/LB19/PS1(1690),PS3(1690),QS1(1690),QS3(1690) COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) C CALL INPUTS CALL PREPAR(IND) IF(IND.NE.1) GO TO 10 CALL INISTR 10 CONTINUE CLOSE(6) CLOSE(8) CLOSE(9) 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,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB2/ISB,IDL,JBA,ICR,NRE,ISTP,MSTP,NDI,FIR COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB18/DWE(20),ERM,ERG,ERC,SCD COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) DIMENSION IDATA(12),IDATAQ(12),PNA(20),FAI(1690) * ,DLT(1690),STII(50,3),TY(1690) DIMENSION IBA(1690),IST(1690),NET(6),NM(1690),MAT(1690) DIMENSION E1(20),P1(20),H1(20),R1(20),C1(20),F1(20),D1(20) 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(11,FILE='M-AVS1.INP',STATUS='UNKNOWN') OPEN(12,FILE='M-AVS2.INP',STATUS='UNKNOWN') OPEN(13,FILE='M-AVS3.INP',STATUS='UNKNOWN') C NN: MAXIMUM NUMBER OF DISPLACEMENT VARIABLES NN=3110 DO 2 I=1,12 2 IDATAQ(I)=0.D0 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 NTR: NUMBER OF ITERATIONS C NCY=3: BEARING CAPACITY PROBLEM, C 5: EARTH PRESSURE & SLOPE STABILITY PROBLEMS C ISB=0:NOT CONSIDER SHEAR BAND C 1: INPUT LRE(L) AT EACH ELEMENT C LRE(L): TYPE OF SHEAR BAND DIRECTION (1 OR 2) C SLIP SURFACE DIRECTION ALF=PAI/4+FAI/2 C JBA=0: NOT CONSIDER HORIZONTAL ACCERALATION, 1:CONSIDER C ICR=0: NOT CONSIDER CRITICAL STATE C 1: CONSIDER FOR CONSTANT CONFINING PRESSURE C 2: CONSIDER AFTER FIRST YIELDING (PERFECTLY PLASTIC) C IAS=0: ASSOCIATED FLOW RULE, 1:NON-ASSOCIATED 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 READ(5,101) NNP,NNE,NMAT,NFL,NWE,NIS,ITER,NTR,NCY READ(5,101) ISB,JBA,ICR,IAS,IPR,IPW READ(5,101) NSTP,MSTP IF(NTR.LE.0) NTR=1 IF(NSTP.LE.0) NSTP=1 IF(MSTP.LE.0) MSTP=NSTP WRITE(6,200) NNP,NNE,NMAT,NFL,NWE,NIS,ITER,NTR,NCY, * ISB,JBA,ICR,IAS,IPR,IPW, * NSTP,MSTP 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): MATERIAL 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,5 IF(IDATA(I).EQ.0.) IDATA(I)=IDATAQ(I) IF(IDATA(I).LT.0.) IDATA(I)=0 IDATAQ(I)=IDATA(I) 14 CONTINUE K2(L)= IDATA(1) MAT(L)=IDATA(2) IBA(L)=IDATA(3) IST(L)=IDATA(4) LRE(L)=IDATA(5) 12 CONTINUE C *** MATERAL PARAMETERS C K: SET NO. C EE(L): YOUNG'S MODULUS FOR K2=1-5, SHEAR RIGIDITY FOR K2=6 C PP(L): POISSON'S RATIO FOR K2=4 & 5, C CROSS AREA FOR K2=1-3, YOUNG'S MODULUS FOR K2=6 C HH(L): 1.0 FOR K2=4 & 5, MOMENT OF INERTIA FOR K2=3 C POISSON'S RATIO FOR K2=6 C ROU(L): DENSITY C CC(L): COHESION, FAI(L): ANGLE OF SHEAR RESISTANCE C DLT(L): ANGLE OF DILATANCY WRITE(6,213) DO 16 M=1,NMAT READ(5,107) K,E1(M),P1(M),H1(M),R1(M),C1(M),F1(M),D1(M) WRITE(6,107) K,E1(M),P1(M),H1(M),R1(M),C1(M),F1(M),D1(M) 16 CONTINUE C DO 18 L=1,NNE M=MAT(L) EE(L)= E1(M) PP(L)= P1(M) HH(L)= H1(M) ROU(L)=R1(M) CC(L)= C1(M) CC0(L)=CC(L) FAI(L)=F1(M) DLT(L)=D1(M) FI=FAI(L)*3.14159D0/180.D0 SFI(L)=DSIN(FI) CFI(L)=DCOS(FI) DL=DLT(L)*3.14159D0/180.D0 C FLOW RULE IF(IAS.LE.0) DL=FI SDL(L)=DSIN(DL) CDL(L)=DCOS(DL) N=K2(L) NET(N)=NET(N)+1 IF(N.EQ.6) ROU(L)=0.D0 18 CONTINUE C *** ELEMENT TYPE WRITE(6,202) (NET(I),I=1,6) C *** EMBANKMENT ELEMENTS DO 20 N=1,NSTP 20 NBA(N)=0 DO 22 L=1,NNE IB=IBA(L) IF(IB.EQ.0) GO TO 22 NBA(IB)=NBA(IB)+1 LB=NBA(IB) LBA(IB,LB)=L 22 CONTINUE WRITE(6,203) DO 24 N=1,NSTP NB=NBA(N) IF(NB.EQ.0) GO TO 24 WRITE(6,101) N,NB WRITE(6,101) (LBA(N,I),I=1,NB) 24 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-DIREC. 2:Y-DIREC. 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,6 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) WRITE(6,107) K1,(STII(K,I),I=1,3) 44 CONTINUE 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(N),N=1,NSTP) WRITE(6,204) (RKH(N),N=1,NSTP) C *** CONVERGENCE CONSTANTS ETC. C SCD: SCALING FACTOR FOR DISPLACEMENT OUTPUT C ERM: CONV. CONST. FOR K2=5 & 6 (=1.0KPA) C ERG: CONV. CONST. FOR K2=3 C ERC: CONV. CONST. FOR CONFINING PRESSURE C FSS: SAFETY FACTOR C READ(5,106) SCD,ERM,ERG,ERC,FSS IF(SCD.LT.0.01D0) SCD=1.D0 IF(FSS.LE.0.D0) FSS=1.D0 WRITE(6,208) SCD,ERM,ERG,ERC,FSS C C *** DEVIDE BY SAFETY FACTOR FSS DO 56 L=1,NNE CC(L)=CC(L)/FSS CC0(L)=CC(L) 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 *** MEMORY OF INITIAL STRESSES DO 88 L=1,NNE DO 88 I=1,3 88 PS0(L,I)=PST(L,I) C C *** MICRO-AVS 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(16I5) 102 FORMAT(I5,2F10.3,3I5) 103 FORMAT(10I5) 104 FORMAT(2I5,5F10.3/10X,5F10.3) 105 FORMAT(3I5,4E10.3,2I5) 106 FORMAT(7E10.3) 107 FORMAT(I5,7E10.3) 108 FORMAT(10F8.3) 200 FORMAT(/'NODES=',I4,' ELEMENTS=',I4,' MATERIALS=',I5 */' LOADS=',I3,' WALL-ELEMENTS=',I3,' INITIAL-STRESSES=',I3 */' ITERATION=',I2,' TRYALS=',I5,' D-EP MATRIX STEPS=',I3 */' SHEAR-BAND=',I2,' HORIZONTAL ACCERARATION=',I2 */' CRITICAL-STATE='I2,' FLOW-RULE=',I2 */' PRINT=',I2,' PORE-WATER PRESSURE=',I2 */' LOADING-STEPS=',I5,' ACTUAL STEPS=',I5) 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 COEFCIENT'/(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 */' CONVERGENCE CONSTANT (K2=5, 6)=',F10.3 */' CONVERGENCE CONSTANT (K2=3)=',F10.3 */' CONVERGENCE CONSTANT FOR CONFINING PRESSURE=',F10.3 */' 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,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB2/ISB,IDL,JBA,ICR,NRE,ISTP,MSTP,NDI,FIR COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) DIMENSION EML(1680),LX(5050) C C *** ARRANGE DISPLACEMENT VECTOR 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 20 L=1,NNP LY(3*L-2)=IX(L) LY(3*L-1)=IY(L) IF(NDF(L).EQ.2) IQ(L)=1 20 LY(3*L)=IQ(L) NX=3*NNP 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 NY=LOC WRITE(6,200) NY IF(NY.LT.NN) GO TO 55 WRITE(*,201) IND=10 RETURN C *** GLOBAL STIFFNESS MATRIX & INVERSE MATRIX 55 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 CALL GSTIFF CALL DINV(GKK,NN,NY,0,DET,IND) IF(IND.EQ.1) GO TO 50 WRITE(*,202) IND C *** NODAL LOADS 50 DO 52 I=1,NY 52 TL(I)=0.D0 TL(NN)=0.D0 IF(NFL.EQ.0) GO TO 60 DO 54 I=1,NFL IT=LOI(I) IT=LY(IT) 54 TL(IT)=TL(IT)+FLO(I)/DFLOAT(NSTP) C *** EMBANKMENT LOADS FOR SEISMIC FORCE 60 CONTINUE WRITE(6,203) NB=NBA(1) IF(NB.LE.0) GO TO 80 DO 62 I=1,NNP 62 EM0(I)=0.D0 IF(JBA.LE.0) GO TO 80 DO 64 I=1,NB LB=LBA(1,I) KLB=K2(LB) GO TO (71,71,71,74,74,64),KLB 71 N1=2 C1=0.5D0 EMI=AES(LB)*ROU(LB)*C1 GO TO 68 74 MM=IJK(LB,4) N1=4 C1=0.25D0 IF(MM.EQ.0) N1=3 IF(MM.EQ.0) C1=0.333D0 EMI=AES(LB)*ROU(LB)*C1 68 DO 66 J=1,N1 II=IJK(LB,J) 66 EM0(II)=EM0(II)+EMI 64 CONTINUE WRITE(6,204) (EM0(I),I=1,NNP) 80 CONTINUE C 200 FORMAT(/'NUMBER OF VARIABLES=',I4) 201 FORMAT(/'* VARIABLE MEMORY=OVER') 202 FORMAT(/'* INDEX=',I3) 203 FORMAT(/'EMBANKMENT LOADS FOR SEISMIC FORCE') 204 FORMAT(10F8.2) RETURN END C * * * * * SUBROUTINE GSTIFF C *** CALCULATE GLOBAL STIFFNES MATRIX IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) C DO 2 I=1,NY DO 2 J=1,NY 2 GKK(I,J)=0.D0 C DO 10 L=1,NNE KOL2=K2(L) GO TO (11,12,11,14,14,16),KOL2 11 NV=4 GO TO 18 12 NV=6 GO TO 18 14 NV=8 IF(IJK(L,4).EQ.0) NV=6 GO TO 18 16 NV=8 18 CONTINUE 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)+EKK(L,I,J) 22 CONTINUE 20 CONTINUE 10 CONTINUE RETURN END C * * * * * SUBROUTINE INISTR C *** APPLY INITIAL STRESS METHOD IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB2/ISB,IDL,JBA,ICR,NRE,ISTP,MSTP,NDI,FIR COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB13/ST0(1690,3),S0P(1690,3),STA(960,3),STAX(1690,3) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB18/DWE(20),ERM,ERG,ERC,SCD COMMON/LB19/PS1(1690),PS3(1690),QS1(1690),QS3(1690) COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) C DO 2 L=1,NNE JFA(L)=0 LFA(L)=0 JTS(L)=0 IFS(L)=0 MFA(L)=0 PS3(L)=0.D0 QS3(L)=0.D0 DO 4 I=1,3 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 8 I=1,3 ST0(L,I)=0.D0 8 S0P(L,I)=0.D0 6 CONTINUE C CALL LOADVE ITR=0 DO 110 I1=1,NTR ITR=ITR+1 CALL DISPLA CALL STRESS DO 20 I=1,NY 20 QTU(I)=PTU(I)+TU(I) DO 26 L=1,NNE DO 22 I=1,6 22 QST(L,I)=PST(L,I)+STR(L,I) IF(K2(L).LE.3) GO TO 26 DO 24 I=1,2 24 QST(L,I)=QST(L,I)-PW(L,I) 26 CONTINUE DO 28 L=1,NNE DO 28 I=1,3 28 STB(L,I)=QST(L,I) CALL YIELD IF(ITER.EQ.0) GO TO 60 IF(NFA.LE.0) GO TO 30 CALL CALSTA DO 32 IE=1,NFA LI=IFA(IE) DO 34 I=1,3 34 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 32 CONTINUE 30 IF(NFAT.GE.1) CALL DEPMAT CALL NOTENS DO 36 L=1,NNE DO 38 I=1,3 38 ST0(L,I)=0.D0 IF(JFA(L).LE.0.AND.JTS(L).LE.0) GO TO 36 DO 40 I=1,3 40 ST0(L,I)=QST(L,I)-STB(L,I) IF(K2(L).EQ.6) ST0(L,1)=0.D0 36 CONTINUE C RMAX=0.D0 DO 50 L=1,NNE IF(K2(L).EQ.3) GO TO 50 DO 52 I=1,3 DI=ST0(L,I)-S0P(L,I) DI=DABS(DI) IF(DI.LE.RMAX) GO TO 52 RMAX=DI LIX=L IXY=I 52 CONTINUE 50 CONTINUE RMAG=0.D0 LIG=0 DO 54 L=1,NNE IF(K2(L).NE.3) GO TO 54 DI=ST0(L,1)-S0P(L,1) DI=DABS(DI) IF(DI.LE.RMAG) GO TO 54 RMAG=DI LIG=L 54 CONTINUE DO 56 L=1,NNE DO 56 I=1,3 56 S0P(L,I)=ST0(L,I) WRITE(6,202) ITR,LIX,IXY,RMAX,LIG,RMAG IF(RMAX.LE.ERM.AND.RMAG.LE.ERG) GO TO 60 110 CONTINUE C OPEN(7,FILE='PSLSFE',STATUS='UNKNOWN') CALL OUTPUT(2) CLOSE(7) CALL MAVS1 CALL MAVS2 CALL MAVS3 RETURN C 60 CALL OUTPUT(1) OPEN(7,FILE='PSLSFE',STATUS='UNKNOWN') CALL OUTPUT(2) CLOSE(7) C IF(ICR.NE.1) GO TO 70 IF(NFAT.LE.0) GO TO 70 DO 62 IL=1,NFAT LI=KFA(IL) IF(LFA(LI).GE.1) GO TO 62 SX=STB(LI,1) SY=STB(LI,2) TA=STB(LI,3) CALL PRINCE(SX,SY,TA,S1,S3,TH) QS1(LI)=S1 QS3(LI)=S3 DI=QS3(LI)-PS3(LI) IF(DABS(DI).GT.ERC) GO TO 62 LFA(LI)=1 CC(LI)=(S1-S3)*0.5D0 SFI(LI)=0.D0 CFI(LI)=1.D0 SDL(LI)=0.D0 CDL(LI)=1.D0 62 CONTINUE 70 DO 64 L=1,NNE PS3(L)=QS3(L) DO 64 I=1,3 64 PST(L,I)=STB(L,I) DO 66 I=1,NY 66 PTU(I)=QTU(I) DO 68 L=1,NNE 68 MFA(L)=JFA(L) C CALL MAVS1 CALL MAVS2 CALL MAVS3 100 CONTINUE C 200 FORMAT(/'* LOADING STEP=',I3) 201 FORMAT('STAGE',I3) 202 FORMAT(I3,' EL',I4,' X-Y',I2,F10.4,' EL-G',I4,F10.4) RETURN END C * * * * * SUBROUTINE LOADVE C *** CALCULATE LOAD VECTOR IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB2/ISB,IDL,JBA,ICR,NRE,ISTP,MSTP,NDI,FIR COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB18/DWE(20),ERM,ERG,ERC,SCD COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) DIMENSION EML(1680) C ITP=ISTP DO 2 I=1,NY 2 TF(I)=TL(I) DO 4 I=1,NNP 4 EML(I)=0.D0 IF(JBA.GE.1) GO TO 38 C IF(ITP.EQ.MSTP) GO TO 18 ITP1=ITP+1 DO 10 IL=ITP1,MSTP NB=NBA(IL) IF(NB.EQ.0) GO TO 10 DO 12 I=1,NB LB=LBA(IL,I) DO 14 J=1,3 14 PST(LB,J)=PS0(LB,J) 12 CONTINUE 10 CONTINUE C 18 NB=NBA(ITP) IF(NB.EQ.0) GO TO 38 DO 20 I=1,NB LB=LBA(ITP,I) CC(LB)=CC0(LB) IF(K2(LB).LE.2) GO TO 20 DO 22 J=1,3 22 PST(LB,J)=PS0(LB,J) 20 CONTINUE C DO 30 I=1,NB LB=LBA(ITP,I) KLB=K2(LB) GO TO (71,71,71,74,74,30),KLB 71 N1=2 C1=0.5D0 EMI=-AES(LB)*ROU(LB)*C1 GO TO 68 74 MM=IJK(LB,4) N1=4 C1=0.25D0 IF(MM.EQ.0) N1=3 IF(MM.EQ.0) C1=0.333D0 EMI=-AES(LB)*ROU(LB)*C1 68 DO 32 J=1,N1 II=IJK(LB,J) 32 EML(II)=EML(II)+EMI 30 CONTINUE C 38 DO 40 L=1,NNP IT=LY(3*L-2) JT=LY(3*L-1) TF(IT)=TF(IT)+EM0(L)*RKH(ITP) 40 TF(JT)=TF(JT)+EML(L) WRITE(6,200) (EML(I),I=1,NNP) 200 FORMAT(/'EMBANKMENT LOADS'/(10F8.2)) RETURN END C * * * * * SUBROUTINE DISPLA C *** CALCULATE DISPLACEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB13/ST0(1690,3),S0P(1690,3),STA(960,3),STAX(1690,3) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) DIMENSION TR(3110),FE(8) C DO 2 I=1,NY 2 TR(I)=0.D0 DO 10 L=1,NNE IF(JFA(L).LE.0.AND.JTS(L).LE.0) GO TO 10 KOL2=K2(L) GO TO (10,10,13,10,15,15),KOL2 13 NV=4 DO 20 I=1,NV 20 FE(I)=BMX(L,1,I)*ST0(L,1) GO TO 18 15 NV=8 IF(IJK(L,4).EQ.0) NV=6 DO 22 I=1,NV FE(I)=0.D0 DO 22 J=1,3 22 FE(I)=FE(I)+BMX(L,J,I)*ST0(L,J) 18 A=AES(L) DO 24 I=1,NV IT=LLL(L,I) 24 TR(IT)=TR(IT)+FE(I)*A 10 CONTINUE DO 30 I=1,NY 30 TR(I)=TF(I)+TR(I) DO 32 I=1,NY TU(I)=0.D0 DO 32 J=1,NY 32 TU(I)=TU(I)+GKK(I,J)*TR(J) TU(NN)=0.D0 RETURN END C * * * * * SUBROUTINE STRESS C *** CALCULATE STRESSES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) 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 C *** FIND YIELD ELEMENTS IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB2/ISB,IDL,JBA,ICR,NRE,ISTP,MSTP,NDI,FIR COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB13/ST0(1690,3),S0P(1690,3),STA(960,3),STAX(1690,3) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) C IL=0 IM=0 DO 10 L=1,NNE KOL2=K2(L) GO TO (10,10,10,14,14,16),KOL2 14 SX=STB(L,1) SY=STB(L,2) TA=STB(L,3) CALL MOHRCO(L,SX,SY,TA,F,S,P2,B0) GO TO 20 16 SG=STB(L,2) TA=STB(L,3) CALL COULOM(L,SG,TA,F,S,0) 20 PI2(L)=S IF(MFA(L).GE.1) GO TO 30 IF(F) 22,24,24 22 JFA(L)=0 GO TO 10 24 IF(JFA(L).GE.1) GO TO 30 AMAX=0. DO 26 I=1,3 S=STAX(L,I) AS=ABS(S) IF(AS.GT.AMAX) AMAX=AS 26 CONTINUE IF(AMAX.GT.0.00001D0) GO TO 30 IM=IM+1 IFA(IM)=L IFS(L)=ISTP 30 IL=IL+1 JFA(L)=IL KFA(IL)=L 10 CONTINUE NFA=IM NFAT=IL RETURN END C * * * * * SUBROUTINE CALSTA C *** CALCULATE YIELD STRESSES IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB13/ST0(1690,3),S0P(1690,3),STA(960,3),STAX(1690,3) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) DIMENSION STD(3) C DO 10 IE=1,NFA LI=IFA(IE) DO 20 I=1,3 20 STD(I)=STB(LI,I)-PST(LI,I) KOL2=K2(LI) GO TO (10,10,10,14,14,16),KOL2 14 SX=STB(LI,1) SY=STB(LI,2) TA=STB(LI,3) CALL MOHRCO(LI,SX,SY,TA,F1,S,P2,B0) SX=PST(LI,1) SY=PST(LI,2) TA=PST(LI,3) CALL MOHRCO(LI,SX,SY,TA,F0,S,P2,B0) C1=-F0/(F1-F0) SX=PST(LI,1)+STD(1)*C1 SY=PST(LI,2)+STD(2)*C1 TA=PST(LI,3)+STD(3)*C1 CALL MOHRCO(LI,SX,SY,TA,F2,S,P2,B0) B2=B0**(-0.5D0) SF=SFI(LI) 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=STB(LI,2) TA=STB(LI,3) CALL COULOM(LI,SG,TA,F1,S,1) SG=PST(LI,2) TA=PST(LI,3) CALL COULOM(LI,SG,TA,F0,S,1) C=-F0/(F1-F0) 18 DO 22 I=1,3 22 STA(IE,I)=PST(LI,I)+STD(I)*C C GO TO (10,10,10,34,34,36),KOL2 34 SX=STA(IE,1) SY=STA(IE,2) TA=STA(IE,3) CALL MOHRCO(LI,SX,SY,TA,F,S,P2,B0) GO TO 38 36 SG=STA(IE,2) TA=STA(IE,3) CALL COULOM(LI,SG,TA,F,S,0) 38 PI1(LI)=S 10 CONTINUE RETURN END C * * * * * SUBROUTINE DEPMAT C *** CALCULATE DEP MATRIX IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB2/ISB,IDL,JBA,ICR,NRE,ISTP,MSTP,NDI,FIR COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB13/ST0(1690,3),S0P(1690,3),STA(960,3),STAX(1690,3) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB18/DWE(20),ERM,ERG,ERC,SCD COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) DIMENSION FD(3),DEP(3,3),T(3,3) * ,EPE(3),PSA(3),FSA(3),D(3,3),DP(3),DPF(3,3),DPFD(3,3) C DO 10 IL=1,NFAT LI=KFA(IL) IS=IFS(LI)-ISTP IF(IS) 12,14,14 12 DO 2 I=1,3 EPP(IL,I)=EPS(LI,I) 2 STC(LI,I)=PST(LI,I) GO TO 10 14 DO 4 I=1,3 EPE(I)=0.D0 DO 4 J=1,3 4 EPE(I)=EPE(I)+DIX(LI,I,J)*(STAX(LI,J)-PST(LI,J)) DO 6 I=1,3 EPP(IL,I)=EPS(LI,I)-EPE(I) 6 STC(LI,I)=STAX(LI,I) 10 CONTINUE C DO 20 IL=1,NFAT LI=KFA(IL) IF(ICR.LE.1) GO TO 18 SFI(LI)=0.D0 CFI(LI)=1.D0 SDL(LI)=0.D0 CDL(LI)=1.D0 18 CONTINUE NCY1=NCY IF(ISB.GE.1) NCY1=1 DO 22 NC=1,NCY1 KOL2=K2(LI) GO TO (20,20,20,24,24,26),KOL2 24 ISB1=ISB+1 GO TO (30,32,32,32),ISB1 32 IF(LRE(LI).LE.0) GO TO 30 SIT=PI3(LI) TFI=SFI(LI)/CFI(LI) TDL=SDL(LI)/CDL(LI) ALF=3.14159D0*0.25D0+DATAN(TFI)*0.5D0 BET=(ALF+SIT)*(-1.D0) IF(LRE(LI).GE.2) BET=ALF-SIT PI4(LI)=BET C E=EE(LI) P=PP(LI) G=E/(2.D0*(1.D0+P)) SI=1.D0 IF(LRE(LI).GE.2) SI=-1.D0 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./(C1*TFI*TDL+G) DPF(1,1)=C1-C2*C2*TFI*TDL*B1 DPF(1,2)=C2-C1*C2*TFI*TDL*B1 DPF(1,3)=-SI*C2*G*TDL*B1 DPF(2,1)=DPF(1,2) DPF(2,2)=C1-C1*C1*TFI*TDL*B1 DPF(2,3)=-SI*C1*G*TDL*B1 DPF(3,1)=-SI*C2*G*TFI*B1 DPF(3,2)=-SI*C1*G*TFI*B1 DPF(3,3)=G-G*G*B1 CB=DCOS(BET) SB=DSIN(BET) T(1,1)=CB*CB T(1,2)=SB*SB T(1,3)=SB*CB T(2,1)=T(1,2) T(2,2)=T(1,1) T(2,3)=-SB*CB T(3,1)=-2.D0*SB*CB T(3,2)=2.D0*SB*CB T(3,3)=CB*CB-SB*SB DO 34 I=1,3 34 FD(I)=EPP(IL,I)/DFLOAT(NCY1) CALL VULT(T,FD,DP,3,3) CALL VULT(DPF,DP,FD,3,3) T(1,1)=CB*CB T(1,2)=SB*SB T(1,3)=2.D0*SB*CB T(2,1)=T(1,2) T(2,2)=T(1,1) T(2,3)=-2.D0*SB*CB T(3,1)=-SB*CB T(3,2)=SB*CB T(3,3)=CB*CB-SB*SB CALL DINV(T,3,3,0,DET,IND) CALL VULT(T,FD,DP,3,3) DO 36 I=1,3 36 STC(LI,I)=STC(LI,I)+DP(I) GO TO 22 C 30 SX=STC(LI,1) SY=STC(LI,2) TA=STC(LI,3) CALL MOHRCO(LI,SX,SY,TA,F,S,P2,B0) B2=B0**(-0.5D0) SD=SDL(LI) PSA(1)=B2*P2-SD PSA(2)=B2*P2*(-1.D0)-SD PSA(3)=B2*4.D0*TA SF=SFI(LI) FSA(1)=B2*P2-SF FSA(2)=B2*P2*(-1.D0)-SF FSA(3)=B2*4.D0*TA GO TO 40 C 26 PSA(1)=0.D0 PSA(2)=-SDL(LI)/CDL(LI) PSA(3)=1.D0 IF(STR(LI,3).LT.0.D0) PSA(3)=-1.D0 FSA(1)=0.D0 FSA(2)=-SFI(LI)/CFI(LI) FSA(3)=1.D0 IF(STR(LI,3).LT.0.) FSA(3)=-1.D0 40 DO 42 I=1,3 DO 42 J=1,3 42 D(I,J)=DMX(LI,I,J) CALL VULT(D,PSA,DP,3,3) DO 44 I=1,3 DO 44 J=1,3 44 DPF(I,J)=DP(I)*FSA(J) CALL MULT(DPF,D,DPFD,3,3,3) DO 46 I=1,3 FD(I)=0.D0 DO 46 J=1,3 46 FD(I)=FD(I)+FSA(J)*D(J,I) FDP=0.D0 DO 48 I=1,3 48 FDP=FDP+FD(I)*PSA(I) R=1.D0/FDP DO 50 I=1,3 DO 50 J=1,3 50 DPFD(I,J)=DPFD(I,J)*R DO 52 I=1,3 DO 52 J=1,3 52 DEP(I,J)=D(I,J)-DPFD(I,J) DO 56 I=1,3 FD(I)=0.D0 DO 56 J=1,3 56 FD(I)=FD(I)+DEP(I,J)*EPP(IL,J)/DFLOAT(NCY1) DO 58 I=1,3 58 STC(LI,I)=STC(LI,I)+FD(I) 22 CONTINUE DO 60 I=1,3 60 STB(LI,I)=STC(LI,I) 20 CONTINUE RETURN END C * * * * * SUBROUTINE NOTENS C *** FIND NO-TENSION ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) COMMON/LB20/RKH(85),EM0(1680),FLO(60),SDL(1690),CDL(1690) C IT=0 DO 10 L=1,NNE JTS(L)=0 KOL2=K2(L) GO TO (10,10,11,10,14,16),KOL2 11 SG=STB(L,1) IF(SG.LE.0.D0) GO TO 10 STB(L,1)=0.D0 IT=IT+1 ITS(IT)=L JTS(L)=IT GO TO 10 14 SX=STB(L,1) SY=STB(L,2) TA=STB(L,3) CALL PRINCE(SX,SY,TA,S1,S3,TH) IF(S3.GE.0.001D0) GO TO 10 IT=IT+1 ITS(IT)=L JTS(L)=IT C=1.D0 IF(TA.LT.0.) C=-1.D0 IF(S1-0.001) 20,20,22 20 STB(L,1)=0.001D0 STB(L,2)=0.001D0 STB(L,3)=0.D0 GO TO 10 22 S=SY-SX IF(S) 30,30,32 30 IF(SX) 50,50,40 40 SY=S1-SX STB(L,2)=SY GO TO 35 32 IF(SY) 52,52,42 42 SX=S1-SY STB(L,1)=SX 35 A=(S1**2-(SX-SY)**2)*0.25D0 IF(A.LT.0.000001D0) A=0.000001D0 STB(L,3)=C*DSQRT(A) GO TO 10 50 C1=1.D0 GO TO 55 52 C1=-1.D0 55 R=(S1-S3)*0.5D0 B2=R*R-TA*TA IF(B2) 56,56,57 56 T=0.5D0*3.14159D0 GO TO 58 57 B=DSQRT(B2) AT=TA/B T=DATAN(AT) 58 S=S1*0.5D0*DCOS(T) STB(L,1)=S1*0.5D0+C1*S STB(L,2)=S1*0.5D0-C1*S STB(L,3)=S1*0.5D0*DSIN(T)*C GO TO 10 16 SG=STB(L,2) IF(SG.GE.0.D0) GO TO 10 IT=IT+1 ITS(IT)=L JTS(L)=IT STB(L,2)=0.001D0 TA=STB(L,3) ATA=DABS(TA) CF=CC(L) IF(ATA.LE.CF) GO TO 10 C=1.D0 IF(TA.LT.0.D0) C=-1.D0 STB(L,3)=C*CF 10 CONTINUE NTS=IT RETURN END C * * * * * SUBROUTINE TRUSS(L) C *** TRUSS ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) 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 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(L,I,J)=TET(I,J) AES(L)=A*EL DO 24 J=1,4 BMX(L,1,J)=0.D0 DO 24 K=1,4 24 BMX(L,1,J)=BMX(L,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(L,I)=LY(IT) RETURN END C * * * * * SUBROUTINE BEAM(L) C *** BEAM (RAHMEN) ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) 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.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(L,I,J)=EKT(I,J) EKK(L,I,J)=TET(I,J) 28 CONTINUE 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(L,I)=LY(IT) RETURN END C * * * * * SUBROUTINE PLANES(L,KOL2) C *** PLANE-STRAIN OR INTERFACE ELEMENT IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) 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(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,14,14,16),KOL2 14 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 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 IE=2*I IO=IE-1 XXL(I)=COD(IO) YYL(I)=COD(IE) 26 CONTINUE C 18 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.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 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=DABS(DA) AES(L)=AE GO TO (80,80,80,64,64,66),KOL2 64 C1=0.25 IF(MM.EQ.0) C1=1.D0 DO 70 I=1,NV DO 70 J=1,NV 70 EKK(L,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 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(L,I,J)=ES(I,J)*AES(L)*0.25D0 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 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) 80 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/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) P1=SX+SY P2=SX-SY A1=P1*SFI(L)+2.D0*CC(L)*CFI(L) 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(L,SG,TA,F,S,IC) C *** COULOMB YIELD CONDITION IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB9/ROU(1690),CC(1690),SFI(1690),CFI(1690) ST=CC(L)+SG*SFI(L)/CFI(L) 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 DINV(AA,N0,N1,N2,DET,IND) C *** INVERSE MATRIX IMPLICIT REAL*8(A-H,O-Z) DIMENSION AA(N0,N0),IPERM(3110),X(3110) 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 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 * * * * * 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.D0) VFUNC=C 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 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 *** OUTPUT FILES & MICRO-AVS DATA IMPLICIT REAL*8(A-H,O-Z) COMMON/LAB1/NNP,NNE,NFL,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB2/ISB,IDL,JBA,ICR,NRE,ISTP,MSTP,NDI,FIR COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB4/STR(1690,6),EPS(1690,6),PS0(1690,6),CC0(1690) COMMON/LAB5/DMX(1690,3,3),DIX(1690,3,3) COMMON/LAB6/EKK(1690,8,8),DBM(1690,6,8),BMX(1690,3,8) COMMON/LAB7/GKK(3110,3110),TF(3110),TU(3110),TL(3110) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB13/ST0(1690,3),S0P(1690,3),STA(960,3),STAX(1690,3) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB18/DWE(20),ERM,ERG,ERC,SCD DIMENSION UX(1680),UY(1680),UM(1680),UZ(1680),WST(1690,6),EF(8), * TT(1690) DIMENSION IYIEL(1690),INOTE(1690),IP(960) 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 UM(I)=QTU(JM)*SCU UZ(I)=0.D0 4 CONTINUE 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*DFLOAT(ISTP)/DFLOAT(NSTP) WRITE(8,220) ISTP,R,UY(NDI) WRITE(IR,202) C *** MICRO-AVS DO 6 I=1,NNP WRITE(9,221) I,UX(I),UY(I),UZ(I) 6 CONTINUE C *** STRESSES DO 10 L=1,NNE KOL2=K2(L) A=1.D0 IF(KOL2.EQ.1.OR.KOL2.EQ.3) A=PP(L) DO 12 I=1,3 WST(L,I)=STB(L,I)*A 12 IF(KOL2.EQ.3) TT(L)=-WST(L,1) DO 14 I=4,6 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.14159D0 10 CONTINUE DO 16 I=1,6 WRITE(IR,203) (WST(L,I),L=1,NNE) 16 WRITE(IR,213) WRITE(IR,204) (PI2(L),L=1,NNE) C *** YIELD ELEMENTS IF(NFAT.LE.0) GO TO 82 DO 20 K=1,NFAT L=KFA(K) UX(L)=0. KOL2=K2(L) IF(KOL2.LE.3.OR.KOL2.GE.6) GO TO 20 DO 22 I=1,3 22 WST(K,I)=EPS(L,I)*100.D0 SX=WST(K,1) SY=WST(K,2) TA=WST(K,3)*0.5D0 CALL PRINCE(SX,SY,TA,S1,S3,TH) WST(K,4)=S1 WST(K,5)=S3 WST(K,6)=TH*180.D0/3.14159D0 20 CONTINUE DO 24 I=1,6 WRITE(IR,213) 24 WRITE(IR,203) (WST(K,I),K=1,NFAT) 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 WRITE(IR,210) (IP(K),K=1,NFAT) 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 WST(K,I)=ST0(L,I) DO 34 I=1,3 34 WRITE(IR,203) (WST(K,I),K=1,NFAT) WRITE(IR,208) DO 36 K=1,NFAT L=KFA(K) DO 36 I=1,3 36 WST(K,I)=STAX(L,I) DO 38 K=1,NFAT DO 40 I=4,6 40 WST(K,I)=0.D0 L=KFA(K) IF(K2(L).LE.3.OR.K2(L).GE.6) GO TO 38 SX=WST(K,1) SY=WST(K,2) TA=WST(K,3) CALL PRINCE(SX,SY,TA,S1,S3,TH) WST(K,4)=S1 WST(K,5)=S3 WST(K,6)=TH*180.D0/3.14159D0 38 CONTINUE DO 42 I=1,6 42 WRITE(IR,201) (WST(K,I),K=1,NFAT) WRITE(IR,209) DO 44 K=1,NFAT DO 44 I=1,3 44 WST(K,I)=EPP(K,I)*100.D0 DO 46 I=1,3 46 WRITE(IR,203) (WST(K,I),K=1,NFAT) 84 WRITE(IR,218) DO 50 K=1,NFAT LI=KFA(K) IP(K)=LRE(LI) UX(K)=PI3(LI)*180.D0/3.14159D0 50 UY(K)=PI4(LI)*180.D0/3.14159D0 WRITE(IR,203) (UX(K),K=1,NFAT) WRITE(IR,203) (UY(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 100 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=SQRT(AP) WRITE(IR,211) AP,APN,APS DO 62 I=1,2 62 WRITE(IR,203) (WST(K,I),K=1,NWE) 100 CONTINUE C *** MICRO-AVS IF(IC.LE.1) GO TO 90 WRITE(9,222) WRITE(9,223) WRITE(9,224) WRITE(9,225) WRITE(9,226) WRITE(9,227) WRITE(9,228) WRITE(9,229) WRITE(9,230) DO 70 L=1,NNE 70 WRITE(9,231) L,(WST(L,I),I=1,5),TT(L),IYIEL(L),INOTE(L) C 200 FORMAT(/' LOADING STEP',I3/'DISPLACEMENTS X, Y') 201 FORMAT(10F8.4) 202 FORMAT('STRESSES SX,SY,TA,S1,S3,ANGLE') 203 FORMAT(10F8.3) 204 FORMAT('SAFETY-FACTOR'/(10F7.3)) 205 FORMAT('YIELD-ELEMENTS'/(10I5)) 206 FORMAT('YIELD-FACTOR'/(10F7.3)) 207 FORMAT('INITIAL-STRESSES') 208 FORMAT('YIELD-STRESSES SX,SY,TA,S1,S3,ANGLE') 209 FORMAT('PLASTIC STRAINS %') 210 FORMAT('FAILURE STEP'/(10I5)) 211 FORMAT('ERATH-THRUST'/' RESULTANT NORMAL SHEAR'/3F10.3) 212 FORMAT('NO-TENSION ELEMENTS'/(10I5)) 213 FORMAT(' ') 214 FORMAT('SHEAR STRAIN %'/(10F7.2)) 215 FORMAT(7E10.3) 216 FORMAT('EQUIVALENT NODAL FORCES'/' EL NODES') 217 FORMAT(16I4) 218 FORMAT('ANGLES OF S1 & SLIP-SURFACE') 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) 90 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,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) DIMENSION SX(2100),SY(2100),XXX(1300),YYY(1300) 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,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB11/LY(5400),NDF(3400),NBA(85),LBA(85,990),LOI(60) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) COMMON/LB18/DWE(20),ERM,ERG,ERC,SCD DIMENSION UX(1680),UY(1680),UM(1680),WST(1690,6),PI9(1690) C QTU(NN)=0. 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,NSTP,NWE,NIS,NFA,NFAT,NY,NN,NTR,ITER,NCY COMMON/LAB3/XX(1680),YY(1680),IX(1680),IY(1680),IQ(1680),ZZ(1680) COMMON/LAB8/AES(1690),EE(1690),PP(1690),HH(1690) COMMON/LB10/IJK(1690,4),K2(1690),LLL(1690,8),LRE(1690),IRE(50) COMMON/LB12/PI1(1690),PI2(1690),PI3(1690),PI4(1690),PI5(1690) COMMON/LB13/ST0(1690,3),S0P(1690,3),STA(960,3),STAX(1690,3) COMMON/LB14/STB(1690,3),STC(1690,3),EPP(960,3),PW(1690,2) COMMON/LB15/IFA(960),JFA(1690),KFA(960),LFA(1690),MFA(1690) COMMON/LB16/NTS,ITS(960),JTS(1690),IFS(1690),IWE(30) COMMON/LB17/QST(1690,6),PST(1690,6),QTU(3110),PTU(3110) DIMENSION WST(6),AX(2100),AY(2100),XXX(1300),YYY(1300) 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 DO 12 I=1,3 12 WST(I)=STB(L,I) 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