a04b.for

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

      SUBROUTINE INITWA (MODEL)
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      GO TO (1,2,3,4,4,6,7,8,8,10,10,12,13,14,14,14),MODEL
C
C
C
    1 RETURN
C
C
C
    2 RETURN
C
C
C
    3 CALL ELT2D3
      RETURN
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
   13 RETURN
C
C
   14 CALL EL2D14
      RETURN
C
C
C
      END
      SUBROUTINE MATRT2 (N,DEN,PROP)
C
C
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 /PRCON/ IDATWR,IPRIC,NPB,IDC,IVC,IAC,IPC,IPNODE(3,15)
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      DIMENSION PROP(*)
C
      EQUIVALENCE (NPAR(15),MODEL),(NPAR(17),NCON),(NPAR(20),IDW)
     1           ,(NPAR(16),NUMMAT),(NPAR(19),ITHERM)
C
C
      IF(IDATWR.LE.1) GO TO 500
      IF (MODEL.EQ.3 .OR. MODEL.EQ.5) GO TO 600
      IF (MODEL.EQ.7 .OR. MODEL.EQ.8 .OR. MODEL.EQ.9) GO TO 600
      IF (MODEL.EQ.10 .OR. MODEL.EQ.11) GO TO 600
      IF (MODEL.EQ.14) GO TO 600
      IF (IDATWR.GT.1 .AND. MODEL.EQ.1) GO TO 15
      IF (IDATWR.GT.1 .AND. MODEL.EQ.2) GO TO 25
C
      IF (IDATWR.GT.1 .AND. MODEL.EQ.4) GO TO 21
      RETURN
C
  500 write(66,2100) N,DEN
C
  600 GO TO (1,2,3,4,5,6,7,8,9,10,10,12,13,14,14,14),MODEL
C
C
C
    1 write(66,2101) (PROP(I), I=1,NCON)
   15 IBUG=0
C
      IF ( PROP(1).LE.0.D0) IBUG=1
      IF (PROP(2).LE.-1.D0 .OR. PROP(2).GE.0.5D0) IBUG=1
      IF (IBUG.EQ.0) RETURN
      write(66,3600) NG,N
      write(66,3650)
      write(66,3700)
      IF ( MODEX.EQ.0 ) RETURN
      STOP
C
C
C
    2 write(66,2102) (PROP(I), I=1,NCON)
   25 IBUG=0
      JBUG=0
      KBUG=0
C
C
      DO 35 I=1,3
   35 IF (PROP(I).LE.0.D0) IBUG=1
      IF (PROP(7).LE.0.D0) IBUG=1
      IF (IBUG.EQ.1) GO TO 36
C
C
      IF (ABS(PROP(4)).GE.(SQRT(ABS(PROP(2)/PROP(1))))) JBUG=1
      IF (ABS(PROP(5)).GE.(SQRT(ABS(PROP(3)/PROP(1))))) JBUG=1
      IF (ABS(PROP(6)).GE.(SQRT(ABS(PROP(3)/PROP(2))))) JBUG=1
C
C
      CHECK1=0.5*(1.-PROP(4)*PROP(4)*PROP(1)/PROP(2)
     1              -PROP(6)*PROP(6)*PROP(2)/PROP(3)
     2              -PROP(5)*PROP(5)*PROP(1)/PROP(3))
      CHECK2=PROP(4)*PROP(5)*PROP(6)*PROP(1)/PROP(3)
C
      IF (CHECK1.LE.CHECK2 .OR. CHECK1.GT.0.5D0) KBUG=1
C
   36 IF ( IBUG.EQ.0 .AND. JBUG.EQ.0 .AND. KBUG.EQ.0 ) RETURN
C
      write(66,3600) NG,N
      write(66,3750)
      write(66,3700)
      IF ( MODEX.EQ.0 ) RETURN
      STOP
C
C
C
    3 IBUG=0
C
      IF (IDATWR.GT.1) GO TO 78
      write(66,2103)
      DO 80 K=1,16
      IP1=K + 16
      IP2=K + 32
      IP3=K + 48
   80 write(66,2104) PROP(K),PROP(IP1),PROP(IP2),PROP(IP3)
      write(66,2105) PROP(65),PROP(66)
   78 CONTINUE
C
      NPTS=INT(PROP(65))
      IF(NPTS.GT.0) GO TO 60
      PROP(65)=1.6D1
      NPTS=16
      GO TO 72
C
   60 IF(NPTS.GE.2 .AND. NPTS.LE.16) GO TO 72
      IBUG=1
      write(66,3600) NG,N
      write(66,3002)
      write(66,3700)
      GO TO 85
C
   72 DO 75 J=2,NPTS
      JJ=J-1
      IF(PROP(J).GT.PROP(JJ)) GO TO 75
      IBUG=1
      write(66,3600) NG,N
      write(66,3003)
      write(66,3700)
      GO TO 85
   75 CONTINUE
C
C
      DO 86 K=1,NPTS
      IP1=K+16
      IP2=K+32
      IF (PROP(IP1).LE.0.D0) IBUG=1
   86 IF (PROP(IP2).LE.-1.D0 .OR. PROP(IP2).GE.0.5D0) IBUG=1
C
      IF ( IBUG.EQ.0 ) RETURN
      write(66,3600) NG,N
      write(66,3850)
      write(66,3650)
      write(66,3700)
   85 IF ( MODEX.EQ.0 ) RETURN
      STOP
C
C
C
    4 ICRACK=INT(PROP(25))
      write(66,2220) ICRACK,(PROP(I),I=26,NCON)
      IP=NCON/4 - 1
      write(66,2200)
      DO 20 I=1,IP
      IPI=I + IP
      IPI2=IPI + IP
      IPI3=IPI2 + IP
   20 write(66,2210) I,PROP(I),PROP(IPI),PROP(IPI2),PROP(IPI3)
      write(66,2211)
   21 IBUG=0
      JBUG=0
      DO 22 I=7,12
      IK1=I+6
      IK2=I+12
   22 IF (PROP(I).LE.0.D0 .OR. PROP(IK1).LE.0.D0 .OR. PROP(IK2).LE.0.D0)
     1   IBUG=1
      IF (PROP(27).LT.0.D0 .OR. PROP(28).LT.0.D0) JBUG=1
      IF (PROP(27).GT.1.D0 .OR. PROP(28).GT.1.D0) JBUG=1
C
      IF (IBUG.EQ.0 .AND. JBUG.EQ.0) RETURN
C
      write(66,3600) NG,N
      write(66,3800)
      write(66,3700)
      IF (MODEX.EQ.0) RETURN
      STOP
C
C
C
    5 IF (PROP(34).EQ.0.D0) PROP(34)=1.D0
      IF (PROP(35).EQ.0.D0) PROP(35)=0.7D0
      IF (PROP(37).EQ.0.D0) PROP(37)=0.1D-3
      IF (PROP(38).EQ.0.D0) PROP(38)=0.5D0
      JWARN=0
      JBUG=0
C
C
       SP1= PROP(9)**2 + PROP(10)**2 + PROP(11)**2 + PROP(12)**2 +
     1     PROP(13)**2 + PROP(14)**2
      SP31=PROP(15)**2 + PROP(16)**2 + PROP(17)**2 + PROP(18)**2 +
     2     PROP(19)**2 + PROP(20)**2
      SP32=PROP(21)**2 + PROP(22)**2 + PROP(23)**2 + PROP(24)**2 +
     3     PROP(25)**2 + PROP(26)**2
      SP33=PROP(27)**2 + PROP(28)**2 + PROP(29)**2 + PROP(30)**2 +
     4     PROP(31)**2 + PROP(32)**2
      SPTOL=SQRT(SP1 + SP31 + SP32 + SP33)
C
      IF (SPTOL.LT.1.D-6) JWARN=3
      IF (SPTOL.GT.3.D1) JWARN=2
C
      IF ( PROP(9).GT.PROP(10) .OR. PROP(10).GT.PROP(11) .OR.
     1    PROP(11).GT.PROP(12) .OR. PROP(12).GT.PROP(13) .OR.
     2    PROP(13).GT.PROP(14)) JBUG=1
      IF (PROP(15).GT.PROP(16) .OR. PROP(16).GT.PROP(17) .OR.
     1    PROP(17).GT.PROP(18) .OR. PROP(18).GT.PROP(19) .OR.
     2    PROP(19).GT.PROP(20)) JBUG=1
      IF (PROP(21).GT.PROP(22) .OR. PROP(22).GT.PROP(23) .OR.
     1    PROP(23).GT.PROP(24) .OR. PROP(24).GT.PROP(25) .OR.
     2    PROP(25).GT.PROP(26)) JBUG=1
      IF (PROP(27).GT.PROP(28) .OR. PROP(28).GT.PROP(29) .OR.
     1    PROP(29).GT.PROP(30) .OR. PROP(30).GT.PROP(31) .OR.
     2    PROP(31).GT.PROP(32)) JBUG=1
C
      IF (JWARN.LT.3) GO TO 55
C
      PROP(9)=0.D0
      PROP(10)=0.25D0
      PROP(11)=0.5D0
      PROP(12)=0.75D0
      PROP(13)=1.D0
      PROP(14)=1.2D0
      PROP(15)=1.D0
      PROP(16)=1.4D0
      PROP(17)=1.7D0
      PROP(18)=2.2D0
      PROP(19)=2.5D0
      PROP(20)=2.8D0
      PROP(21)=1.3D0
      PROP(22)=1.5D0
      PROP(23)=2.0D0
      PROP(24)=2.3D0
      PROP(25)=2.7D0
      PROP(26)=3.2D0
      PROP(27)=1.25D0
      PROP(28)=1.45D0
      PROP(29)=1.95D0
      PROP(30)=2.25D0
      PROP(31)=2.65D0
      PROP(32)=3.15D0
   55 IF (IDATWR.GT.1) GO TO 56
C
      write(66,2230) (PROP(I),I=1,8)
      IP1=8
      write(66,2235) (PROP(IP1 + J),J=1,24)
      IF (JBUG.EQ.1) write(66,2244)
      IF (JWARN.EQ.2) write(66,2242)
      IF (JWARN.EQ.3) write(66,2243)
      write(66,2240) (PROP(J),J=33,38)
      write(66,2241) PROP(39)
      write(66,2211)
C
C
   56 IBUG=0
      IF (PROP(1).LE.0.D0) IBUG=1
      IF (PROP(2).LE.-1.D0 .OR. PROP(2).GE.0.5D0) IBUG=1
      IF (PROP(4).LT.0.D0) IBUG=1
      IF (PROP(5).GE.0.D0) IBUG=1
      IF (PROP(6).GE.0.D0) IBUG=1
      IF (PROP(7).GE.0.D0) IBUG=1
      IF (PROP(8).GE.0.D0) IBUG=1
C
      ESEC=PROP(5)/PROP(6)
      IF (ESEC.GT.PROP(1)) IBUG=1
      IF (PROP(7).LT.PROP(5)) IBUG=1
      IF (PROP(8).GT.PROP(6)) IBUG=1
C
      IF (PROP(33).LT.0.D0 .OR. PROP(33).GT.1.D0) IBUG=1
      IF (PROP(34).LT.0.D0) IBUG=1
      IF (PROP(35).LT.0.D0 .OR. PROP(35).GT.1.D0) IBUG=1
C
      IF (PROP(37).LT.0.D0 .OR. PROP(37).GT.1.D0) IBUG=1
      IF (PROP(38).LT.0.D0 .OR. PROP(38).GT.1.D0) IBUG=1
      IF (IBUG.EQ.0 .AND. JBUG.EQ.0) GO TO 54
      IF (IBUG.EQ.0) GO TO 57
      write(66,3600) NG,N
      write(66,3500)
      write(66,3700)
   57 IF (MODEX.EQ.0) GO TO 54
      STOP
   54 CONTINUE
C
      IF (PROP(34).GT.6.D0) JWARN=1
      IF (PROP(36).GT.0.D0) JWARN=1
C
      IF (JWARN.EQ.1) write(66,3501)
      RETURN
C
C
    6 RETURN
C
C
C
    7 IBUG=0
      IF (PROP(1).GT.0.D0 .AND. PROP(2).GE.0.D0) GO TO 142
      IBUG=1
  142 IF (PROP(3).GE.0.D0 .AND. PROP(4).GT.0.D0) GO TO 143
      IBUG=1
  143 IF (PROP(5).LT.0.D0 .AND. PROP(6).LT.0.D0) GO TO 144
      IBUG=1
  144 IF (PROP(7).GE.0.D0) GO TO 145
      IBUG=1
  145 IF (PROP(8).LE.0.D0) GO TO 141
      IBUG=1
  141 IF (IDATWR.LE.1) write(66,2110) (PROP(I),I=1,NCON)
      IF (IBUG.EQ.0) RETURN
      write(66,3600) NG,N
      write(66,3410)
      write(66,3700)
      IF (MODEX.EQ.0) RETURN
      STOP
C
C
C
    8 IF (NCON.GT.4) GO TO 200
C
      IBUG=0
      IF (PROP(1).LE.0.D0) IBUG=1
      IF (PROP(2).LE.-1.D0 .OR. PROP(2).GE.0.5D0) IBUG=1
      IF (PROP(4).LT.0.D0) IBUG=1
      IF (PROP(3).GT.0.D0) GO TO 150
      IBUG=1
  150 IF (PROP(4).LT.PROP(1)) GO TO 152
      IBUG=1
  152 CONTINUE
C
      IF (IDATWR.LE.1) write(66,2106) (PROP(I),I=1,NCON)
      IF (IBUG.EQ.0) RETURN
      write(66,3600) NG,N
      write(66,3401)
      write(66,3700)
      IF (MODEX.EQ.0) RETURN
      STOP
C
  200 IBUG=0
      EPSY=0.01
      IF (PROP(1).GT.0.0) EPSY = PROP(3)/PROP(1)
      IF (PROP(4).EQ.0.0) PROP(4) = EPSY
      IF ((PROP(4) - EPSY) .LT. EPSY*1.D-9) GO TO 163
      IBUG=1
      write(66,3600) NG,N
      write(66,3399)
      write(66,3700)
  163 IF (IDATWR.GT.1) GO TO 160
      write(66,2111) (PROP(I),I=1,3)
      write(66,2112) PROP(3),PROP(4)
C
  160 CONTINUE
      IF (PROP(1).LE.0.D0) IBUG=1
      IF (PROP(2).LE.-1.D0 .OR. PROP(2).GE.0.5D0) IBUG=1
      IF (IBUG.EQ.1) GO TO 162
      IF (PROP(3).GT.0.D0) GO TO 161
      IBUG=1
  162 write(66,3600) NG,N
      write(66,3401)
      write(66,3700)
  161 ICP=4
      DO 165 I=1,6
      IF (PROP(ICP).EQ.0.D0) GO TO 165
      ICP2=ICP+2
      IF (PROP(ICP).NE.PROP(ICP2)) GO TO 165
      IBUG=1
      IF (IDATWR.LE.1) write(66,2114) (PROP(K),K=5,ICP2)
      write(66,3600) NG,N
      write(66,3404) ICP,ICP2
      write(66,3700)
  165 ICP=ICP+2
C
      IF (IBUG.EQ.0) GO TO 167
      IF (MODEX.EQ.0) RETURN
      STOP
C
  167 ETOLD=PROP(1)
      DO 210 J=6,NCON,2
      ET=(PROP(J - 1) - PROP(J - 3))/(PROP(J) - PROP(J - 2))
      IF (IDATWR.LE.1) write(66,2113) PROP(J-1),PROP(J),ET
C
      IF (ET.LT.0.D0) GO TO 169
      IF (ET.GE.PROP(1)) GO TO 169
      IF (ET.LE.ETOLD) GO TO 168
  169 write(66,3600) NG,N
      write(66,3405)
      write(66,3700)
      IF (MODEX.EQ.0) RETURN
      STOP
  168 ETOLD=ET
C
  210 CONTINUE
      write(66,2211)
      RETURN
C
C
C
    9 IF (NCON.GT.4) GO TO 220
C
      IBUG=0
      IF (PROP(1).LE.0.D0) IBUG=1
      IF (PROP(2).LE.-1.D0 .OR. PROP(2).GE.0.5D0) IBUG=1
      IF (PROP(4).LT.0.D0) IBUG=1
      IF (PROP(3).GT.0.D0) GO TO 154
      IBUG=1
  154 IF (PROP(4).LT.PROP(1)) GO TO 156
      IBUG=1
  156 CONTINUE
C
      IF (IDATWR.LE.1) write(66,2106) (PROP(I),I=1,NCON)
      IF (IBUG.EQ.0) RETURN
      write(66,3600) NG,N
      write(66,3401)
      write(66,3700)
      IF (MODEX.EQ.0) RETURN
      STOP
C
  220 IF (IDATWR.GT.1) GO TO 170
      write(66,2111) (PROP(I),I=1,3)
      write(66,2112) PROP(3),PROP(4)
C
  170 IBUG=0
      IF (PROP(1).LE.0.D0) IBUG=1
      IF (PROP(2).LE.-1.D0 .OR. PROP(2).GE.0.5D0) IBUG=1
      IF (IBUG.EQ.1) GO TO 172
      IF (PROP(3).GT.0.D0) GO TO 171
      IBUG=1
  172 write(66,3600) NG,N
      write(66,3401)
      write(66,3700)
  171 ICP=4
      DO 175 I=1,6
      IF (PROP(ICP).EQ.0.D0) GO TO 175
      ICP2=ICP+2
      IF (PROP(ICP).NE.PROP(ICP2)) GO TO 175
      IBUG=1
      IF (IDATWR.LE.1) write(66,2114) (PROP(K),K=5,ICP2)
      write(66,3600) NG,N