a25b.for

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

 2120 FORMAT (1H ,4X,4(E14.6,7X))
 2125 FORMAT (//,4X,46HNUMBER OF TEMPERATURE POINTS. . .(PROP(65))..=,
     1        E14.6,/,
     2        4X,46HREFERENCE TEMPERATURE. . . . . . (PROP(66))..=,
     3        E14.6,//)
 2130 FORMAT (/ 42H    INCORRECT NUMBER OF TEMPERATURE POINTS)
 2140 FORMAT (/ 35H    TEMPERATURE POINTS OUT OF ORDER)
 2111 FORMAT (1H ,4X,29HE ............( PROP(1) ).. =,E14.6,/,
     1        1H ,4X,29HVNU ..........( PROP(2) ).. =,E14.6,/,
     2        1H ,4X,29HYIELD ........( PROP(3) ).. =,E14.6,//)
 2112 FORMAT (1H ,4X,36HPIECEWISE-LINEAR STRESS-STRAIN CURVE,/,
     1        1H ,6X,6HSTRESS,10X,6HSTRAIN,12X,2HET,//,
     2        6X,E14.6,2X,E14.6)
 2113 FORMAT (6X,3(E14.6,2X))
 2114 FORMAT (6X,E14.6,2X,E14.6,2X)
 2400 FORMAT (24H ORTHOTROPIC PROPERTIES ,//,
     1        16H MODULUS  EAA = ,E15.6,/,
     2        16H MODULUS  EAB = ,E15.6,/,
     3        16H MODULUS  EBB = ,E15.6,/,
     4        16H MODULUS  GAB = ,E15.6,/,
     5        16H MODULUS  GT  = ,E15.6,//)
 3000 FORMAT (//37H * * *     STOP OF SOLUTION     * * *,//,
     1          35H INPUT ERROR IN MATERIAL PROPERTIES,/,
     2          19H ELEMENT GROUP NO =,I5,/,
     3          27H MATERIAL PROPERTY SET NO =,I5)
 3050 FORMAT (/ 35H    YOUNG*S MODULUS MUST BE .GT. 0.,/,
     1          49H    POISSON*S RATIO MUST BE .GT. -1. AND .LT. 0.5)
 3100 FORMAT (/ 40H SEE USER*S MANUAL FOR INPUT DESCRIPTION,//,
     1          37H * * *   END OF ERROR MESSAGE   * * *,//)
 3150 FORMAT (/ 44H   FOR ALL TEMPERATURES THE MODEL REQUIRES :)
 3401 FORMAT (/ 35H    YOUNG*S MODULUS MUST BE .GT. 0.,/,
     1          49H    POISSON*S RATIO MUST BE .GT. -1. AND .LT. 0.5,/,
     2          50H    HARDENING MODULUS MUST BE .LT. YOUNG*S MODULUS,/,
     *          37H    HARDENING MODULUS MUST BE .GE. 0.,/,
     3          40H    INITIAL YIELD STRESS MUST BE .GT. 0.)
 3402 FORMAT (//50H INPUT ERROR DETECTED IN (MATRIT/SHELL)            //
     1          19H ELEMENT GROUP NO = ,I5/
     2          27H MATERIAL PROPERTY SET NO = ,I5/
     3          44H HARDENING MODULUS (ET) GREATER OR EQUAL TO  ,
     4          44H YOUNG*S MODULUS (E) IS NOT ALLOWED          //)
 3403 FORMAT (//50H INPUT ERROR IN MATERIAL PROPERTIES                //
     1          15H  *** STOP ***   //)
 3404 FORMAT (/ 48H    FOR THE MULTI-LINEAR ELASTIC-PLASTIC MODEL  ,/,
     1           9H    PROP(,I2,16H) EQUAL TO PROP(,I2,16H) IS NOT ALLOW
     2ED,/,
     3          35H    BECAUSE ET CANNOT BE CALCULATED)
 3405 FORMAT (/ 46H    FOR THE MULTI-LINEAR ELASTIC-PLASTIC MODEL,/,
     1          24H    WE MUST HAVE THAT : ,/,
     2          51H       E .GT. ET1 .GE. ET2 .GE. ET3 ....... .GE. 0.)
 3410 FORMAT (//39H INPUT ERROR DETECTED IN (MATRIT/SHELL),//,
     1        19H ELEMENT GROUP NO = ,I5,/,
     2        27H MATERIAL PROPERTY SET NO = ,I5,/,
     3        22H ET.LT.0.0 NOT ALLOWED  )
 3412 FORMAT (//39H INPUT ERROR DETECTED IN (MATRIT/SHELL),//,
     1        19H ELEMENT GROUP NO = ,I5,/,
     2        27H MATERIAL PROPERTY SET NO = ,I5,/,
     3        32H STIFFENING MATERIAL NOT ALLOWED  )
 3450 FORMAT (//39H INPUT ERROR DETECTED IN (MATRIT/SHELL),//,
     1          20H ELEMENT GROUP NO = ,I5,/,
     2          28H MATERIAL PROPERTY SET NO = ,I5,/,
     3          45H ORTHOTROPIC PROPERTIES INCORRECTLY SPECIFIED,/,
     4          17H SEE USERS MANUAL ,//)
C
C
      END
      SUBROUTINE RSTNOD(COSXY,VNI,VNT,V1,ANG,IVCOD)
C
C
C     NODES
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      DIMENSION COSXY(*),VNI(*),VNT(*),V1(*),ANG(*)
C
      DO 2 I=1,6
    2 COSXY(I)=0.D0
C
      VTOL=1.0D-8
C
      IF (IVCOD.EQ.2) GO TO 40
C
C
C
      TEMP=ABS(VNI(2)) - 1.0
      TEMP=ABS(TEMP)
      IF (TEMP.GT.VTOL) GO TO 10
C
      VNT(1)=-ANG(1)*ABS(VNI(2))
      VNT(2)=VNI(2)
      VNT(3)= ANG(2)*VNI(2)
      IF (IVCOD) 60,5,7
    5 DO 6 L=1,3
    6 VNI(L)=VNT(L)
    7 TEMP=1.0 - ABS(VNT(2))
      IF (TEMP.GT.VTOL) GO TO 50
C
      COSXY(3)=VNT(2)
      COSXY(4)=VNT(2)
      IF (VNT(2).LT.0.D0) COSXY(4)=-VNT(2)
      RETURN
C
C     REGULAR CASE
C
   10 DUM=VNI(1)*VNI(1) + VNI(3)*VNI(3)
      DUM=SQRT(DUM)
      COSXY(1)=VNI(3)/DUM
      COSXY(2)=0.D0
      COSXY(3)=-VNI(1)/DUM
      TEMP1=COSXY(3)*VNI(2)
      TEMP2=-COSXY(3)*VNI(1) + COSXY(1)*VNI(3)
      TEMP3=-COSXY(1)*VNI(2)
      DUM=SQRT(TEMP1*TEMP1 + TEMP2*TEMP2 + TEMP3*TEMP3)
      COSXY(4)=TEMP1/DUM
      COSXY(5)=TEMP2/DUM
      COSXY(6)=TEMP3/DUM
C
      VNT(1)=VNI(1) - COSXY(4)*ANG(1) + COSXY(1)*ANG(2)
      VNT(2)=VNI(2) - COSXY(5)*ANG(1)
      VNT(3)=VNI(3) - COSXY(6)*ANG(1) + COSXY(3)*ANG(2)
      TEMP=1.0 - ABS(VNT(2))
      TEMP=ABS(TEMP)
      IF (TEMP.GT.VTOL) GO TO 18
C
      DO 16 I=1,3
   16 VNT(I)=VNI(I)
      RETURN
C
   18 IF(IVCOD) 25,20,25
C
C
   20 DO 22 L=1,3
   22 VNI(L)=VNT(L)
   25 RETURN
C
C
   40 COSXY(1)=V1(1)
      COSXY(2)=V1(2)
      COSXY(3)=V1(3)
C
      V21=VNT(2)*V1(3)-VNT(3)*V1(2)
      V22=VNT(3)*V1(1)-VNT(1)*V1(3)
      V23=VNT(1)*V1(2)-VNT(2)*V1(1)
      DUM=SQRT(V21*V21+V22*V22+V23*V23)
      DUMI=1./DUM
      COSXY(4)=V21*DUMI
      COSXY(5)=V22*DUMI
      COSXY(6)=V23*DUMI
      RETURN
C
   50 DUM=SQRT(VNT(1)*VNT(1) + VNT(3)*VNT(3))
      COSXY(1)=VNT(3)/DUM
      COSXY(2)=0.D0
      COSXY(3)=-VNT(1)/DUM
      TEMP1=COSXY(3)*VNT(2)
      TEMP2=-COSXY(3)*VNT(1) + COSXY(1)*VNT(3)
      TEMP3=-COSXY(1)*VNT(2)
      DUM=SQRT(TEMP1*TEMP1 + TEMP2*TEMP2 + TEMP3*TEMP3)
      COSXY(4)=TEMP1/DUM
      COSXY(5)=TEMP2/DUM
      COSXY(6)=TEMP3/DUM
C
   60 RETURN
C
      END
      SUBROUTINE ROTBTH (ROTB,X,LM,NDOPT,IGLOB,IELD,NEQ,NEQT,KKK)
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
C     THIS SUBROUTINE
C
C
      DIMENSION ROTB(*),X(*),LM(*),NDOPT(*),IGLOB(*)
      DIMENSION RB(3)
C
      IL=0
      IR=0
      NNMS=0
      DO 50 I=1,IELD
      IL=IL+3
      IF (NDOPT(I).GT.0) GO TO 50
      NNMS=NNMS + 1
      NROTN=2
      IF (IGLOB(NNMS).EQ.0) NROTN=3
      DO 22 J=1,NROTN
   22 RB(J)=ROTB(IR+J)
      IR=IR + NROTN
      DO 40 J=1,NROTN
      II=LM(IL+J)
      IF (II.EQ.0 .OR. II.GT.NEQT) GO TO 40
      IF (II.LT.0) II=NEQ-II
C
      IF (KKK.EQ.1) X(II)=X(II) - RB(J)
      IF (KKK.EQ.2) X(II)=X(II) + RB(J)
C
   40 CONTINUE
      IL=IL + NROTN
   50 CONTINUE
C
      RETURN
C
      END
      SUBROUTINE ANCAL (VN,V1,ANG,RSDCOS,X,LM,NEGSKS,ISKEW,IVCOD)
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
C
      DIMENSION VN(*),V1(*),ANG(*),LM(*),RSDCOS(9,*),X(*)
      DIMENSION ROT(3),AUX(3)
C
      COMMON /SOL/ NUMNP,NEQ,NWK,NWM,NWC,NUMEST,MIDEST,MAXEST,NSTE,MA
      COMMON /ISUBST/ ISUBC,NSUBST,NSUB,NTUSE,NEGLS,NEGNLS,NUMNPS,
     1                NODCON,NODRET,IDOFS(12),NDOFS,NEQS,NWKS,MAXESC,
     1                MAS,NSTAPE,ILOA(13),KRSIZM,NEQC
      COMMON /DISCON/ NDISCE,NIDM
C
      NEQT = NEQ + NDISCE
      IF (ISUBC.GT.0) NEQT = NEQS
C
C
      VNX=VN(1)
      VNY=VN(2)
      VNZ=VN(3)
      CALL V1CAL (VNX,VNY,VNZ,V1,V1X,V1Y,V1Z,IVCOD)
      V2X=VNY*V1Z - V1Y*VNZ
      V2Y=VNZ*V1X - V1Z*VNX
      V2Z=VNX*V1Y - V1X*VNY
      DO 10 L=1,3
      II=LM(1+L)
      IF (II.LT.0) II=NEQ-II
      ROT(L)=0.D0
      IF (II.GT.0 .AND. II.LE.NEQT) ROT(L) = X(II)
   10 CONTINUE
      IF (NEGSKS.EQ.0) GO TO 50
      IF (ISKEW.LE.0) GO TO 50
      DO 15 I=1,3
   15 AUX(I)=0.D0
      DO 30 I=1,3
      MM=I
      DO 20 J=1,3
      AUX(I)=AUX(I) + RSDCOS(MM,ISKEW)*ROT(J)
   20 MM=MM + 3
   30 CONTINUE
      DO 40 I=1,3
   40 ROT(I)=AUX(I)
   50 ANG(1)=V1X*ROT(1) + V1Y*ROT(2) + V1Z*ROT(3)
      ANG(2)=V2X*ROT(1) + V2Y*ROT(2) + V2Z*ROT(3)
      RETURN
      END
      SUBROUTINE V1CAL (VNX,VNY,VNZ,V1,V1X,V1Y,V1Z,IVCOD)
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
C
      DIMENSION V1(*)
C
      IF (IVCOD.EQ.2) GO TO 50
C
      VTOL=1.D-08
      VNYA=ABS(VNY)
      TEMP=ABS(VNYA - 1.)
      IF (TEMP.GT.VTOL) GO TO 30
C
C     SPECIAL CASE
C
      IF (VNY.LT.0.D0) GO TO 20
      V1X=0.D0
      V1Y=0.D0
      V1Z=1.D0
      RETURN
   20 V1X=0.D0
      V1Y=0.D0
      V1Z=-1.D0
      RETURN
C
C     STANDARD CASE
C
   30 V1X=VNZ
      V1Y=0.D0
      V1Z=-VNX
      DUM=SQRT(V1X*V1X+V1Z*V1Z)
      DUMI=1./DUM
      V1X=V1X*DUMI
      V1Y=V1Y*DUMI
      V1Z=V1Z*DUMI
      RETURN
   50 V1X=V1(1)
      V1Y=V1(2)
      V1Z=V1(3)
      RETURN
      END
      SUBROUTINE CAUSHL
C
C .                                                                   .
C .                                                                   .
C .                                                                   .
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLTN,KPLOTE
      COMMON /SHELL1/ DISD(9),IELD,IELP,NPT,IDW,NDROT
      COMMON /SHELMT/ D(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
C
C
      F11=DISD(1) + 1.
      F12=DISD(4)
      F13=DISD(5)
      F21=DISD(6)
      F22=DISD(2) + 1.
      F23=DISD(7)
      F31=DISD(8)
      F32=DISD(9)
      F33=DISD(3) + 1.
C
C
C     USE DET(F)=1.0
C
      DET=1.D0
      S11=STRESS(1)
      S22=STRESS(2)
      S33=STRESS(3)
      S12=STRESS(4)
      S13=STRESS(5)
      S23=STRESS(6)
C
      PKFT1=S11*F11 + S12*F12 + S13*F13
      PKFT2=S12*F11 + S22*F12 + S23*F13
      PKFT3=S13*F11 + S23*F12 + S33*F13
      STRESS(1)= DET*(F11*PKFT1 + F12*PKFT2 + F13*PKFT3)
      STRESS(4)= DET*(F21*PKFT1 + F22*PKFT2 + F23*PKFT3)
      STRESS(5)= DET*(F31*PKFT1 + F32*PKFT2 + F33*PKFT3)
C
      PKFT1=S11*F21 + S12*F22 + S13*F23
      PKFT2=S12*F21 + S22*F22 + S23*F23
      PKFT3=S13*F21 + S23*F22 + S33*F23
      STRESS(2)= DET*(F21*PKFT1 + F22*PKFT2 + F23*PKFT3)
C
      PKFT1=S11*F31 + S12*F32 + S13*F33
      PKFT2=S12*F31 + S22*F32 + S23*F33
      PKFT3=S13*F31 + S23*F32 + S33*F33
      STRESS(3)= DET*(F31*PKFT1 + F32*PKFT2 + F33*PKFT3)
      STRESS(6)= DET*(F21*PKFT1 + F22*PKFT2 + F23*PKFT3)
C
      RETURN
      END
      SUBROUTINE SIGROT (STR,NRX,ISTR)
C
C     COORDINATE SYSTEM
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /SHROT/ XJ(3,3),DCA(3,3)
      COMMON /SHELL4/ XJI(3,3),P(3,32),H(32)
      COMMON /SHELL1/ DISD(9),IELD,IELP,NPT,IDW,NDROT
      COMMON /GT/ XJT(3,3),EPS(4)
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1           ,NDOFDM,KLIN,IEIG,IMASSN,IDAMPN
C
      DIMENSION SIG(3,3),DSIG(3,3),DJ(3,3),STR(*)
      DIMENSION DC(3)
C
      IF (IELD.NE.4 .OR. NPAR(3).NE.2) GO TO 01
C
      TNORM=0.D0
      DO 62 J=1,3
   62 TNORM=TNORM + XJT(3,J)*XJT(3,J)
      TNORM=SQRT(TNORM)
      DO 63 J=1,3
   63 DCA(J,3)=XJT(3,J)/TNORM
C
      DC(1)=XJT(2,2)*XJT(3,3) - XJT(2,3)*XJT(3,2)
      DC(2)=XJT(2,3)*XJT(3,1) - XJT(2,1)*XJT(3,3)
      DC(3)=XJT(2,1)*XJT(3,2) - XJT(2,2)*XJT(3,1)
      TNORM=0.D0
      DO 60 J=1,3
   60 TNORM=TNORM + DC(J)*DC(J)
      TNORM=SQRT(TNORM)
      DO 64 J=1,3
   64 DCA(J,1)=DC(J)/TNORM
C
      DCA(1,2)=DCA(2,3)*DCA(3,1) - DCA(2,1)*DCA(3,3)
      DCA(2,2)=DCA(1,1)*DCA(3,3) - DCA(1,3)*DCA(3,1)
      DCA(3,2)=DCA(1,3)*DCA(2,1) - DCA(1,1)*DCA(2,3)
   01 CONTINUE
C
C
      FAC=1.D0
      IF (ISTR.EQ.2) FAC=2.D0
      SIG(1,1)=STR(1)
      SIG(1,2)=STR(4)/FAC
      SIG(1,3)=STR(5)/FAC
      SIG(2,2)=STR(2)
      SIG(2,3)=STR(6)/FAC
      SIG(3,3)=STR(3)
      SIG(2,1)=SIG(1,2)
      SIG(3,1)=SIG(1,3)
      SIG(3,2)=SIG(2,3)
C
      IF (NRX-1) 100,2,5
C
C
    2 DO 3 I=1,3
      DO 3 J=1,3
    3 DJ(I,J)=DCA(I,J)
      GO TO 10
C
C
    5 DO 6 I=1,3
      DO 6 J=1,3
    6 DJ(I,J)=DCA(J,I)
C
C
   10 DO 20 I=1,3
      DO 20 J=1,3
      TEMP=0.D0
      DO 22 L=1,3
   22 TEMP=TEMP + SIG(I,L)*DJ(L,J)
   20 DSIG(I,J)=TEMP
C
      DO 30 I=1,3
      DO 30 J=I,3
      TEMP=0.D0
      DO 32 L=1,3
   32 TEMP=TEMP + DJ(L,I)*DSIG(L,J)
   30 SIG(I,J)=TEMP
C
C
      DO 50 I=1,3
   50 STR(I)=SIG(I,I)
      STR(4)=SIG(1,2)*FAC
      STR(5)=SIG(1,3)*FAC
      STR(6)=SIG(2,3)*FAC
C
  100 RETURN
      END
      SUBROUTINE ZEROWA (MODEL)
C .                                                                   .
C .                                                                   .
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      GO TO (1,4,3,2),MODEL
C
C
C
    1 RETURN
C
C
    2 CALL SHMAT2
      RETURN
C
C
    3 RETURN
C
C
    4 RETURN
      END
      SUBROUTINE SHBASE (NINTR,NINTS,NINTRS)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      COMMON /TRANG/ TRLW4(4,3),TRLW7(7,3),TRLWD(13,3),
     1               XGRS(36,2),WGTRS(36)
      COMMON /SHELL5/ ISHAPE
C
      NINTRS=NINTR*NINTS
C
      IF (ISHAPE) 60,60,70
C
C
   60 K=0
      DO 65 I=1,NINTR
      DO 65 J=1,NINTS
      K=K + 1
      XGRS(K,1)=XG(I,NINTR)
      XGRS(K,2)=XG(J,NINTS)
   65 WGTRS(K)=WGT(I,NINTR)*WGT(J,NINTS)
      RETURN
C
C
   70 IF (NINTRS.GT.1) GO TO 75
      XGRS(1,1)=-1.0D0/3.0D0
      XGRS(1,2)=0.D0
      WGTRS(1)=2.0D0