a04b.for

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

 2010 FORMAT(45H0***ERROR  MATERIAL PROPERTIES NOT ADMISSABLE  )
C
      END
      SUBROUTINE POSINV (A,NMAX,NDD)
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      DIMENSION A(NDD,NDD)
C
      DO 200 N=1,NMAX
C
      D=A(N,N)
      DO 100 J=1,NMAX
  100 A(N,J)=-A(N,J)/D
C
      DO 150 I=1,NMAX
      IF(N-I) 110,150,110
  110 DO 140 J=1,NMAX
      IF(N-J) 120,140,120
  120 A(I,J)=A(I,J)+A(I,N)*A(N,J)
  140 CONTINUE
  150 A(I,N)=A(I,N)/D
C
      A(N,N)=1.0/D
C
  200 CONTINUE
C
      RETURN
      END
      SUBROUTINE STSTN (XX,PROP,DISD,IDW,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 /MATMOD/ STRESS(4),STRAIN(4),C(4,4),IPT,NEL,IPS,ISVE
C
      DIMENSION WA(IDW,*),XX(2,*),PROP(*),DISD(*),DN(4)
C
      EQUIVALENCE (NPAR(5),ITYP2D),(NPAR(15),MODEL),(NPAR(3),INDNL)
     1           ,(NPAR(17),NCON)
C
C
      IST=3
      IF (ITYP2D.EQ.0) IST=4
C
C
C
C
C
      STRAIN(1)=DISD(1)
      STRAIN(2)=DISD(2)
      STRAIN(3)=DISD(3) + DISD(4)
      STRAIN(4)=0.D0
      IF (ITYP2D.EQ.0) STRAIN(4)=DISD(5)
      IF (INDNL.LE.1) GO TO 80
C
C
      DN(1)=0.5*(DISD(1)*DISD(1) + DISD(4)*DISD(4))
      DN(2)=0.5*(DISD(2)*DISD(2) + DISD(3)*DISD(3))
      DN(3)=DISD(3)*DISD(1) + DISD(4)*DISD(2)
      IF (IST.EQ.4) DN(4)=0.5*DISD(5)*DISD(5)
C
      IF (INDNL.EQ.3) GO TO 60
C
C
      DO 20 I=1,IST
   20 STRAIN(I)=STRAIN(I) + DN(I)
      GO TO 80
C
C
   60 IF (MODEL.NE.1) GO TO 80
      DO 40 I=1,IST
   40 STRAIN(I)=STRAIN(I) - DN(I)
C
C
C
C
C
   80 GO TO (1,2,3,4,4,6,7,8,8,10,10,12,13,14,14,14) ,MODEL
C
C
C
C
    1 A1=C(1,1)
      B1=C(1,2)
C
      STRESS(1) = A1*STRAIN(1) + B1*STRAIN(2)
      STRESS(2) = B1*STRAIN(1) + A1*STRAIN(2)
      STRESS(3) = C(3,3)*STRAIN(3)
      STRESS(4) = 0.D0
C
      A2=-PROP(2)/(1. - PROP(2))
      IF (ITYP2D.GE.2) STRAIN(4)=A2*(STRAIN(1) + STRAIN(2))
      IF (ITYP2D.GE.2) GO TO 110
C
      STRESS(4) = B1*(STRAIN(1)+STRAIN(2))
      IF (ITYP2D.EQ.1) GO TO 110
C
C
      STRESS(1) = STRESS(1) + B1*STRAIN(4)
      STRESS(2) = STRESS(2) + B1*STRAIN(4)
      STRESS(4) = STRESS(4) + A1*STRAIN(4)
C
  110 RETURN
C
C
C
C
    2 IF (PROP(3).EQ.0.D0) GO TO 1
C
      DO 102 I=1,4
  102 STRESS(I)=0.D0
C
C
      DO 103 J=1,IST
      DO 103 I=1,IST
  103 STRESS(I) = STRESS(I) + C(I,J)*STRAIN(J)
      IF (ITYP2D.GE.2) STRAIN(4)=
     1  -(C(4,1)*STRAIN(1)+C(4,2)*STRAIN(2)+C(4,3)*STRAIN(3))/C(4,4)
      IF (ITYP2D.NE.1) GO TO 120
C
C
      STRESS(4)=C(4,1)*STRAIN(1) + C(4,2)*STRAIN(2) + C(4,3)*STRAIN(3)
C
  120 RETURN
C
C
C
    3 CALL ELT2D3
      RETURN
C
C
C
C
    4 CALL ELT2D4
      RETURN
C
C
C
    6 CALL ELT2D6
      RETURN
C
C
C
    7 CALL ELT2D7
      RETURN
C
C
C
    8 CALL ELT2D8
      RETURN
C
C
C
   10 CALL EL2D10
      RETURN
C
C
C
   12 CALL EL2D12
      RETURN
C
C
C
   13 CALL EL2D13
      RETURN
C
C
C
   14 CALL EL2D14
      RETURN
C
C
      END
      SUBROUTINE CAUCHY
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 /DISDER/ DISD(5)
      COMMON /MATMOD/ STRESS(4),STRAIN(4),D(4,4),IPT,NEL,IPS,ISVE
      EQUIVALENCE (NPAR(5),ITYP2D)
C
C
      IF (ITYP2D-1) 710,720,730
  710 X33=1.+DISD(5)
      GO TO 750
  720 X33=1.D0
      GO TO 750
  730 X33=SQRT (1.+2.*STRAIN(4))
C
  750 X11=1.+DISD(1)
      X12=   DISD(3)
      X22=1.+DISD(2)
      X21=   DISD(4)
C
      DET=X33*(X11*X22 - X12*X21)
      CALL DETCH (NG,NEL,DET,1)
C
  760 DET=1./DET
      S1=STRESS(1)
      S2=STRESS(2)
      SS=STRESS(3)
C
      STRESS(1)=DET * (S1*X11*X11 + 2.*SS*X11*X12 + S2*X12*X12)
      STRESS(2)=DET * (S1*X21*X21 + 2.*SS*X21*X22 + S2*X22*X22)
      STRESS(3)=DET * (S1*X11*X21 + SS*(X11*X22+X12*X21) + S2*X12*X22)
      IF (ITYP2D.GE.2) RETURN
      STRESS(4)=DET * (STRESS(4)*X33*X33)
C
      RETURN
C
C
C
      END
      SUBROUTINE CGDT2 (YZ,EDIS,NDPN,NOD5,E1,E2,IDEATH,EDISB,IEL,
     1                  ITYP2D)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /PSTCH/ STRCH(3),RDCS(3)
      COMMON /MATMOD/ STRESS(4),STRAIN(4),D(4,4),IPT,NEL,IPS,ISVE
C
      DIMENSION CG(3,3),F(3,3),RLMN(3,3),S(3,3),IPRM(3),ICOL(3),
     1          PV(3),DISD(5),B(4,18),YZ(*),EDISB(*),EDIS(*),NOD5(*),
     2          XX(27),DUMMY1(27),DUMMY2(378)
C
      DATA IPRM /2,3,1/
C
C
      ND=NDPN*IEL
      DO 500 I=1,ND
  500 XX(I)=YZ(I)
      IF (IDEATH.NE.1) GO TO 510
C
      DO 505 I=1,ND
  505 XX(I)=XX(I) + EDISB(I)
C
  510 ND=2*IEL
      IF (NDPN.EQ.3) CALL PLST3D (YZ,XX,DUMMY1,EDIS,DUMMY2,IEL,3)
C
      CALL DERIQ (NEL,XX,B,DET,E1,E2,XBAR,NOD5)
C
      DO 515 I=1,5
  515 DISD(I)=0.D0
C
      DO 520 J=2,ND,2
      I=J - 1
      DISD(1)=DISD(1) + B(1,I)*EDIS(I)
      DISD(2)=DISD(2) + B(2,J)*EDIS(J)
      DISD(3)=DISD(3) + B(3,I)*EDIS(I)
  520 DISD(4)=DISD(4) + B(3,J)*EDIS(J)
C
      IF (ITYP2D.NE.0) GO TO 550
      DO 525 I=1,ND,2
  525 DISD(5)=DISD(5) + B(4,I)*EDIS(I)
C
C
  550 IF (ITYP2D - 1) 1,2,3
C
    1 F(3,3)=DISD(5) + 1.0
      GO TO 4
    2 F(3,3)=1.D0
      GO TO 4
    3 F(3,3)=SQRT(1.0 + 2.0*STRAIN(4))
C
    4 F(1,1)=DISD(1) + 1.0
      F(1,2)=DISD(3)
      F(1,3)=0.D0
      F(2,1)=DISD(4)
      F(2,2)=DISD(2) + 1.0
      F(2,3)=0.D0
      F(3,1)=0.D0
      F(3,2)=0.D0
C
C
C
      DO 5 I=1,3
      DO 5 J=1,3
    5 CG(I,J)=0.D0
C
      DO 15 I=1,3
      DO 15 J=I,3
      DO 10 M=1,3
   10 CG(I,J)=CG(I,J) + F(M,I)*F(M,J)
   15 CG(J,I)=CG(I,J)
C
C
C
C
C
      STRAV=(CG(1,1) + CG(2,2) + CG(3,3))/3.0
      DO 18 I=1,3
      DO 18 J=1,3
   18 S(I,J)=CG(I,J)
C
      S(1,1)=S(1,1) - STRAV
      S(2,2)=S(2,2) - STRAV
      S(3,3)=S(3,3) - STRAV
C
C
      SBAR=0.D0
      DO 20 I=1,3
      DO 20 J=1,3
   20 SBAR=SBAR + S(I,J)*S(I,J)
      SBAR=SQRT(SBAR/2.)
      SBTOL=ABS(STRAV)*1.D-8
      IF (SBAR.LE.SBTOL) SBAR=0.D0
      IF (SBAR.EQ.0.D0) GO TO 31
C
C
      RJ3=0.D0
      DO 30 I=1,3
      DO 30 J=1,3
      DO 30 K=1,3
   30 RJ3=RJ3 + S(I,J)*S(J,K)*S(K,I)
      RJ3=RJ3/3.
C
C
      TEMP=-3.D0*SQRT(3.D0)/2.D0
      TEMP=TEMP*RJ3/SBAR**3
   31 IF (SBAR.EQ.0.D0) TEMP=0.D0
      IF (ABS(TEMP) .LE. 1.D0) GO TO 32
      IF (ABS(TEMP) .LE. 1.0001D0) GO TO 34
      write(66,2000)
      STOP
C
   34 IF (TEMP .LT. (-1.D0)) TEMP=-1.D0
      IF (TEMP .GT. 1.D0) TEMP=1.D0
   32 PI=4.D0*ATAN(1.D0)
      TEMPCS=SQRT(1.-TEMP*TEMP)
      PHI=ATAN2(TEMP,TEMPCS)
      IF (PHI.GT.PI) PHI=PHI - 2.*PI
      PHI=PHI/3.0
C
C
      A1=2.*SBAR/SQRT(3.D0)
      PV(1)   =A1*SIN(PHI + 2.*PI/3.) + STRAV
      PV(2)   =A1*SIN(PHI) + STRAV
      PV(3)   =A1*SIN(PHI + 4.*PI/3.) + STRAV
C
C
      TOL=0.1D-1
      IND=0
      SPREAD=PV(1) - PV(3)
      ROERR=ABS(PV(1))*0.000001
      IF (PV(1).EQ.0.D0) ROERR=ABS(PV(3))*0.1D-5
      IF (SPREAD.LE.ROERR) GO TO 82
      DIF1=PV(1) - PV(2)
      DIF2=PV(2) - PV(3)
      IF (DIF1.LT.SPREAD*TOL) IND=3
      IF (DIF2.LT.SPREAD*TOL) IND=1
C
      N=3
      I1=IND
      IF (IND.EQ.0) I1=1
      NEIG=2
      DO 78 NX=1,NEIG
C
      DO 35 J1=1,N
      ICOL(J1)=J1
      S(J1,J1)=CG(J1,J1) - PV(I1)
   35 RLMN(J1,I1)=0.D0
C
C
      DO 65 J=1,N
C
      PIVI=0.D0
      DO 40 I=J,N
      DO 40 K=J,N
      IF (ABS(PIVI).GT.ABS(S(I,K))) GO TO 40
      IMAX=I
      KMAX=K
      PIVI=S(I,K)
   40 CONTINUE
      IF (KMAX.EQ.J) GO TO 47
C
C
      ISAVE=ICOL(KMAX)
      ICOL(KMAX)=ICOL(J)
      ICOL(J)=ISAVE
      DO 45 JJ=1,N
      SAVE=S(JJ,KMAX)
      S(JJ,KMAX)=S(JJ,J)
   45 S(JJ,J)=SAVE
C
   47 PIVI=1./PIVI
C
C
      DO 50 K=J,N
      IF (IMAX.EQ.J) GO TO 50
      SAVE=S(J,K)
      S(J,K)=S(IMAX,K)
      S(IMAX,K)=SAVE
   50 S(J,K)=S(J,K)*PIVI
C
      IF (J.EQ.(N-1)) GO TO 70
      I2=J + 1
      DO 65 K2=I2,N
      DO 65 J2=I2,N
   65 S(K2,J2)=S(K2,J2) - S(K2,J)*S(J,J2)
C
   70 N1=N - 1
      RLMN(ICOL(N),I1)=1.D0
      DO 75 J=1,N1
      IA=N
      DO 75 K=1,J
      RLMN(ICOL(N-J),I1)=RLMN(ICOL(N-J),I1)- S(N-J,IA)*RLMN(ICOL(IA),I1)
   75 IA=IA - 1
C
C
      RMAX=0.D0
      DO 74 L1=1,3
   74 IF (ABS(RLMN(L1,I1)).GT.RMAX) RMAX=ABS(RLMN(L1,I1))
      RMAX=RMAX*0.0001
      DO 76 L1=1,3
   76 IF (ABS(RLMN(L1,I1)).LE.RMAX) RLMN(L1,I1)=0.D0
      XLN=RLMN(1,I1)**2 + RLMN(2,I1)**2 + RLMN(3,I1)**2
      XLN=1.0/SQRT(XLN)
      IF (RLMN(3,I1).LT.0.D0) XLN=-XLN
      DO 77 J1=1,3
   77 RLMN(J1,I1)=RLMN(J1,I1)*XLN
C
      IF (I1.EQ.1) GO TO 100
      I1=3
      IF (IND.EQ.3) I1=1
   78 CONTINUE
C
C
      IND=3
      I2=IPRM(IND)
      I1=IPRM(I2)
      X=RLMN(2,IND)*RLMN(3,I2) - RLMN(3,IND)*RLMN(2,I2)
      Y=RLMN(3,IND)*RLMN(1,I2) - RLMN(1,IND)*RLMN(3,I2)
      Z=RLMN(1,IND)*RLMN(2,I2) - RLMN(2,IND)*RLMN(1,I2)
      XLN=1./SQRT(X*X + Y*Y + Z*Z)
      IF (Z.LT.0.D0) XLN=-XLN
      RLMN(1,I1)=X*XLN
      RLMN(2,I1)=Y*XLN
      RLMN(3,I1)=Z*XLN
      GO TO 84
C
   82 DO 83 I1=1,3
      DO 83 J1=1,3
      RLMN(J1,I1)=0.D0
   83 IF (I1.EQ.J1) RLMN(J1,I1)=1.D0
   84 CONTINUE
C
C
  100 STRCH(1)=SQRT(PV(1))
      STRCH(2)=SQRT(PV(2))
      STRCH(3)=SQRT(PV(3))
C
      DO 110 I=1,3
  110 RDCS(I)=RLMN(I,1)
C
      RETURN
C
 2000 FORMAT (///,106H   ERROR   UNABLE TO CALCULATE THE EIGENVALUES OF
     1THE CAUCHY-GREEN DEFORMATION TENSOR   (SUBROUTINE CGDT2))
C
C*FILE END
      END