random.for

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

FROGRAM  RANDOM
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
        STOCHASTIC STRUCTURE NUMERICAL ANALYSIS
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          DIMENSION D1(3,3), DA(3,3,400), UV(6), 
      ﹡  SGI(3), UCI(4,600), SCI(2,3,400), CI(4,400)
          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
          COMMON /TRCOV/COV(150,150), DIB3(150), 
      ﹡  PSI(150,150), DIA3(150)
      ﹡  /TRCOV2/P(4), H(4), SP(4), DSP(2,4)
      ﹡  /TRCOV3/INEL(40,2)
      ﹡  /SMCOV/SMEAN(3,400), COVAR(3,3,400)
      ﹡  /SKC/SK(12000), TK(12000)
      ﹡  /NSAM/ SA(150)
      ﹡  /SUIJI5/DE(3,3), DALTA(3,3)
      ﹡  /RAN/XU, RI
          DATA  P/-0.86111363, -0.33988, 0,33988, 0.86111363/
          DATA  P/0.3478548, 0.6521452, 0,6521452, 0.3478548/
          WRITE(﹡,2000)
2000      FORMAT(1X,'START')
          OPEN (15, FILE='IPT3.DAT', STATUS='OLD')
          OPEN (16, FILE='OPT3.DAT', STATUS='NEW')
          CALL READD
          WRITE(﹡,2001)
2001      FORMAT(1X,'READD END')
          NN=NX*2
          CALL FMDA
          WRITE(﹡,2002)
2002      FORMAT(1X,'FMDA END')
          CALL DIAVAR
          WRITE(﹡,2004)
2003      FORMAT(1X,'DIAVAR END')
          WRITE(﹡,3003)
3003      FORMAT(1X,'EIGNPROBLEM START')
          CALL TRED2(NEA+NEL,NEA+NEL)
          CALL IMTQL2(NEA+NEL,NEA+NEL,IERR)
          WRITE(﹡,4004)
4004      FORMAT(1X,'EIGNPROBLEM END')
          DO 110 I=1,3
          DO 110 J-1,3
          DO 120 K=1,2
          D(I,J,K)=0.0
120       DD2EG(I,J,K)=0.0
          DO 110 K=1,4
          DD(I,J,K)=0.0
110       DD2(I,J,K)=0.0
          CALL  TDEG
          WRITE(﹡,2003)
2004      FORMAT(1X,'ELASTIC MATRIX END')
          IF(NERL.NE.1)  GOTO 130
          MX=NEL+NEA
          XY=FLOAT(MX)/2.
          MQ=IFIX(XY)
          DO 131 I=1,MQ
          K=MX-I+1
          AX=DIA3(I)
          DIA3(I)=DIA3(K)
          DIA3(K)=AX
          DO 131 J=1,MX
          AX=PSI(J,I)
          PSI(J,I)=PSI(J,K)
          PSI(J,K)=AX
131       CONTINUE
          CALL ERLE
          GOTO 8888
130       NEB=NEA+NEL
8888      WRITE(16,2005)
2005      FORMAT(/1X,27 ('﹡'), 'OPERATION RESULTS',27('﹡')/)
          DO 90 I=1,MN
90        SK(I)=0.0
          DO 30 IE=1,NE
          LT=INE(IE,4)
          DO 35 I=1,3
          DO 35 J=1,3
35        D1(I,J)=D(I,J,LT)
          CALL BSK(IE)
          CALL ESK(IE,D1)
          CALL FORMK(IE)
30        CONTINUE
          WRITE(﹡,4006)
4006      FORMAT(1X,'ASSEMBLE END')
          CALL FORF
          CALL BOOND
          CALL SLOVE
          CALL BACK
          DO 100 J=1,NN
100       VB(J)=F(J)
          WRITE(16,9097)
9097      FORMAT(1X,22('﹡','THE DETERMINATE   
       ﹡ DISPLACEMENT',22('﹡')/)
          WRITE(16,7001) (VB(J),J=1,NN)
7001      FORMAT(6F12.7)
          WRITE(﹡,4007)
4007      FORMAT(1X,'VB( ) END')
          NYC=0
          CALL DDPJE
          WRITE(﹡,4008)
4008      FORMAT(1X,'VC(,) END')
          IF(LNY.NE.2.AND.LNY.NE.3) GOTO 4010
          CALL D2DPJE
          WRITE(﹡,4009)
4009      FORMAT(1X,'V2( ) END')
4010      CONTINUE
          WRITE(﹡,3009)
3009      FORMAT(1X,'ANALYSIS START')
          CALL DISPL
          CALL STOCH
          CALL SETOUT
          STOP '...END...'
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
     INPUT THE BASIC INFORMATION ABOUT THE CALCULATION 
                  OF FINITE ELEMENT
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE READD
          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)
          WRITE(16,9001)
9001      FORMAT(1X, '...PLEASE CHECK THE INPUTED DATA:'//)
          WRITE(16,9002)
9002      FORMAT(1X,20('﹡'),'NX,NE,MT,NR,MP,
       ﹡ MTX,LNY,NERL,NYC',20('﹡')/)
          READ(15,110,END=11)NX,NE,MT,NR,MP,MTX,LNY,NERL,NYC
110       FORMAT(9I3)
11        WRITE(16,120)NX,NE,MT,NR,MP,MTX,LNY,NERL,NYC
120       FORMAT(13X,9I5)
          WRITE(16,40)
40        FORMAT(/1X,33('﹡'),'MY(,)',33('﹡')/)
          READ(15,20,END=25) ((MY(I,J),J=1,2),I=1,2)
20        FORMAT(4I3)
25        WRITE(16,30) ((MY(I,J),J=1,2),I=1,2)
30        FORMAT(26X,4I5)
          WRITE(16,9003)
9003      FORMAT(/1X,33('﹡'),'AE(,)',33('﹡')/)
          READ(15,130,END=13) ((AE(N,I),I=1,8),N=1,MT)
13        WRITE(16,131) ((AE(N,I),I=1,8),N=1,MT)
          WRITE(16,9004)
9004      FORMAT(/1X,33('﹡'),'X-Y',33('﹡')/)
          READ(15,140,END=14) (X(I),I=1,NX),(Y(I),I=1,NX)
14        WRITE(16,141) (X(I),I=1,NX),(Y(I),I=1,NX)
          WRITE(16,9005)
9005      FORMAT(/1X,33('﹡'),'INE(,)',33('﹡')/)
          READ(15,150,END=15) ((INE(N,I),I=1,5),N=1,NE)
15        WRITE(16,151)  ((INE(N,I),I=1,5),N=1,NE)
          WRITE(16,9006)
9006      FORMAT(/1X,33('﹡'),'RR( )',33('﹡')/) 
          READ(15,160,END=16) (RR(I),I=1,NR)
16        WRITE(16,161)  (RR(I),I=1,NR)
          WRITE(16,9007)
9007      FORMAT(/1X,33('﹡'),'NPP( )',33('﹡')/)
          READ(15,170,END=17) (NPP(I),I=1,MP)
17        WRITE(16,171) (NPP(I),I=1,MP)
          MP2=2*MP
          WRITE(16,9008)
9008      FORMAT(/1X,33('﹡'),'PP( )',33('﹡')/)
          READ(15,180,END=18) (PP(I),I=1,MP2)
18        WRITE(16,181) (PP(I),I=1,MP2)
130       FORMAT(8E8.5)
140       FORMAT(6F5.2)
150       FORMAT(5I3)
160       FORMAT(2F10.8)
170       FORMAT(2I3)
180       FORMAT(4F7.2)
131       FORMAT(1X,8E9.3)
141       FORMAT(12F6.2)
151       FORMAT(6X,10I6)
161       FORMAT(6F12.8)
171       FORMAT(11X,10I5)
181       FORMAT(10F7.2)
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
                CALCULATE THE ELASTIC MATRIX
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE TDEG
          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)
          /SUIJI2/D(3,3,2),DD(3,3,4),DD2(3,3,4),
       ﹡ DD2EG(3,3,2),EK(6,6)
          DO 40 LT=1,MT
          E=AE(LT,1)
          G=AE(LT,2)
          EGT=4.*G-E
          EGS=E-2.*G
          D(1,1,LT)=4.*G*G/EGT
          D(2,2,LT)=D(1,1,LT)
          D(1,2,LT)=2.*G*EGS/EGT
          D(2,1,LT)=D(1,2,LT)
          D(3,3,LT)=G
          MR=(LT-1)*2+1
          MS=(LT-1)*2+2
          DO 50 I=1,2
          DO 50 J=1,2
50        DD(I,J,MR)=4.*G*G*E/EGT/EGT*MY(LT,1)
          DD(1,1,MS)=-8.*G*G*EGS/EGT/EGT*MY(LT,2)
          DD(2,2,MS)=DD(1,1,MS)
          DD(1,2,MS)=-2.*G*EGS*EGS/EGT/EGT*MY(LT,2)
          DD(2,1,MS)=DD(1,2,MS)
          DD(3,3,MS)=G*MT(LT,2)
          DO 60 I=1,2
          DO 60 I=1,2
          DD2(I,J,MR)=8.*G*G*E*E/EGT/EGT*MY(LT,1)
          DD2(I,J,MR)=8.*G*G*E*E/EGT/EGT*MY(LT,2)
60        DD2EG(I,J,LT)=-8.*G*G*E*E/EGT/EGT*
       ﹡ MY(LT,1)*MY(LT,2)
40        CONTINUE
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
                 CALCULATE THE STRAIN MATRIX
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE BSK(E)
          DIMENSION B(3,6)
          COMMON/FONDA/AE(2,8),X(300),Y(300),
       ﹡ INE(400,5),RR(40),NPP(30),PP(60)
       ﹡ /SUIJI3/BB(3,6,400),AREA(400)
          II=INE(IE,1)
          JJ=INE(IE,2)
          KK=INE(IE,3)
          CJ=X(II)-X(KK)
          CK=X(JJ)-X(II)
          BJ=Y(KK)-Y(II)
          BK=Y(II)-Y(JJ)
          AREA(IE)=O.5*(BJ*CK-BK*CJ)
          IF(AREA(IE).LE.0.0)  GOTO 450
          B(1,1)=-0.5*(BJ+BK)
          B(1,2)=0.0
          B(1,3)=O.5*BJ
          B(1,4)=0.0
          B(1,5)=0.5*BK
          B(1,6)=0.0
          B(2,1)=0.0
          B(2,2)=-0.5*(CJ+CK)
          B(2,3)=0.0
          B(2,4)=0.5*CJ
          B(2,5)=0.0
          B(2,6)=0.5*CK
          B(3,1)=B(2,2)
          B(3,2)=-0.5*(BJ+BK)
          B(3,3)=B(2,4)
          B(3,4)=B(1,3)
          B(3,5)=B(2,6)
          B(3,6)=B(1,5)
          DO 900 M=1,3
          DO 900 N=1,6
900       BB(M,N,IE)=B(M,N)
          RETURN
450       WRITE(﹡,460)
460       FORMAT(5X,'THE AREA OF FINITE ELEMENT',
       ﹡ I3,'IS ZERO OR NEGATIVE')
          STOP
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
         CALCULATE THE ELEMENT RIGID MATRIX
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE ESK(IE,S)
          DIMENSION  S(3,3)
          COMMON/SUIJI3/BB(3,6,400),AREA(400)
       ﹡ /SUIJI2/D(3,3,2),DD(3,3,4),DD2(3,3,4)
       ﹡ DD2EG(3,3,2),EK(6,6)
          DO 10 I=1,6
          DO 10 J=1,6
          EK(I,J)=0.0
          DO 10 K=1,3
          DO 10 L=1,3
10        EK(I,J)=EK(I,J)+BB(K,I,IE)*S(K,L)*
       ﹡ BB(L,J,IE)/AREA(IE)
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          FORM THE TOTAL RIGID MATRIX
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE FORMK(IE)
          DIMENSION IN(6)
          COMMON/FONDA/AE(2,8),X(300),Y(300),
       ﹡ INE(400,5),RR(40),NPP(30),PP(60)
       ﹡ /SUIJI2/D(3,3,2),DD(3,3,4),DD2(3,3,4),
       ﹡ DD2EG(3,3,2),EK(6,6)
       ﹡ /SUIJI4/MDA(600),MXB,MN
       ﹡ /SKC/SK(12000),TK(12000)
          IN(1)=INE(IE,1)*2-1
          IN(2)=IN(1)+1
          IN(3)=INE(IE,2)*2-1
          IN(4)=IN(3)+1
          IN(5)=INE(IE,3)*2-1
          IN(6)=IN(5)+1
          DO 620 I=1,6
          I1=IN(I)
          DO 620 J=1,6
          J1=IN(J)
          IF(I1-J1)620,610,610
610       K=MDA(I1)-I1+J1
          SK(K)=SK(K)+EK(I,J)
620       CONTINUE
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          THE CALCULATION OF THE MAXIMUM BIND WIDETH
                AND THE ELEMENT AMOUNT IN SK
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE FMDA
          DIMENSION MF(3)
          COMMON/CONTR/NX,NE,MT,NR,MP,MTX,NN,LNY,NERL,NYC
          COMMON/FONDA/AE(2,8),X(300),Y(300),
       ﹡ INE(400,5),RR(40),NPP(30),PP(60)
       ﹡ /SUIJI4/MDA(600),MXB,MN
          DO 510 I=1,100
510       MDA(I)=0
          MDA(1)=1
          MDA(2)=2
          DO 540 I=1,NE
          II=INE(I,1)
          JJ=INE(I,2)
          KK=INE(I,3)
          MF(1)=II
          MF(2)=JJ
          MF(3)=kk
          DO 540 J=1,3
          DO 540 K=1,3
          IF(MF(J)-MF(K)) 540,540,520
520       I1=(MF(J)-MF(K))*2
          L1=MF(K)*2
          IF(MDA(L1)-(I1+2)) 530,540,540
530       MDA(L1)=I1+2
          MDA(L1-1)=I1+1
540       CONTINUE
          MXB=MDA(1)
          DO 550 I=2,NN
          IF(MDA(I).GT.MXB) MXB=MDA(I)
550       MDA(I)=MDA(I)+MDA(I-1)
          MN=MDA(NN)
          WRITE(16,570)
570       FORMAT(/1X,34('﹡'),'MDA',34('﹡')/)
          WRITE(16,580)  (I,MDA(I),I=1,NN)
580       FORMAT((1X,6('(',I3,'),I5,2X)))
          WRITE(16,590) MXB
590       FORMAT(/1X,'THE MAXINUM BIND WIDETH=',I5/)
          WRITE(16,600) MN
600       FORMAT(/1X,'THE ELEMENT AMOUNT IN SK=',I5/)
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
                     QUOTE BOUNDARY CONDITION
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE BOOND
          COMMON/CONTR/NX,NE,MT,NR,MP,MTX,NN,LNY,NERL,NYC
      ﹡  /SKC/SK(12000), TK(12000)
      ﹡  /FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
      ﹡  /SUIJI4/MDA(600),MXB,MN
          DO 540 I=1,NR
          A=RR(I)*1000.0
          I1=IFIX(A)
          A=(A-FLOAT(I1)*10.0
          K=IFIX(A)
          A=(A-FLOAT(K))*10.0
          L=IFIX(A)
          IF(K) 510,510,520
510       I2=MDA(2*I1-1)
          SK(I2)=1.7E20
520       IF(L) 530,530,540
530       I3=MDA(2*I1)
          SK(I3)=1.7E20
540       CONTINUE
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
                CALCULATE LOAD MATRIX
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE FORF
          COMMON/CONTR/NX,NE,MT,NR,MP,MTX,NN,LNY,NERL,NYC
       ﹡ /SUIJI1/VB(600), VC(600,150), V2(600), F(600), QF(600)
       ﹡ /FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
       ﹡ RR(40), NPP(30), PP(60)
          DO 610 I=1,NN
610       F(I)=0.0
          DO 630 J=1,MP
          L=NPP(J)*2
          IF(L) 630,630,620
620       K=L-1