a21.for

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

      SUBROUTINE EL3D12
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      COMMON A(1)
      REAL A
      DIMENSION IA(1)
      EQUIVALENCE (A(1),IA(1))
C
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,
     1            ISTAT,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /DIMEL/ N101,N102,N103,N104,N105,N106,N107,N108,N109,
     1               N110,N111,N112,N113,N114,N120,N121,N122,N123,
     2               N124,N125
      COMMON /MTMD3D/ C(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /DPR/ ITWO
C
      EQUIVALENCE (NPAR(10),NINT),(NPAR(11),NINTZ),(NPAR(17),NCON),
     1            (NPAR(20),JDW),(NPAR(7),MXNODS)
C
C
      IDW=JDW*ITWO
      NPT=NINT*NINT*NINTZ
C
C
      MATP=IA(N107+NEL-1)
C
C
      NM=N111+(MATP-1)*NCON*ITWO
C
C
      NN=N112+(NEL-1)*(NPT*IDW+MXNODS)
C
      IF (IND.NE.0) GO TO 100
C
C
      CALL IUSER3 (A(NN),A(NN+IDW*NPT),A(NN),A(NM),IDW,JDW,
     1             A(N6A+ITWO))
      GO TO 599
C
C
C
C
  100 NS=NN+(IPT-1)*IDW
      NS1=NS+6*ITWO
      NS2=NS+12*ITWO
      NS3=NS+(JDW-2)*ITWO
C
C
      KK=NN+IDW*NPT
C
C
      ND9DIM=MXNODS-8
      LL=N108+(NEL-1)*ND9DIM
C
C
      CALL USER3 (A(NM),A(NS),A(NS1),A(NS2),A(NS3),A(KK),A(LL),
     1            A(N6A+ITWO),A(N6B+ITWO))
C
C
  599 CONTINUE
      RETURN
      END
      SUBROUTINE IUSER3 (WA,IIWA,IWA,PROP,IIDW,JDW,TEMPV1)
C
C
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /DPR/ ITWO
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,
     1            ISTAT,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,IREF,
     A             IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /ELSTP/ TIME,IDTHF
      COMMON /CONST/ DT,DTA,CONS(21),DTOD,IOPE
      COMMON /EM3D/ NOD(27),NODM(27),NOD9M(19)
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /THREED/ DE,IELD,IELX,IEL,NPT,IDW,NND9,ISOCOR
      COMMON /MTMD3D/ C(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
C
      DIMENSION WA(JDW,*),IWA(IIDW,*),PROP(16,*),DUMMY(23),IDUMMY(2),
     A          IIWA(*),TEMPV1(*),H(27),XDM1(3,27),XDM2(3,3),XDM3(3,1)
      DIMENSION EPS(6),DEPS(6),CTD(5),CTI(8),SCP(5)
      DIMENSION CTDD(5),DTHSTR(6),THSTR1(6),THSTR2(6)
C
      EQUIVALENCE (NPAR(10),NINT),(NPAR(11),NINTZ)
C
C
      II=0
      DO 5 K=1,27
      IF ( NODM(K).EQ.0 ) GO TO 5
      II=II+1
      IIWA(II)=NODM(K)
    5 CONTINUE
C
C
      IINTP=1
      IPT=0
      DO 30 LX=1,NINT
      E1=XG(LX,NINT)
      DO 30 LY=1,NINT
      E2=XG(LY,NINT)
      DO 30 LZ=1,NINTZ
      E3=XG(LZ,NINTZ)
      IPT=IPT+1
C
      CALL FUNCT (E1,E2,E3,H,XDM1,NOD9M,XDM2,XDUM,XDM3,IINTP)
C
      TEMP1=0.D0
      DO 35 K=1,IEL
      KK=IIWA(K)
   35 TEMP1=TEMP1+H(K)*TEMPV1(KK)
   30 WA(13,IPT)=TEMP1
C
C
      NNG=NG
      NNEL=NEL
      IIPT=IPT
      DO 10 K=1,5
   10 SCP(K)=PROP(K+11,8)
      CALL INTMAT (NNG,NNEL,IIPT,PROP,TEMP1,ALFA,CTD,CTI,1)
C
C
      KEY=1
      TIM=TIME
      DDT=DT
      TEMP2=TEMP1
      DO 11 K=1,23
   11 DUMMY(K)=0.D0
      DO 12 K=1,2
   12 IDUMMY(K)=0
      INTER=INT( PROP(11,8) )
      KR=INTER
C
      CALL CUSER3 (NNG,NNEL,IIPT,STRESS,EPS,STRAIN,DEPS,DTHSTR,
     A             THSTR1,THSTR2,KR,INTER,SCP,DUMMY,IDUMMY,C,ALFA,CTD,
     B             ALFAA,CTDD,CTI,TEMP1,TEMP2,TIM,DDT,KEY)
C
      NPT=NINT*NINT*NINTZ
      JJDW=JDW-2
C
      DO 25 J=1,NPT
C
C
      DO 15 I=1,12
   15 WA(I,J)=0.D0
C
      DO 20 II=14,JJDW
   20 WA(II,J)=DUMMY(II-13)
C
      KJ=JJDW*ITWO+1
      KJJ=KJ + ITWO
      IWA(KJ,J)=IDUMMY(1)
   25 IWA(KJJ,J)=IDUMMY(2)
C
C
      RETURN
      END
      SUBROUTINE INTMAT (NG,NEL,IPT,PROP,TMP,ALFA,CTD,CTI,KKK)
C
C
C     TEMPERATURE
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      DIMENSION PROP(16,*),CTD(5),CTI(8)
C
C
      TTOL=1.D-10
      NPTS=INT(PROP(9,8))
      L=0
      DO 20 K=2,NPTS
      L=L+1
      TDIFF=TMP - PROP(K,1)
      IF (ABS(TDIFF).LT.TTOL) GO TO 25
      IF (TDIFF.LT.0.D0) GO TO 25
   20 CONTINUE
C
      write(66,3000) NG,NEL,IPT,TMP
      STOP
C
   25 XRATIO=(TMP-PROP(L,1))/(PROP(L+1,1)-PROP(L,1))
C
C
      IF (XRATIO.GT.1.D0) XRATIO=1.D0
      IF (XRATIO.LT.0.D0) XRATIO=0.D0
C
C
      ALFA=PROP(L,2) + XRATIO*( PROP(L+1,2) - PROP(L,2) )
C
      IF (KKK.EQ.0) RETURN
C
      DO 30 M=3,7
   30 CTD(M-2)=PROP(L,M) + XRATIO*( PROP(L+1,M) - PROP(L,M) )
C
C
      DO 31 N=1,8
   31 CTI(N)=PROP(N,8)
C
      RETURN
C
 3000 FORMAT (//,15H *** ERROR *** ,
     A         /,30H ELEMENT GROUP NUMBER         ,I5,
     B         /,30H ELEMENT NUMBER               ,I5,
     C         /,30H INTEGRATION POINT NUMBER     ,I5,
     D         /,30H INTERPOLATED TEMPERATURE     ,E14.6,
     E         /,51H TEMPERATURE LIES OUTSIDE RANGE USED IN THE MATERIA,
     F           51HL CONSTANTS DEFINITION                             ,
     G         /,13H *** STOP ***                                      )
      END
      SUBROUTINE USER3 (PROP,SIG,EPS,ARRAY,IARRAY,NDS,NOD9M,TEMPV1,
     1                  TEMPV2)
C
C
C .                                                                    .
C .                                                                    .
C .                                                                    .
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /CONST/ DT,DTA,CONS(21),DTOD,IOPE
      COMMON /THREED/ DE,IELD,IELX,IEL,NPT,IDW,NND9,ISOCOR
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,
     1            ISTAT,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,IREF,
     1             IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /MTMD3D/ C(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
      COMMON /DISDR/ DISD(9)
      COMMON /ELSTP/ TIME,IDTHF
      COMMON /PORT/ INPORT,JNPORT,NPUTSV,IBPORT,IFPORT,JDC,JVC,JAC
C
      EQUIVALENCE (NPAR(20),JDW),(NPAR(3),INDNL)
C
      DIMENSION PROP(16,*),SIG(*),EPS(*),ARRAY(*),IARRAY(*),
     1          H(27),XDM1(3,27),XDM2(3,3),XDM3(3,1),TEMPV2(*),
     2          TEMPV1(*),NDS(*),NOD9M(*),DEPS(6),DEPSS(6),DUMMY(23),
     3          IDUMMY(2),CTD(5),CTI(8),SCP(5)
      DIMENSION CTDD(5),DTHSTR(6),THSTR1(6),THSTR2(6)
C
      TREF=PROP(10,8)
      INTER=INT( PROP(11,8) )
      XINTER=FLOAT(INTER)
      DO 5 I=1,5
    5 SCP(I)=PROP(I+11,8)
C
C
      IINTP=1
      CALL FUNCT (E1,E2,E3,H,XDM1,NOD9M,XDM2,XDUM,XDM3,IINTP)
C
C     SOLUTION STEP
C
      TEMP1=0.D0
      TEMP2=0.D0
      DO 10 K=1,IEL
      KK=NDS(K)
      TEMP1=TEMP1+H(K)*TEMPV1(KK)
   10 TEMP2=TEMP2+H(K)*TEMPV2(KK)
C
      TMPOLD=ARRAY(1)
      CTEMP=TEMP2
      IF (IEQUIT.EQ.1 .AND. KPRI.NE.0) CTEMP=TEMP1
      DTMP=( CTEMP - TMPOLD ) / XINTER
C
C
C
      DO 100 I=1,6
  100 STRESS(I)=SIG(I)
      NNG=NG
      NNEL=NEL
      IIPT=IPT
      TIM=TIME
      DDT=DT
      JJDW=JDW - 14
      JDW1=JJDW - 1
C
      DO 130 I=1,JDW1
  130 DUMMY(I)=ARRAY(I+1)
      DO 135 I=1,2
  135 IDUMMY(I)=IARRAY(I)
C
      DO 150 I=1,6
  150 DEPSS(I)=( STRAIN(I) - EPS(I) ) / XINTER
      DO 155 I=4,6
  155 DEPS(I)=DEPSS(I)
C
C
C
      KEY=2
      TMP1=TMPOLD
C
      DO 200 K=1,INTER
      KR=K
      TMP2=TMP1 + DTMP
C
      CALL INTMAT (NNG,NNEL,IIPT,PROP,TMP2,ALFAA,CTDD,CTI,1)
      EPST2=ALFAA*( TMP2 - TREF )
      CALL INTMAT (NNG,NNEL,IIPT,PROP,TMP1,ALFA,CTD,CTI,1)
      EPST1=ALFA*( TMP1 - TREF )
      DEPST=EPST2 - EPST1
C
      DO 160 I=1,3
  160 DEPS(I)=DEPSS(I) - DEPST
C
      DO 165 J=1,6
      THSTR1(J)=0.D0
      THSTR2(J)=0.D0
  165 DTHSTR(J)=0.D0
C
      DO 170 J=1,3
      THSTR1(J)=EPST1
      THSTR2(J)=EPST2
  170 DTHSTR(J)=EPST2 - EPST1
C
      CALL CUSER3 (NNG,NNEL,IIPT,STRESS,EPS,STRAIN,DEPS,DTHSTR,
     A             THSTR1,THSTR2,KR,INTER,SCP,DUMMY,IDUMMY,C,ALFA,CTD,
     B             ALFAA,CTDD,CTI,TMP1,TMP2,TIM,DDT,KEY)
C
  200 TMP1=TMP1 + DTMP
C
C
C
C
      IF (IUPDT.NE.0) GO TO 210
C
      DO 220 J=1,6
      SIG(J)=STRESS(J)
  220 EPS(J)=STRAIN(J)
C
      ARRAY(1)=TMP2
C
C
      DO 230 I=2,JJDW
  230 ARRAY(I)=DUMMY(I-1)
C
      IARRAY(1)=IDUMMY(1)
      IARRAY(2)=IDUMMY(2)
C
C
C
  210 IF (KPRI.EQ.0) GO TO 240
      IF (ICOUNT.EQ.3) RETURN
C
C
      KEY=3
      CALL INTMAT (NNG,NNEL,IIPT,PROP,TEMP2,ALFA,CTD,CTI,1)
      CALL CUSER3 (NNG,NNEL,IIPT,STRESS,EPS,STRAIN,DEPS,DTHSTR,
     A             THSTR1,THSTR2,KR,INTER,SCP,DUMMY,IDUMMY,C,ALFA,CTD,
     B             ALFAA,CTDD,CTI,TMP1,TMP2,TIM,DDT,KEY)
C
      RETURN
C
C
C
  240 IF (IPRI.NE.0 .OR. IPS.EQ.0) GO TO 255
C
      KEY=4
C
C
      IF (INDNL.EQ.2) CALL CAUCH3
      CALL CUSER3 (NNG,NNEL,IIPT,STRESS,EPS,STRAIN,DEPS,DTHSTR,
     A             THSTR1,THSTR2,KR,INTER,SCP,DUMMY,IDUMMY,C,ALFA,CTD,
     B             ALFAA,CTDD,CTI,TMP1,TMP2,TIM,DDT,KEY)
C
C
  255 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 260
      IF (ISVE.EQ.0) GO TO 260
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-3',NEL,IPT,(STRESS(I),I=1,6),
     2                       (STRAIN(I),I=1,6)
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9000) 'OUTPUT-3',NEL,IPT,(STRESS(I),I=1,6),
     2                       (STRAIN(I),I=1,6)
C
 9000 FORMAT ( A,/,2I10,/,(4E20.13) )
C
C
  260 CONTINUE
      RETURN
C
C
C
      END
      SUBROUTINE XJ123 (STRESS,XI1,XI2,XI3,XJ1,XJ2,XJ3)
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
C
C
C
      DIMENSION STRESS(6)
C
      S11 = STRESS(1)
      S22 = STRESS(2)
      S33 = STRESS(3)
      S12 = STRESS(4)
      S31 = STRESS(5)
      S23 = STRESS(6)
C
      SMEAN = ( S11+S22+S33 )/3.0
C
      D11 = S11 - SMEAN
      D22 = S22 - SMEAN
      D33 = S33 - SMEAN
      D12 = S12
      D31 = S31
      D23 = S23
C
      XI1 = 3.0*SMEAN
      XJ1 = 0.D0
C
      XI2 = 0.5*( S11*S11 + S22*S22 + S33*S33 ) + S23*S23
     1     + S12*S12 + S31*S31
      XJ2 = 0.5*( D11*D11 + D22*D22 + D33*D33 ) + D23*D23
     1     + D12*D12 + D31*D31
C
      XI3 = (1./3.)*( S11*S11*S11 + S22*S22*S22 + S33*S33*S33 )
     1    + S11*( S12*S12 + S31*S31 ) + S22*( S12*S12 + S23*S23 )
     2    + S33*( S23*S23 + S31*S31 ) + 2.*S12*S31*S23
      XJ3 = (1./3.)*( D11*D11*D11 + D22*D22*D22 + D33*D33*D33 )
     1    + D11*( D12*D12 + D31*D31 ) + D22*( D12*D12 + D23*D23 )
     2    + D33*( D23*D23 + D31*D31 ) + 2.*D12*D31*D23
C
      RETURN
      END
      SUBROUTINE CUSER3 (NG,NEL,IPT,STRESS,EPS,STRAIN,DEPS,DEPST,
     A                   THSTR1,THSTR2,KTR,INTER,SCP,ARRAY,IARRAY,D,
     B                   ALFA,CTD,ALFAA,CTDD,CTI,TMP1,TMP2,TIME,DT,KEY)
C*I
C*I
C*I
C*I
C*I
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      DIMENSION STRESS(6),STRAIN(6),DEPS(6),ARRAY(23),IARRAY(2),D(6,6),
     A          CTD(5),CTI(8),EPS(6),SCP(5)
      DIMENSION CTDD(5),DEPST(6),THSTR1(6),THSTR2(6)
C
      GO TO (1,2,3,4), KEY
C*I
C*I
C*I    K E Y  =  1
C*I
C*I
    1 CONTINUE
C*I
C*I
C*I
      RETURN
C*I
C*I
C*I   K E Y = 2
C*I
C*I
    2 CONTINUE
C*I
C*I
C*I
      RETURN
C*I
C*I
C*I
C*I  K E Y = 3
C*I
    3 CONTINUE
C*I
C*I
C*I
      RETURN
C*I
C*I
C*I
C*I   K E Y = 4
C*I
    4 CONTINUE
C*I
C*I
C*I
      RETURN
C*FILE END
      END