random.for

随机有限元分析的经典程序

          F(L)=PP(2*J)
          F(K)=PP(2*J-1)
630       CONTINUE
          DO 640 I=1,NN
640       QF(I)=F(I)
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
             RESOLVE THE TORAL RIGID MATRIX
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE SLOVE 
          COMMON/CONTR/NX,NE,MT,NR,MP,MTX,NN,LNY,NERL,NYC
      ﹡  /SUIJI4/MDA(600),MXB,MN
      ﹡  /SKC/SK(12000), TK(12000)
          DO 100 J=1,NN
          K2=MDA(J)
          B=0.0
          IF(J.EQ.1) GOTO 20
          M3=MDA(J-1)
          DO 10 K=1,J-1
          M1=MDA(J)-J+K
          IF(M1.LE.M3) GOTO 10
          B=B+SK(M1)*SK(M1)
10        CONTINUE
20        SK(K2)=SQRT(SK(K2)-3)
          DO 80 I=J+1,NN
          K1=MDA(I)-I+J
          A=0.0
          K3=MDA(I-1)
          IF(K1.LE.K3) GOTO 80
          IF(J.EQ.1) GOTO 5
          DO 30 K=1,J-1
          M1=MDA(I)-I+K
          M2=MDA(I)-J+K
          IF(M1.LE.K3.OR.M2.LE.M3) GOTO 30
          A=A+SK(M1)*SK(M2)
30        CONTINUE
5         SK(K1)=(SK(K1)-A)/SK(K2)
80        CONTINUE
100       CONTINUE
          DO 110 I=1,MN
110       TK(I)=SK(I)
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
       CALCULATE THE AVERAGE VALUE OF VELOCITY
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE BACK
          COMMON/CONTR/NX,NE,MT,NR,MP,MTX,NN,LNY,NERL,NYC
      ﹡  /SUIJI1/VB(600), VC(600,150), V2(600), F(600), QF(600)
      ﹡  /SUIJI4/MDA(600),MXB,MN
      ﹡  /SKC/SK(12000), TK(12000)
          F(1)=F(1)/TK(1)
          DO 20 I=2,NN
          K=I-1
          K1=MDA(I)
          K2=MDA(K)
          A=0.0
          DO 10 J=1,K
          M=MDA(I)-I+J
          IF(M.LE.K2) GOTO 10
          A=A+TK(M)*F(J)
10        CONTINUE
          F(I)=(F(I)-A)/TK(K1)
20        CONTINUE
          F(NN)=F(NN)/TK(MN)
          L=0
          DO 40 I=1,NN-1
          L=L+1
          L1=NN-L
          A=0.0
          DO 30 J=L1+1,NN
          K=MDA(J)-J+L1
          IF(K.LE.MDA(J-1)) GOTO 30
          A=A+TK(K)*F(J)
30        CONTINUE
          M=MDA(L1)
          F(L1)=(F(L1)-A)/TK(M)
40        CONTINUE
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡ 
           CALCULATE TEH SENSITIVITY DEVIATION]
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE DDPJE
          DIMENSION D1(3,3)
          COMMON/CONTR/NX,NE,MT,NR,MP,MTX,NN,LNY,NERL,NYC
      ﹡  /CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2), NEL
      ﹡  /FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
      ﹡  /SUIJI1/VB(600), VC(600,150), V2(600), F(600), QF(600)
      ﹡  /SUIJI2/D(3,3,2), DD(3,3,4), DD2(3,3,4), 
      ﹡  DD2EG(3,3,2), EK(6,6)
      ﹡  /SUIJI3/BB(3,6,400), AREA(400)
      ﹡  /SUIJI4/MDA(600), MXB, MN
      ﹡  /TRCOV/COV(150,150), DIB3(150), 
      ﹡  PSI(150,150), DIA3(150)
      ﹡  /TRCOV3/INEL(40,2)
      ﹡  /SKC/SK(12000), TK(12000)
          DO 40 JE=1,NEB
          WRITE(﹡,1) JE
1         FORMAT(1X,'JE=',I4)
          DO 5 I=1,MN
5         SK(I)=0.0
          DO 10 I=1,NN
10        F(I)=0.0
          IF(NEL.GT.1) GOTO 60
          DO 50 I=1,NN
50        F(I)=QF(I)*PSI(NEA+1,JE)
          GOTO 90
60        DO 80 KE=NEA+1,NEA+NEL 
          K=KE-NEA
          DO 80 I=1,2
          DO 80 I=1,2
          M=2*(INEL(K,I)-1)+J
80        F(M)=F(M)+QF(M)*PSI(KE,JE)*0.5
          IF(LNY.NE.2.AND.LNY.NE.3) GOTO 26
90        DO 25 KE=1,NEA
          DO 20 IE=1,NE
          CALL IEKE(IE,KE,KR,MR)
          IF(KR.EQ.0) GOTO 20
          DO 12 I=1,3
          DO 12 J=1,3
12        D1(I,J)=DD(I,J,MR)
          CALL ESK(IE,D1)
          DO 15 I=1,6
          DO 15 J=1,6
15        EK(I,J)=EK(I,J)*PSI(KE,JE)
          CALL FKRMK(IE)
20        CONTINUE
25        CONTINUE
          CALL ART(VB)
26        CONTINUE
          CALL BACK
          DO 30 I=1,NN
30        VC(I,JE)=F(I)
40        CONTINUE
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          CALCULATE THE SECOND ORDER DEVIATION
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE D2DPJE
          DIMENSION V(600),D1(3,3)
          COMMON /CONTR/NX, NE, MT, NR, MP, MTX,
      ﹡  NN, LNY, NERL, NYC
      ﹡  /CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2), NEL
      ﹡  /FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
      ﹡  /SUIJI1/VB(600), VC(600,150), V2(600), F(600), QF(600)
      ﹡  /SUIJI2/D(3,3,2), DD(3,3,4), DD2(3,3,4), 
      ﹡  DD2EG(3,3,2), EK(6,6)
      ﹡  /SUIJI3/BB(3,6,400), AREA(400)
      ﹡  /SUIJI4/MDA(600), MXB, MN
      ﹡  /TRCOV/COV(150,150), DIB3(150), 
      ﹡  PSI(150,150), DIA3(150)
      ﹡  /SKC/SK(12000), TK(12000)
          DO 10 I=1,NN
10        F(I)=0.0
          DO 70 JE=1,NEB
          DO 5 J=1,MN
5         SK(J)=0.0
          DO 30 KE=1,NEA
          DO 20 IE=1,NE
          CALL IEKE(IE,KE,KR,MR)
          IF(KR.EQ.0) GOTO 20
          DO 12 I=1,3
          DO 12 J=1,3
12        D1(I,J)=DD(I,J,MR)
          CALL ESK(IE,D1)   
          DO 15 I=1,6
          DO 15 J=1,6
15        EK(I,J)=EK(I,J)*PSI(KE,JE)*DIA3(JE)
          CALL FORMK(IE)
20        CONTINUE
30        CONTINUE
          DO 40 I=1,NN
40        V(I)=VC(I,JE)
          CALL ART(V)
          DO 41 I=1,MN
41        SK(I)=0.0
          DO 60 KE=1,NEA
          DO 50 IE=1,NE
          CALL IEKE(IE,KE,KR,MR)
          IF(KR.EQ.0) GOTO 50
          DO 42 I=1,3
          DO 42 I=1,3
42        D1(I,J)=DD2(I,J,MR)
          CALL ESK(IE,D1)
          DO 45 I=1,6
          DO 45 J=1,6
45        EK(I,J)=EK(I,J)*PSI(KE,JE)**2*DIA3(JE)
          CALL FORMK(IE)]
50        CONTINUE
60        CONTINUE
          CALL ART(VB)
          DO 90 I=1,MN
90        SK(I)=0.0
          DO 100 KE=1,NEA
          DO 110 IE=1,NE
          INT=INE(IE,5)
          MTY=INE(IE,4)
          INT=INT+(MTY-1)*NEADA(1)
          INS=(MY(MTY,1)+MTY(MTY,2)-1)*(INT+NEAD(MTY))
          IF(INS.EQ.0) GOTO 110
          IF(KE.EQ.INT.OR.KE.EQ.INS) GOTO 120
          GOTO 110
120       CONTINUE
          DO 130 I=1,3
          DO 130 J=1,3
130       D1(I,J)=DD2EG(I,J,MTY)
          CALL ESK(IE,D1)
          DO 140 I=1,6
          DO 140 J=1,6
140       EK(I,J)=EK(I,J)*PSI(INT,JE)*PSI(INS,JE)*DIA3(JE)
          CALL FORMK(IE)
110       CONTINUE
100       CONTINUE
          CALL ART(VB)
70        CONTINUE
          CALL BACK
          DO 80 I=1,NN
80        V2(I)=F(I)
          RETURN
          END
         
          SUBROUTINE IEKE(IE,KE,KR,MR)
          COMMON/CONTR1/NEAD(2), NEADA(2), 
      ﹡  NEA, NEB, MY(2,2), NEL
      ﹡  /FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
          MTY=INE(IE,4)
          INT=(MTY-1)*NEADA(1)+INE(IE,5)
          INS=(MY(MTY,1)+MY(MTY,2)-1)*(INT+NEAD(MTY))
          IF(KE.EQ.INT) GOTO 10
          IF(KE.EQ.INS) GOTO 20
          KR=0
          RETURN
10        MR=(MTY-1)*2+1
          IF(MY(MTY,1).EQ.0) MR=MR+1
          GOTO 30
20        MR=(MTY-1)*2+2
30        KR=INT
          RETURN
          END

          SUBROUTINE ART(V)
          DIMENSION V(NN)
          COMMON /CONTR/NX, NE, MT, NR, MP, MTX,
      ﹡  NN, LNY, NERL, NYC
      ﹡  /SUIJI1/VB(600), VC(600,150), V2(600), F(600), QF(600)    
      ﹡  /SUIJI4/MDA(600), MXB, MN
      ﹡  /SKC/SK(12000), TK(12000)
          F(1)=F(1)-SK(1)*V(1)
          DO 20 I=2,NN
          DO 10 J=1,I
          M=MDA(I)-I+J
          IF(M.LE.MDA(I-1)) GOTO 10
          F(I)=F(I)-SK(M)*V(J)
10        CONTINUE
20        CONTINUE
          DO 40 I=1,NN-1
          DO 30 J=I+1,NN
          M=MDA(J)-J+1
          IF(M.LE.MDA(J-1)) GOTO 30
          F(I)=F(I)-SK(M)*V(J)
30        CONTINUE
40        CONTINUE
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
      INPUT THE INFORMATION ABOUT LOCAL AVERAGE ELEMENT AND 
        CALCULATE COVARIANT OF LOCAL AVERAGE RANDOM FIELD
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE DIAVAR
          COMMON/FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
      ﹡  /CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2), NEL
      ﹡  /TRCOV/COV(150,150), DIB3(150), 
      ﹡  PSI(150,150), DIA3(150)
      ﹡  /TRCOV1/VAR(4), SCA(10), INEA(150,5), AERA(150)
      ﹡  /TRCOV2/P(4), H(4), SP(4), DSP(2,4)
          WRITE(16,9002)
9002      FORMAT(/1X,33('﹡'),'VAR( )',33('﹡')/)
          READ(15,110,END=11) (VAR(I),I=1,4)
110       FORMAT(4F12.8)
11        WRITE(16,111) (VAR(I),I=1,4)
111       FORMAT(4F18.8)
          WRITE(16,9010)
9010      FORMAT(/1X,33('﹡'),'SCA( )',33('﹡')/)
          READ(15,112,END=12) (SCA(I),I=1,10)
112       FORMAT(5E14.4)
12        WRITE(16,113) (SCA(I),I=1,10)
113       FORMAT(1X,5E14.4)
          READ(15,120,END=13) (NEAD(I),I=1,2)
120       FORMAT(2I3)
13        WRITE(16,121) (NEAD(I),I=1,2)
121       FORMAT(/1X,'MERTIRAL 1# LOCAL AVERAGE MESH NUMBER=',
      ﹡  I3/1X,'MERTIRAL 2# LOCAL AVERAGE MESH NUMBER=',I3/)
          DO 10 I=1,2
10        NEADA(I)=(MY(I,1)+MY(I,2))*NEAD(I)
          NEA=NEADA(1)+NEADA(2)
          WRITE(16,3000)
3000      FORMAT(/1X,'LOCAL AVERAGE ELEMENTS AMOUNT=',I3/)
          NP=NEAD(1)+NEAD(2)
          WRITE(16,9003)
9003      FORMAT(/1X,32('﹡'),'INEA(,)',32('﹡')/)
          READ(15,150,END=9) ((INEA(N,1),I=1,5),N=1,NP)
9         WRITE(16,151) ((INEA(N,1),I=1,5),N=1,NP)
150       FORMAT(5I3)
151       FORMAT(11X,10I5)
          DO 2000 I=1,NP
          AERA(I)=0.0
          DO 2000 J=1,4
          N=J+1
          IF(N.GT.4) N=1
          M=INEA(I,J)
          K=INEA(I,N)
2000      AERA(I)=AERA(I)+(X(M)-X(K))*(Y(M)+Y(K))/2.0
          DO 160 M=1,NEA
          CALL DVM(M,LV1,LM)
          I1=INEA(LM,1)
          I2=INEA(LM,2)
          I3=INEA(LM,3)
          DO 160 N=1,M
          CALL DVM(N,LV2,LN)
          LX=LV1
          IF(LV2.LT.LV1) LX=LV2
          GOTO (152,153,154,154) LX
152       LY=LV1*LV2
          GOTO 156
153       LY=LV1+LV2+1
          GOTO 156
154       LY=LV1+LV2+2
156       SCB=SCA(LY)
          IF (SCB.EQ.0.0) GOTO 160
          IF(INEA(LM,5).EQ.0.AND.INEA(LN,5).EQ.0) GOTO 2
          CALL UDIAVAR(M,N,LM,LN,LV1,LV2,SCB) GOTO 160
2         J1=INEA(LN,1)
          J2=INEA(LN,2)
          J3=INRA(LN,3)
          TX1=X(I2)-X(I1)
          TY1=Y(I3)-Y(I2)
          TX2=X(J2)-X(J1)
          TY2=Y(J3)-Y(J2)
          T10=ABS(X(I2)-X(J1))
          T11=ABS(X(I1)-X(J1))
          T12=ABS(X(I1)-X(J2))
          T13=ABS(X(I2)-X(J2))
          T20=ABS(Y(I3)-Y(J2))
          T21=ABS(Y(I2)-Y(J2))
          T22=ABS(Y(I2)-Y(J3))
          T23=ABS(Y(I3)-Y(J3))
          G10=1.0-T10/SCB/3.0
          IF(T10-SCB) 172,172,171
171       G10=(1.0-SCB/T10/3.0)*SCB/T10
172       G11=1.0-T11/SCB/3.0
          IF(T11-SCB) 174,174,173
173       G11=(1.0-SCB/T11/3.0)*SCB/T11
174       G12=1.0-T12/SCB/3.0
          IF(T12-SCB) 176,176,175
175       G12=(1.0-SCB/T12/3.0)*SCB/T12
176       G13=1.0-T13/SCB/3.0
          IF(T13-SCB) 178,178,177
177       G13=(1.0-SCB/T13/3.0)*SCB/T13
178       G20=1.0-T20/SCB/3.0
          IF(T20-SCB) 180,180,179
179       G20=(1.0-SCB/T20/3.0)*SCB/T20
180       G21=1.0-T21/SCB/3.0
          IF(T21-SCB) 182,182,181
181       G21=(1.0-SCB/T21/3.0)*SCB/T21
182       G22=1.0-T22/SCB/3.0
          IF(T22-SCB) 184,184,183
183       G22=(1.0-SCB/T22/3.0)*SCB/T22
184       G23=1.0-T23/SCB/3.0
          IF(T21-SCB) 186,186,185
185       G23=(1.0-SCB/T23/3.0)*SCB/T23
186       C00=T10*T10*T20*T20*G10*G20
          C01=T10*T10*T21*T21*G10*G21
          C02=T10*T10*T22*T22*G10*G22
          C03=T10*T10*T23*T23*G10*G23
          C10=T11*T11*T20*T20*G11*G20
          C11=T11*T11*T21*T21*G11*G21
          C12=T11*T11*T22*T22*G11*G22
          C13=T11*T11*T23*T23*G11*G23
          C20=T12*T12*T20*T20*G12*G20
          C21=T12*T12*T21*T21*G12*G21
          C22=T12*T12*T22*T22*G12*G22
          C23=T12*T12*T23*T23*G12*G23
          C30=T13*T13*T20*T20*G13*G20
          C31=T13*T13*T21*T21*G13*G21
          C32=T13*T13*T22*T22*G13*G22
          C00=T13*T13*T23*T23*G13*G23