a06.for

ADINA84有限元编程学习的好例子

      SUBROUTINE ELT2D4
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1           ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /DIMEL/ N101,N102,N103,N104,N105,N106,N107,N108,N109,N110,
     1               N111,N112,N113,N114,N120,N121,N122,N123,N124,N125
      COMMON /MATMOD/ STRESS(4),STRAIN(4),D(4,4),IPT,NEL,IPS,ISVE
      COMMON /CRACK/ GAMMA,STIFAC,SHEFAC,SIGMAT,SIGMAC,EPSC,SIGMAU,EPSU,
     1               RAM5,RBM5,RCM5,ICRACK,MOD45
      COMMON /CONCR2/ BETA,GAMA,RKAPA,ALFA,SIGP(4),TEP(4),EP(4),YP(3),
     1                E,VNU,RK,G,E12,E14,E24,EPSCP,SIGCP,FALSTR,ILFSET
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /DPR/ ITWO
      COMMON A(1)
      REAL A
      DIMENSION IA(1)
C
      EQUIVALENCE (NPAR(7),MXNODS),(NPAR(10),NINT),(NPAR(17),NCON)
      EQUIVALENCE (A(1),IA(1))
C
C
      IDW=15*ITWO
      NPT=NINT*NINT
      MATP=IA(N107 + NEL - 1)
      NM=N109 + (MATP - 1)*NCON*ITWO
      NNOD=N110 + (NEL - 1)*(IDW*NPT + MXNODS) + IDW*NPT
      N6A1=N6A + ITWO
      N6B1=N6B + ITWO
C
C
      IF (NPAR(15).EQ.5) GO TO 30
      NM25=NM + 24*ITWO
      ICRACK=IA(NM25)
      GAMMA=DOUBLE (A(NM25 + ITWO))
      STIFAC=DOUBLE (A(NM25 + 2*ITWO))
      SHEFAC=DOUBLE (A(NM25 + 3*ITWO))
      GO TO 50
C
   30 SIGMAT=DOUBLE (A(NM + 3*ITWO))
      SIGMAC=DOUBLE (A(NM + 4*ITWO))
      EPSC  =DOUBLE (A(NM + 5*ITWO))
      SIGMAU=DOUBLE (A(NM + 6*ITWO))
      EPSU  =DOUBLE (A(NM + 7*ITWO))
      BETA=DOUBLE (A(NM + 32*ITWO))
      GAMA=DOUBLE (A(NM + 33*ITWO))
      RKAPA=DOUBLE (A(NM + 34*ITWO))
      ALFA=DOUBLE (A(NM + 35*ITWO))
      STIFAC=DOUBLE (A(NM + 36*ITWO))
      SHEFAC=DOUBLE (A(NM + 37*ITWO))
C
C
C
C
   50 NDM=2*MXNODS
      NO=N102 + (NEL - 1)*NDM*ITWO
      ND5DIM=MXNODS - 4
      NOO=N111 + (NEL - 1)*ND5DIM
C
      IF (IND.NE.0) GO TO 100
C
      NN=N110 + (NEL - 1)*IDW*NPT
      IF (NPAR(19).GT.0) NN=NN + (NEL - 1)*MXNODS
C
      CALL ICDMOD (NEL,NPT,A(NM),A(NN),A(NO),A(NOO),IA(NNOD),A(N6A1))
C
      GO TO 599
C
C
C
C
  100 IP=6*ITWO
      NMI=NM + IP
      IF (NPAR(15).EQ.5) NMI=NM + 8*ITWO
      NMI2=NMI + IP
      NMI3=NMI2 + IP
      NMI4=NMI3 + IP
C
      NN=N110 + (NEL - 1)*IDW*NPT + (IPT - 1)*IDW
      IF (NPAR(19).GT.0) NN=NN + (NEL - 1)*MXNODS
      NN4=NN + 4*ITWO
      NN8=NN + 8*ITWO
      NN9=NN8 + ITWO
      NN10=NN9 + ITWO
      NN11=NN10 + ITWO
      NN12=NN11 + ITWO
C
      CALL CDMOD (NEL,A(NM),A(NMI),A(NMI2),A(NMI3),A(NMI4),A(NN),A(NN4),
     1            A(NN8),A(NN9),A(NN10),A(NN11),A(NN12),STRESS,STRAIN,D,
     2            IPT,IPS,ISVE,IA(NNOD),A(N6A1),A(N6B1),A(NO),
     3            A(NOO),A(NN))
C
  599 CONTINUE
      RETURN
C
      END
      SUBROUTINE ICDMOD (NEL,NPT,PROP,WA,YZ,NOD5,NODS,TEMPV1)
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1           ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /TODIM/ BET,THIC,DE,IEL,NND5,ISOCOR
      COMMON /CRACK/ GAMMA,STIFAC,SHEFAC,SIGMAT,SIGMAC,EPSC,SIGMAU,EPSU,
     1               RAM5,RBM5,RCM5,ICRACK,MOD45
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      COMMON /EM2D/ S(378),XM(27),B(4,18),RE(27),EDIS(27),EDISI(27),
     1              XX(27),NOD(9),NODM(9),NOD5M(5)
      DIMENSION PROP(*),WA(15,*),YZ(*),NOD5(*),NODS(*),TEMPV1(*)
      DIMENSION H(9),P(2,9),XJ(2,2),PGRV(6)
C
      EQUIVALENCE (NPAR(17),NCON),(NPAR(10),NINT),(NPAR(15),MODEL),
     1            (NPAR(19),ITHERM)
C
C
C     INITIALIZE WA
C
      DO 11 J=1,NPT
      DO 10 I=1,15
   10 WA(I,J)=0.D0
   11 WA(12,J)=1.D3
      IF (ITHERM.EQ.0) GO TO 18
      II=0
      DO 15 K=1,9
      IF (NODM(K).EQ.0) GO TO 15
      II=II + 1
      NODS(II)=NODM(K)
   15 CONTINUE
      GO TO 20
   18 IF (MODEL.EQ.5) RETURN
      IF (ICRACK.LT.1) RETURN
C
C
C
C
      IPOINT=NCON/4 - 1
      PGRV(1)=0.D0
      DO 8 K=2,IPOINT
      DEV=PROP(K) - PROP(K-1)
      BULKK=(PROP(IPOINT+K) + PROP(IPOINT+K-1))/2.0
    8 PGRV(K)=PGRV(K-1) + BULKK*DEV
C
C
C
C
   20 IINTP=1
      DO 100 LX=1,NINT
      R1=XG(LX,NINT)
      DO 100 LZ=1,NINT
      S1=XG(LZ,NINT)
      IPT=(LX-1)*NINT + LZ
C
      CALL FUNCT2 (R1,S1,H,P,NOD5,XJ,DET,YZ,NEL,IINTP)
C
      IF (ITHERM.EQ.0) GO TO 25
      TMPOLD=0.D0
      DO 23 K=1,IEL
      KK=NODS(K)
   23 TMPOLD=TMPOLD + H(K)*TEMPV1(KK)
      WA(10,IPT)=TMPOLD
      IF (MODEX.EQ.2) GO TO 100
      TREF=PROP(39)
      XTOL=1.D-5
      TSUP=TREF*(1.+XTOL)
      TINF=TREF*(1.-XTOL)
      IF (TMPOLD.LE.TSUP .AND. TMPOLD.GE.TINF) GO TO 100
      write(66,2010) TMPOLD,TREF
      IF (MODEX .EQ. 0) GO TO 100
      STOP
C
   25 KK=0
      ZDEPTH=0.D0
      DO 30 K=1,IEL
      KK=KK + 2
   30 ZDEPTH=ZDEPTH + H(K)*YZ(KK)
      PGRAV=-GAMMA*ZDEPTH
      WA(11,IPT)=PGRAV
C
  100 CONTINUE
      IF (ITHERM.GT.0) RETURN
C
C
C
C
      DO 55 IPT=1,NPT
      PGRAV=WA(11,IPT)
      DO 35 L=2,IPOINT
      J=L
      IF (PGRAV.LT.PGRV(L)) GO TO 40
   35 CONTINUE
      write(66,2002) NEL
      STOP
   40 CONTINUE
      I=J - 1
      DD=PROP(J) - PROP(I)
      C1=PROP(IPOINT+I)
      C2=(PROP(IPOINT+J) - PROP(IPOINT+I))/(2.0*DD)
      IF (C2.GT.1.D-8) GO TO 41
      DEV=(PGRAV - PGRV(I))/C1
      GO TO 50
   41 FAC=C1*C1 + 4.0*C2*(PGRAV - PGRV(I))
      FAC=SQRT(FAC)
      DEV1=(-C1 + FAC)/(2.0*C2)
      DEV2=(-C1 - FAC)/(2.0*C2)
      I1=0
      I2=0
      IF (DEV1.GE.0.D0 .AND. DEV1.LE.DD) I1=1
      IF (DEV2.GE.0.D0 .AND. DEV2.LE.DD) I2=1
      IF (I1.NE.I2) GO TO 45
      write(66,2003)
      STOP
   45 IF (I1.EQ.1) DEV=DEV1
      IF (I2.EQ.1) DEV=DEV2
   50 EVGRAV=PROP(I) + DEV
      WA(10,IPT)=EVGRAV
   55 CONTINUE
C
C
      RETURN
C
C
 2002 FORMAT(55H **STOP - GRAVITATIONAL STRAIN TOO LARGE FOR MATERIAL C
     1       ,10HURVE INPUT,12H FOR ELEMENT,I5)
 2003 FORMAT(52H **STOP - ERROR IN PRESSURE-VOLUMETRIC STRAIN CALC.  )
 2010 FORMAT (44H ** STOP **-ERROR IN TEMPERATURES DEFINITION,/,
     150H INITIAL TEMP. MUST BE EQUAL TO TREF. FOR CONCRETE,/,
     19HT.INIT. =,E14.6,10X,6HTREF =,E14.6)
C
      END
      SUBROUTINE CDMOD (NEL,EKK,RKLD,RKUN,GLD,SP33,SIG,EPS,EVMAX,
     1                  EVGRAV,PGRAV,ANGLE,CRKSTR,STRESS,STRAIN,C,IPT,
     2                  IPS,ISVE,NODS,TEMPV1,TEMPV2,YZ,NOD5,WA)
C
C .                                                                   .
C .                                                                   .
C .                                                                   .
C .                                                                   .
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1           ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /CRACK/ GAMMA,STIFAC,SHEFAC,SIGMAT,SIGMAC,EPSC,SIGMAU,EPSU,
     1               RAM5,RBM5,RCM5,ICRACK,MOD45
      COMMON /CONCR2/ BETA,GAMA,RKAPA,ALFA,SIGP(4),TEP(4),EP(4),YP(3),
     1                E,VNU,RK,G,E12,E14,E24,EPSCP,SIGCP,FALSTR,ILFSET
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),F1,F2,F3
      COMMON /SIGPRN/ IPRNT
      COMMON /PORT/ INPORT,JNPORT,NPUTSV,IBPORT,IFPORT,JDC,JVC,JAC
C
      DIMENSION EKK(*),RKLD(*),RKUN(*),GLD(*),SP33(*),SIG(4),EPS(4),
     1          CRKSTR(3),STRESS(4),STRAIN(4),C(4,4),NODS(*),TEMPV1(*),
     2          TEMPV2(*),YZ(*),NOD5(*),WA(*),DUMWA(15)
      EQUIVALENCE (SIGP(1),P1),(SIGP(2),P2),(SIGP(4),P3)
C
C
      EQUIVALENCE (NPAR(5),ITYP2D),(NPAR(17),NCON),(NPAR(10),NINT)
     1           ,(NPAR(15),MODEL),(NPAR(19),ITHERM),(NPAR(3),INDNL)
C
C
      IF (MODEL.NE.5 .OR. ITHERM.NE.1) GO TO 300
      CALL THERM2 (NODS,NOD5,EVGRAV,TEMPV2,YZ,EKK,ITYP2D,1)
  300 CONTINUE
      IF (IUPDT.EQ.0) GO TO 2
      DO 1 I=1,15
    1 DUMWA(I)=WA(I)
C
C
    2 IF (MODEL.EQ.5) GO TO 3
      FALSTR=PGRAV
      EVTOT=EVMAX + EVGRAV
      IF (ITYP2D.LT.2) GO TO 3
      write(66,2005)
      STOP
    3 NPT=NINT*NINT
      IPOINT=6
      MOD45=1
      NUMCRK=0
      ANGPRI=1.D3
      ILFSET=0
C
C
C
C
      TMM=(SIG(1)+SIG(2)+SIG(4))/3.
      T11=SIG(1) - TMM
      T22=SIG(2) - TMM
      T33=SIG(4) - TMM
      T12=SIG(3)
C
      EVV=-(EPS(1)+EPS(2)+EPS(4))
      EMM=-EVV/3.
      E11=EPS(1) - EMM
      E22=EPS(2) - EMM
      E33=EPS(4) - EMM
      E12=EPS(3)
      SBAR=((SIG(1)-SIG(2))**2+(SIG(1)-SIG(4))**2+(SIG(2)-SIG(4))**2)/6.
      IF (MODEL.EQ.5) EVV=SQRT(SBAR + SIG(3)**2) + 3.*ALFA*TMM
C
C
C
C
      IKAS=1
      IF (ABS(EVV-EVMAX).GT.ABS(EVMAX)*1.D-8) IKAS=-1
      IK=1
      ITEM=1
C
C
    5 IF (MODEL.EQ.4) GO TO 10
C
C
      E=EKK(1)
      VNU=EKK(2)
      IF (ITYP2D.GE.2) STRAIN(4)=(STRAIN(1) + STRAIN(2))*VNU/(VNU - 1.)
C
C
      IF (PGRAV.NE.1.D2) GO TO 12
      IKAS=1
      E=E*STIFAC
      SIGCP=CRKSTR(1)
      EPSCP=CRKSTR(2)
      IF (EVV.EQ.0.D0) GO TO 45
      E=EKK(1)
C
C     LOADING
C
   12 IF (ITEM.NE.1) GO TO 7
      IF (IKAS.GT.0 .AND. (ANGLE.GT.3.61D2 .OR. PGRAV.EQ.1.D2))
     1 CALL PRNCPL (SIG,EPS,ANG,RKLD,RKUN,GLD,SP33,TEP,PGRAV,MODEL,1)
C
      IF (ANGLE.LT.3.61D2 .AND. PGRAV.NE.1.D2)
     1  CALL CRAKID (SIG,EPS,PGRAV,CRKSTR,RKLD,RKUN,GLD,SP33,
     2               ANGLE,TEP,NUMCRK,MODEL,1)
C
C
      IF (IKAS.GT.0 .AND. (ANGLE.GT.3.61D2 .OR. PGRAV.EQ.1.D2))
     1 CALL PRNCPL (SIG,STRAIN,ANG,RKLD,RKUN,GLD,SP33,EP,PGRAV,MODEL,2)
C
      IF (ANGLE.LT.3.61D2 .AND. PGRAV.NE.1.D2)
     1  CALL CRAKID (SIG,STRAIN,PGRAV,CRKSTR,RKLD,RKUN,GLD,SP33,
     2               ANGLE,EP,NUMCRK,MODEL,2)
      IF (IKAS.LT.0) GO TO 45
C
C
      NUMINT=3
    7 IF (ITEM.EQ.2) NUMINT=1
      DENM=ABS(P1) + ABS(P2) + ABS(P3)
      IF (DENM.LE.0.00001*SIGMAT) GO TO 45
C
      RP=EPSU/EPSC
      ES=SIGCP/EPSCP
      EU=(SIGMAU*SIGCP/SIGMAC)/(EPSU*EPSCP/EPSC)
      RAM5=EKK(1)/EU + (RP - 2.)*RP*RP*EKK(1)/ES - (2.*RP+1)*(RP-1.)**2
      RAM5=RAM5/(RP*(RP - 1.)**2)
      RBM5=2.*EKK(1)/ES - 3. - 2.*RAM5
      RCM5=2. - EKK(1)/ES + RAM5
C
      N1=1
      N2=4
      IF (PGRAV.NE.1.D2) GO TO 9
      N1=4
      IF (SIGP(2).LT.SIGP(4)) N1=2
      N2=N1
    9 DO 6 J=N1,N2
      IF (J.EQ.3) GO TO 6
      I=J
      IF (J.EQ.4) I=3
      YP(I)=E
      IF (SIGP(J).GE.0.001*SIGCP) GO TO 6
C
      YP(I)=0.D0
      DO 4 L=1,NUMINT
      E1=XG(L,NUMINT)
      DE=TEP(J) + (1 + E1)*(EP(J) - TEP(J))/2.
      DE=DE/EPSCP
      TY=E
      IF (DE.LE.0.D0) GO TO 4
      TY=E*(1. - RBM5*DE*DE - 2.*RCM5*DE**3)
      TY=TY/(1. + RAM5*DE + RBM5*DE*DE + RCM5*DE**3)**2
    4 YP(I)=YP(I) + 0.5*WGT(L,NUMINT)*TY
    6 CONTINUE
C
C
      IF (PGRAV.NE.1.D2) GO TO 8
      E=YP(I)
      IF (E.GT.0.D0) E=0.D0
      GO TO 45
C
C
    8 MOD45=1
      DE=RKAPA*SIGCP
      IF (P1.LT.DE .OR. P2.LT.DE .OR. P3.LT.DE) MOD45=2
      E=(ABS(P1)*YP(1) + ABS(P2)*YP(2) + ABS(P3)*YP(3))/DENM
      IF (MOD45.EQ.1) GO TO 45
C
C
      E12=E
      DENM=ABS(P1) + ABS(P2)
      IF (DENM.NE.0.D0) E12=(ABS(P1)*YP(1) + ABS(P2)*YP(2))/DENM
      E14=E
      DENM=ABS(P1) + ABS(P3)
      IF (DENM.NE.0.D0) E14=(ABS(P1)*YP(1) + ABS(P3)*YP(3))/DENM
      E24=E
      DENM=ABS(P2) + ABS(P3)
      IF (DENM.NE.0.D0) E24=(ABS(P2)*YP(2) + ABS(P3)*YP(3))/DENM
      GO TO 45
C
C
   10 DO 15 L=IK,IPOINT
      J=L
      IF (EVTOT.LT.EKK(J)) GO TO 28
   15 CONTINUE
      write(66,2015) NG,NEL,EVTOT,EKK(IPOINT)
      STOP
C
   28 I=J-1
C
      DELEV=EKK(J) - EKK(I)
      DELEI=EVTOT - EKK(I)
      RATIO=DELEI / DELEV
C
      IF (IKAS) 30,35,35
C
   30 RKUNLO=RKUN(I) + RATIO * (RKUN(J) - RKUN(I))
   35 RK =RKLD(I) + RATIO * (RKLD(J) - RKLD(I))
      G    =GLD(I) + RATIO * ( GLD(J) - GLD(I) )
C
      IF (IKAS) 40,45,45
C
   40 G=G * (RKUNLO/RK)
      RK=RKUNLO
C
   45 IF (ITEM.EQ.2) GO TO 100
C
C
C
C
      IF (MODEL.EQ.5) GO TO 46
      IF (ICRACK.LT.1) GO TO 50
      IF (ANGLE.GT.3.61D2) GO TO 50
      CALL CRAKID (SIG,EPS,PGRAV,CRKSTR,RKLD,RKUN,GLD,SP33,
     1             ANGLE,TEP,NUMCRK,MODEL,1)
C
      CALL CRAKID (SIG,STRAIN,PGRAV,CRKSTR,RKLD,RKUN,GLD,SP33,
     1             ANGLE,EP,NUMCRK,MODEL,2)
      GO TO 47
C
   46 RK=E/(3.*(1. - 2.*VNU))
      G=E/(2.*(1. + VNU))
      IF (PGRAV.EQ.1.D2) GO TO 50
      IF (ANGLE.GT.3.61D2 .AND. MOD45.EQ.1) GO TO 50
C
   47 CALL DCRACK (C,SIG,ANGLE,MODEL,ITYP2D,NUMCRK,1,1)