a04b.for

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

   80 CONTINUE
C
      DO 85 I=1,NDOFL
   85 DUM(I)=EDIS(I)
      DO 90 I=1,IEL
      I2=2*I
      I3=3*I
      EDIS(I3-2)=T(1,1)*DUM(I2-1) + T(2,1)*DUM(I2)
      EDIS(I3-1)=T(1,2)*DUM(I2-1) + T(2,2)*DUM(I2)
      EDIS(I3  )=T(1,3)*DUM(I2-1) + T(2,3)*DUM(I2)
   90 CONTINUE
C
      RETURN
C
C
  100 NDIM3=NDOFG*(NDOFG+1)/2
      DO 102 I=1,NDIM3
  102 DUM(I)=0.D0
      I=1
      KG1=1
      KG2=1
      KG3=2
      DO 120 L=1,IEL
      IP=I+1
      I1=I-1
      J =(I1+I1)/3 + 1
      J1=J-1
      ND2=NDOFL*J1 - (J1-1)*J1/2 + 1
      ND3=NDOFG*I1 - (I1-1)*I1/2 + 1
      NDG=NDOFG*IP - IP*I/2 + 1
      LD2=0
      LD3=0
      LDG=0
      DO 130 K=L,IEL
      DUM(ND3+LD3)=S(ND2+LD2) + SGNL(KG1)
      DUM(ND3+LD3+1)=S(ND2+LD2+1)
      DUM(NDG+LDG)=SGNL(KG1)
      KG1=KG1+1
      LDG=LDG+3
      LD2=LD2+2
  130 LD3=LD3+3
      L1=L-1
      LP=L+1
      MD3=ND3 + NDOFG - 3*L1
      MD2=ND2 + NDOFL - L1 - L1
      DUM(MD3)=S(MD2) + SGNL(KG2)
      IF (LP.GT.IEL) GO TO 120
      KG2=KG2 + IEL - L1
      LD2=1
      LD3=2
      DO 140 K=LP,IEL
      DUM(MD3+LD3)=S(MD2+LD2)
      DUM(MD3+LD3+1)=S(MD2+LD2+1) + SGNL(KG3)
      KG3=KG3+1
      LD2=LD2+2
  140 LD3=LD3+3
      KG3=KG3+1
  120 I=I+3
C
C
      IR=0
      DO 160 K=1,IEL
  160 ILSK(K)=1
      DO 180 I=1,3
      DO 182 J=1,3
  182 TMA(J+IR)=T(J,I)
  180 IR=IR+3
C
      CALL ATKA (TMA,DUM,ILSK,IEL,3)
C
      DO 190 I=1,NDIM3
  190 S(I)=DUM(I)
C
C
      RETURN
C
      END
      SUBROUTINE QUADS (ND,B,S,YZ,PROP,RE,EDIS,EDISI,IDW,WA,NOD5)
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 /TODIM/ BET,THIC,DE,IEL,NND5,ISOCOR
      COMMON /MATMOD/ STRESS(4),STRAIN(4),D(4,4),IPT,NEL,IPS,ISVE
      COMMON /DISDER/ DISD(5)
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      COMMON /GNLSTF/ SGNL(45)
C
      DIMENSION B(4,*),S(*),YZ(*),RE(*),EDIS(*),PROP(*),WA(*),NOD5(*)
      DIMENSION DB(4),XX(18),BS(4,18),DI(4,4),EDISI(*)
C
      EQUIVALENCE (NPAR(10),NINT),(NPAR(5),ITYP2D),(NPAR(3),INDNL)
      EQUIVALENCE (NPAR(15),MODEL)
C
C
      NPT=NINT*NINT
      IST=4
      IF (ITYP2D.NE.0) IST=3
      KST=IST-1
C
      IF (IND.GE.4) GO TO 100
C
C
C
C
      CALL STSTL (NEL,YZ,PROP,D)
C
      DO 10 LX=1,NINT
      E1=XG(LX,NINT)
      DO 10 LY=1,NINT
      E2=XG(LY,NINT)
      WT=WGT(LX,NINT)*WGT(LY,NINT)
C
C
      CALL DERIQ (NEL,YZ,B,DET,E1,E2,XBAR,NOD5)
C
C
      IF (IST.EQ.3) XBAR=THIC
      FAC=WT*XBAR*DET
C
      KL=1
      DO 50 J=1,ND,2
      DO 52 K=1,3
      DB(K)=D(K,1)*B(1,J) + D(K,3)*B(3,J)
   52 DB(K)=DB(K)*FAC
      DO 51 I=J,ND,2
      S(KL)=S(KL) + B(1,I)*DB(1) + B(3,I)*DB(3)
      KL=KL + 1
      S(KL)=S(KL) + B(2,I+1)*DB(2) + B(3,I+1)*DB(3)
   51 KL=KL + 1
   50 KL=KL + ND - J
C
      KL=ND + 1
      DO 54 J=2,ND,2
      DO 56 K=1,3
      DB(K)=D(K,2)*B(2,J) + D(K,3)*B(3,J)
   56 DB(K)=DB(K)*FAC
      KS=KL
      DO 55 I=J,ND,2
      S(KS)=S(KS) + B(2,I)*DB(2) + B(3,I)*DB(3)
   55 KS=KS + 2
      IF (J-ND) 57,54,54
   57 K=J + 1
      KS=KL + 1
      DO 58 II=K,ND,2
      S(KS)=S(KS) + B(1,II)*DB(1) + B(3,II)*DB(3)
   58 KS=KS + 2
   54 KL=KL + 2*ND - 2*J + 1
C
      IF (IST.EQ.3) GO TO 10
      KL=1
      DO 60 J=1,ND,2
      DB(1)=D(1,4)*B(4,J)*FAC
      DB(2)=D(2,4)*B(4,J)*FAC
      DB(3)=D(3,4)*B(4,J)*FAC
      DB(4)=D(4,1)*B(1,J) + D(4,3)*B(3,J) + D(4,4)*B(4,J)
      DB(4)=DB(4)*FAC
      DO 61 I=J,ND,2
      S(KL)=S(KL) + B(1,I)*DB(1) + B(3,I)*DB(3) + B(4,I)*DB(4)
      KL=KL + 1
      S(KL)=S(KL) + B(2,I+1)*DB(2) + B(3,I+1)*DB(3)
   61 KL=KL + 1
   60 KL=KL + ND - J
      KL=ND + 1
      DO 59 J=2,ND,2
      DB(4)=D(4,2)*B(2,J) + D(4,3)*B(3,J)
      DB(4)=DB(4)*FAC
      DO 62 I=J,ND
      S(KL)=S(KL) + B(4,I)*DB(4)
   62 KL=KL + 1
   59 KL=KL + ND - J
C
   10 CONTINUE
C
      RETURN
C
C
C
C
C
  100 IF (INDNL.LE.2) GO TO 122
      IF (ITYP2D.LE.2) GO TO 118
      DO 116 I=1,ND
  116 XX(I)=YZ(I)
      GO TO 122
  118 DO 120 J=1,ND
  120 XX(J) = YZ(J) + EDIS(J)
C
C
  122 IF (MODEL.GT.2) GO TO 125
      IF (INDNL.LE.2) GO TO 124
      CALL STSTL (NEL,XX,PROP,D)
      GO TO 125
  124 CALL STSTL (NEL,YZ,PROP,D)
C
C
C
C
  125 DO 300 LX=1,NINT
      E1=XG(LX,NINT)
      DO 300 LY=1,NINT
      E2=XG(LY,NINT)
      WT=WGT(LX,NINT)*WGT(LY,NINT)
      IPT=(LX-1)*NINT + LY
      IF (INDNL.EQ.3) GO TO 200
C
C
C
C
C
      CALL DERIQ (NEL,YZ,B,DET,E1,E2,XBAR,NOD5)
C
C
      DO 130 I=1,5
  130 DISD(I)=0.D0
      DO 140 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)
  140 DISD(4)=DISD(4) + B(3,J)*EDIS(J)
      IF (IST.EQ.3) GO TO 160
      DO 150 I=1,ND,2
  150 DISD(5)=DISD(5) + B(4,I)*EDIS(I)
C
C
  160 CALL STSTN (YZ,PROP,DISD,IDW,WA)
C
      IF (INDNL.LE.1) GO TO 221
C
C
      DO 164 J=2,ND,2
      I=J - 1
      BS(1,I)=B(1,I) + B(1,I)*DISD(1)
      BS(1,J)=B(1,I)*DISD(4)
      BS(2,I)=B(2,J)*DISD(3)
      BS(2,J)=B(2,J) + B(2,J)*DISD(2)
      BS(3,I)=B(3,I) + B(3,I)*DISD(1) + B(3,J)*DISD(3)
  164 BS(3,J)=B(3,J) + B(3,I)*DISD(4) + B(3,J)*DISD(2)
      IF (IST.EQ.3) GO TO 167
      DO 166 I=1,ND,2
      J=I + 1
      BS(4,J)=0.D0
  166 BS(4,I)=B(4,I) + B(4,I)*DISD(5)
C
C
  167 IF (IST.EQ.3) XBAR=THIC
      FAC=WT*XBAR*DET
      TAU11=STRESS(1)*FAC
      TAU22=STRESS(2)*FAC
      TAU12=STRESS(3)*FAC
      TAU33=STRESS(4)*FAC
      DO 170 I=1,ND
  170 RE(I)=RE(I) + BS(1,I)*TAU11 + BS(2,I)*TAU22 + BS(3,I)*TAU12
      IF (IST.EQ.3) GO TO 176
      DO 174 J=1,ND,2
  174 RE(J)=RE(J) + BS(4,J)*TAU33
C
  176 IF (ICOUNT - 2) 178,178,300
  178 IF (IREF) 300,179,300
C
C
  179 DO 183 I=1,IST
      DO 183 J=I,IST
      DI(I,J)=D(I,J)*FAC
  183 DI(J,I)=DI(I,J)
      KL=0
      DO 180 J=1,ND
      DO 182 K=1,IST
      DB(K)=0.D0
      DO 184 L=1,IST
  184 DB(K)=DB(K) + DI(K,L)*BS(L,J)
  182 CONTINUE
C
      DO 180 I=J,ND
      KL=KL + 1
      DUM=0.D0
      DO 186 K=1,IST
  186 DUM=DUM + BS(K,I)*DB(K)
  180 S(KL)=S(KL) + DUM
C
      GO TO 365
C
C
C
C
C
  200 CALL DERIQ (NEL,XX,B,DET,E1,E2,XBAR,NOD5)
C
C
C
      IF (MODEL.GT.1) GO TO 215
      DO 210 I=1,5
  210 DISD(I)=0.D0
      DO 212 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)
  212 DISD(4)=DISD(4) + B(3,J)*EDIS(J)
      IF (IST.EQ.3) GO TO 216
      DO 214 I=1,ND,2
  214 DISD(5)=DISD(5) + B(4,I)*EDIS(I)
      GO TO 216
C
C
  215 DO 217 I=1,5
  217 DISD(I)=0.D0
      DO 218 J=2,ND,2
      I=J - 1
      DISD(1)=DISD(1) + B(1,I)*EDISI(I)
      DISD(2)=DISD(2) + B(2,J)*EDISI(J)
      DISD(3)=DISD(3) + B(3,I)*EDISI(I)
  218 DISD(4)=DISD(4) + B(3,J)*EDISI(J)
      IF (IST.EQ.3) GO TO 216
      DO 219 I=1,ND,2
  219 DISD(5)=DISD(5) + B(4,I)*EDISI(I)
C
C
  216 CALL STSTN (XX,PROP,DISD,IDW,WA)
C
  221 IF (ITYP2D.EQ.0) GO TO 222
      XBAR=THIC
      IF (INDNL.LE.1 .OR. ITYP2D.EQ.1) GO TO 222
      IF (MODEL.GT.1) GO TO 223
      EXT=1.0 - 2.0*STRAIN(4)
      XBAR=XBAR/SQRT(EXT)
      GO TO 222
C
  223 XBAR=THIC*EXP(STRAIN(4))
C
  222 FAC=WT*XBAR*DET
C
C
      TAU11=STRESS(1)*FAC
      TAU22=STRESS(2)*FAC
      TAU12=STRESS(3)*FAC
      TAU33=STRESS(4)*FAC
      DO 340 J=2,ND,2
      I=J - 1
      RE(I)=RE(I) + B(1,I)*TAU11 + B(3,I)*TAU12
  340 RE(J)=RE(J) + B(2,J)*TAU22 + B(3,J)*TAU12
      IF (IST.EQ.3) GO TO 350
      DO 345 J=1,ND,2
  345 RE(J)=RE(J) + B(4,J)*TAU33
C
  350 IF (ICOUNT-2) 220,220,300
  220 IF (IREF) 300,230,300
C
C
  230 DO 232 I=1,IST
      DO 232 J=I,IST
      DI(I,J)=D(I,J)*FAC
  232 DI(J,I)=DI(I,J)
      KL=1
      DO 250 J=1,ND,2
      DO 252 K=1,3
  252 DB(K)=DI(K,1)*B(1,J) + DI(K,3)*B(3,J)
      DO 251 I=J,ND,2
      S(KL)=S(KL) + B(1,I)*DB(1) + B(3,I)*DB(3)
      KL=KL + 1
      S(KL)=S(KL) + B(2,I+1)*DB(2) + B(3,I+1)*DB(3)
  251 KL=KL + 1
  250 KL=KL + ND - J
      KL=ND + 1
C
      DO 254 J=2,ND,2
      DO 256 K=1,3
  256 DB(K)=DI(K,2)*B(2,J) + DI(K,3)*B(3,J)
      KS=KL
      DO 255 I=J,ND,2
      S(KS)=S(KS) + B(2,I)*DB(2) + B(3,I)*DB(3)
  255 KS=KS + 2
      IF (J-ND) 257,254,254
  257 K=J + 1
      KS=KL + 1
      DO 258 II=K,ND,2
      S(KS)=S(KS) + B(1,II)*DB(1) + B(3,II)*DB(3)
  258 KS=KS + 2
  254 KL=KL + 2*ND - 2*J + 1
C
      IF (IST.EQ.3) GO TO 365
      KL=1
      DO 260 J=1,ND,2
      DB(1)=DI(1,4)*B(4,J)
      DB(2)=DI(2,4)*B(4,J)
      DB(3)=DI(3,4)*B(4,J)
      DB(4)=DI(4,1)*B(1,J) + DI(4,3)*B(3,J) + DI(4,4)*B(4,J)
      DO 261 I=J,ND,2
      S(KL)=S(KL) + B(1,I)*DB(1) + B(3,I)*DB(3) + B(4,I)*DB(4)
      KL=KL + 1
      S(KL)=S(KL) + B(2,I+1)*DB(2) + B(3,I+1)*DB(3)
  261 KL=KL + 1
  260 KL=KL + ND - J
      KL=ND + 1
      DO 259 J=2,ND,2
      DB(4)=DI(4,2)*B(2,J) + DI(4,3)*B(3,J)
      DO 262 I=J,ND
      S(KL)=S(KL) + B(4,I)*DB(4)
  262 KL=KL + 1
  259 KL=KL + ND - J
C
C
C
C
C
  365 IF (INDNL.EQ.1) GO TO 300
      IF (ITYP2D.EQ.3) GO TO 500
C
      KL=1
      DO 400 J=1,ND,2
      DB1=TAU11*B(1,J) + TAU12*B(3,J)
      DB2=TAU12*B(1,J) + TAU22*B(3,J)
C
      KS=KL
      DO 401 I=J,ND,2
      KSS=KS + ND - J + 1
      DUM=B(1,I)*DB1 + B(3,I)*DB2
      S(KS)=S(KS) + DUM
      S(KSS)=S(KSS) + DUM
  401 KS=KS + 2
  400 KL=KL + 2*ND - 2*J + 1
C
      IF (IST.EQ.3) GO TO 300
      KL=1
      DO 420 J=1,ND,2
      DB3=TAU33*B(4,J)
      DO 421 I=J,ND,2
      S(KL)=S(KL) + DB3*B(4,I)
  421 KL=KL + 2
  420 KL=KL + ND - J
      GO TO 300
C
  500 KS=1
      DO 510 J=1,ND,2
      DB1=TAU11*B(1,J) + TAU12*B(3,J)
      DB2=TAU12*B(1,J) + TAU22*B(3,J)
      DO 512 I=J,ND,2
      SGNL(KS)=SGNL(KS) + B(1,I)*DB1 + B(3,I)*DB2
  512 KS=KS+1
  510 CONTINUE
C
  300 CONTINUE
C
C
      RETURN
      END
      SUBROUTINE QUADM (NEL,ND,XM,CM,XX,NOD5)
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 /TODIM/ BET,THIC,DE,IEL,NND5,ISOCOR
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      DIMENSION CM(*),XM(27),D(18),XX(2,9),NOD5(*)
      DIMENSION H(9),P(2,9),XJ(2,2)
C
      EQUIVALENCE (NPAR(5),ITYP2D)
C
C
C
      IINTP=0
      IF (IMASS.EQ.1) GO TO 9
      DO 08 I=1,378
    8 CM(I)=0.D0
    9 DO 7 I=1,ND
    7 XM(I)=0.D0
C
      DO 100 LX=1,3
      R=XG(LX,3)