a04b.for

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

      DO 100 LY=1,3
      S=XG(LY,3)
      WT=WGT(LX,3)*WGT(LY,3)
C
C
C
      CALL FUNCT2 (R,S,H,P,NOD5,XJ,DET,XX,NEL,IINTP)
C
C
      IF (ITYP2D.EQ.0) GO TO 40
      IF (ITYP2D.GT.0) XBAR=THIC
      GO TO 60
   40 XBAR=0.D0
      DO 50 K=1,IEL
   50 XBAR=XBAR + H(K)*XX(1,K)
C
   60 FAC=WT*XBAR*DET*DE
C
C     CONSISTENT MASS
C
      NDPN=2
      IF(ITYP2D.EQ.3)NDPN=3
      IF (IMASS.LT.2) GO TO 320
      DO 200 I = 1,IEL
      D(2*I-1) = H(I)
  200 D(2*I) = H(I)
      KL=1
      ND1=2*IEL
      NDPN2=NDPN-2
      NDPN1=NDPN-1
      DO 300 I=1,ND1,2
      DO 301 J=I,ND1,2
      CM(KL)=CM(KL) + D(I)*D(J)*FAC
  301 KL=KL + NDPN
  300 KL=KL+(ND-(I+((I-1)/2)*NDPN2))*NDPN1-NDPN2
      GO TO 100
C
C     LUMPED MASS
C
  320 FACM=FAC/IEL
      DO 325 I=1,ND,NDPN
  325 XM(I)=XM(I) + FACM
C
  100 CONTINUE
C
      IF (IMASS.EQ.1) GO TO 335
      IF(NDPN.EQ.3)GO TO 410
C
      KL=1
      DO 401 I=1,ND,2
      KS=KL + ND - I + 1
      DO 400 J=I,ND,2
      CM(KS)=CM(KL)
      KS=KS + 2
  400 KL=KL + 2
  401 KL=KL + ND - I
C
      RETURN
C
C
C
  410 KL=1
      DO 451 I=1,ND,3
      KS1=KL+ND-I+1
      KS2=KS1+ND-I
      DO 450 J=I,ND,3
      CM(KS1)=CM(KL)
      CM(KS2)=CM(KL)
      KL=KL+3
      KS1=KS1+3
      KS2=KS2+3
  450 CONTINUE
      KL=KL+2*(ND-I)-1
  451 CONTINUE
C
      RETURN
C
C
C
C
C
  335 DO 340 I=1,ND,NDPN
  340 XM(I+1)=XM(I)
      IF (NDPN.LT.3) RETURN
      DO 350 I=1,ND,3
  350 XM(I+2)=XM(I)
C
      RETURN
C
C
      END
      SUBROUTINE DERIQ (NEL,XX,B,DET,R,S,X1BAR,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
      DIMENSION     XX(2,*),B(4,*),NOD5(*)
      DIMENSION H(9),P(2,9),XJ(2,2),XJI(2,2)
C
      EQUIVALENCE (NPAR(5),ITYP2D)
C
C
C
      IINTP=0
      CALL FUNCT2 (R,S,H,P,NOD5,XJ,DET,XX,NEL,IINTP)
C
C
C
      DUM = 1.0/DET
      XJI(1,1) = XJ(2,2)* DUM
      XJI(1,2) =-XJ(1,2)* DUM
      XJI(2,1) =-XJ(2,1)* DUM
      XJI(2,2) = XJ(1,1)* DUM
C
C
      DO 130 K=1,IEL
      K2=K*2
      B(1,K2-1) = 0.D0
      B(1,K2  ) = 0.D0
      B(2,K2-1) = 0.D0
      B(2,K2  ) = 0.D0
      DO 120 I=1,2
      B(1,K2-1) = B(1,K2-1) + XJI(1,I) * P(I,K)
  120 B(2,K2  ) = B(2,K2  ) + XJI(2,I) * P(I,K)
      B(3,K2  ) = B(1,K2-1)
  130 B(3,K2-1) = B(2,K2  )
C
C
      IF (ITYP2D.GT.0) RETURN
C
C
      X1BAR = 0.D0
      DO 50 K=1,IEL
   50 X1BAR = X1BAR + H(K)* XX(1,K)
C
C
      IF(X1BAR.GT.0.1D-7) GO TO 150
C
C
      ND=2*IEL
      DO 140 K=1,ND
  140 B(4,K)=B(1,K)
      RETURN
C
C     NON-ZERO RADIUS
C
  150 DUM = 1.0/X1BAR
      DO 160 K=1,IEL
      K2=K*2
      B(4,K2  ) = 0.D0
  160 B(4,K2-1) = H(K) * DUM
C
      RETURN
C
C
      END
      SUBROUTINE FUNCT2 (R,S,H,P,NOD5,XJ,DET,XX,NEL,IINTP)
C
C
C .                                                                   .
C .   P R O G R A M                                                   .
C .                                                                   .
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /VAR/   NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,IREF,
     1               IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,
     1            ISTAT,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /TODIM/ BET,THIC,DE,IEL,NND5,ISOCOR
      DIMENSION H(*),P(2,*),NOD5(*),IPERM(4),XJ(2,2),XX(2,*)
C
      EQUIVALENCE (NPAR(8),IDEGEN)
C
      DATA IPERM/2,3,4,1/
C
      RP = 1.0 + R
      SP = 1.0 + S
      RM = 1.0 - R
      SM = 1.0 - S
      R2 = 1.0 - R*R
      S2 = 1.0 - S*S
C
C
C
C     4-NODE ELEMENT
C
      H(1) = 0.25* RP* SP
      H(2) = 0.25* RM* SP
      H(3) = 0.25* RM* SM
      H(4) = 0.25* RP* SM
      P(1,1)=0.25*SP
      P(1,2)=-P(1,1)
      P(1,3)=-0.25*SM
      P(1,4)=-P(1,3)
      P(2,1)=0.25*RP
      P(2,2)=0.25*RM
      P(2,3)=-P(2,2)
      P(2,4)=-P(2,1)
C
      IF (IEL.EQ.4) GO TO 80
C
C
      I=0
    2 I=I + 1
      IF (I.GT.NND5) GO TO 40
      NN=NOD5(I) - 4
      GO TO (5,6,7,8,9),NN
C
    5 H(5) = 0.50* R2* SP
      P(1,5)=-R*SP
      P(2,5)=0.50*R2
      GO TO 2
    6 H(6) = 0.50* RM* S2
      P(1,6)=-0.50*S2
      P(2,6)=-RM*S
      GO TO 2
    7 H(7) = 0.50* R2* SM
      P(1,7)=-R*SM
      P(2,7)=-0.50*R2
      GO TO 2
    8 H(8) = 0.50* RP* S2
      P(1,8)=0.50*S2
      P(2,8)=-RP*S
      GO TO 2
    9 H(9)=R2*S2
      P(1,9)=-2.*R*S2
      P(2,9)=-2.*S*R2
C
C
   40 IH=0
   41 IH=IH + 1
      IF (IH.GT.NND5) GO TO 50
      IN=NOD5(IH)
      IF (IN.EQ.9) GO TO 46
      I1=IN - 4
      I2=IPERM(I1)
      H(I1)=H(I1) - 0.5*H(IN)
      H(I2)=H(I2) - 0.5*H(IN)
      H(IH + 4)=H(IN)
      DO 45 J=1,2
      P(J,I1)=P(J,I1) - 0.5*P(J,IN)
      P(J,I2)=P(J,I2) - 0.5*P(J,IN)
   45 P(J,IH + 4)=P(J,IN)
      GO TO 41
   46 DO 47 I=1,4
      H(I)=H(I) - 0.25*H(9)
      P(1,I)=P(1,I) - 0.25*P(1,9)
   47 P(2,I)=P(2,I) - 0.25*P(2,9)
      N51=NND5 - 1
      IH=0
   48 IH=IH + 1
      IF (IH.GT.N51) GO TO 49
      IN=NOD5(IH)
      I1=IN - 4
      I2=IPERM(I1)
      H(IH+4)=H(IH+4) - 0.5*H(9)
      H(I1  )=H(I1  ) + 0.25*H(9)
      H(I2  )=H(I2  ) + 0.25*H(9)
      DO 52 J=1,2
      P(J,IH+4)=P(J,IH+4) - 0.50*P(J,9)
      P(J,I1  )=P(J,I1  ) + 0.25*P(J,9)
   52 P(J,I2  )=P(J,I2  ) + 0.25*P(J,9)
      GO TO 48
   49 H(NND5+4)=H(9)
      P(1,NND5+4)=P(1,9)
      P(2,NND5+4)=P(2,9)
C
C
   50 IF (IDEGEN.LE.0) GO TO 80
      IF (ISOCOR.LE.0) GO TO 80
C
      DH2D=R2*S2
      H(2)=H(2) + 0.125*DH2D
      H(3)=H(3) + 0.125*DH2D
      H(6)=H(6) - 0.25*DH2D
C
      P(1,2)=P(1,2) - 0.25*R*S2
      P(2,2)=P(2,2) - 0.25*S*R2
      P(1,3)=P(1,3) - 0.25*R*S2
      P(2,3)=P(2,3) - 0.25*S*R2
      P(1,6)=P(1,6) + 0.5*R*S2
      P(2,6)=P(2,6) + 0.5*S*R2
C
C
   80 IF (IINTP.GT.0) RETURN
      DO 100 I=1,2
      DO 100 J=1,2
      DUM = 0.D0
      DO 90 K=1,IEL
   90 DUM = DUM + P(I,K)* XX(J,K)
  100 XJ(I,J) = DUM
C
C
      DET = XJ(1,1)* XJ(2,2) - XJ(2,1)* XJ(1,2)
      CALL DETCH (NG,NEL,DET,1)
      RETURN
C
      END
      SUBROUTINE MAXMIN (STRESS,P1,P2,AG)
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      DIMENSION STRESS(*)
C
C
      CC = (STRESS(1)+STRESS(2)) * 0.5
      BB = (STRESS(1)-STRESS(2)) * 0.5
      CR =  SQRT(BB**2 + STRESS(3)**2)
      P1 =  CC+CR
      P2 =  CC-CR
      AG=4.5D1
      IF(ABS(BB).LT.0.1D-7) RETURN
C
      AG =  28.648* ATAN2(STRESS(3),BB)
C
      RETURN
C
      END
      SUBROUTINE STSTL (NEL,XX,PROP,C)
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
C
      DIMENSION XX(2,*),PROP(*),C(4,*),D(4,4),T(4,4)
C
      EQUIVALENCE (NPAR(5),ITYP2D),(NPAR(15),MODEL)
C
C
      GO TO (1,2), MODEL
C
C
C
C
    1 YM=PROP(1)
      PV=PROP(2)
      C1=YM/(1+PV)
      B1=C1*PV/(1.-2.*PV)
      A1=B1+C1
      C(1,1)=A1
      C(1,2)=B1
      C(1,3)=0.D0
      C(2,1)=B1
      C(2,2)=A1
      C(2,3)=0.D0
      C(3,1)=0.D0
      C(3,2)=0.D0
      C(3,3)=C1/2.
C
      IF (ITYP2D.EQ.1) RETURN
C
      C(1,4)=B1
      C(2,4)=B1
      C(3,4)=0.D0
      C(4,1)=B1
      C(4,2)=B1
      C(4,3)=0.D0
      C(4,4)=A1
C
      IF (ITYP2D.LT.2) RETURN
C
C
      GO TO 95
C
C
C
C
    2 IF (PROP(3).EQ.0.D0) GO TO 1
      PI=4.D0*ATAN(1.D0)
C
C
      DX = XX(1,2) - XX(1,1)
      DY = XX(2,2) - XX(2,1)
      XL = DX**2 + DY**2
      IF(XL.GT.0.1D-7) GO TO 102
      write(66,2000) NEL
      STOP
  102 XL = SQRT(XL)
      SA = ABS(DY/XL)
      CA = SQRT(1.-SA*SA)
      AL = ATAN2(SA,CA)
      IF(DX.GE.0.D0 .AND. DY.GE.0.D0) P12 =         AL
      IF(DX.LT.0.D0 .AND. DY.GE.0.D0) P12 = PI    - AL
      IF(DX.LT.0.D0 .AND. DY.LT.0.D0) P12 = PI    + AL
      IF(DX.GE.0.D0 .AND. DY.LT.0.D0) P12 = PI*2.0- AL
C
C
      PI2=PI*2.
      IF(ABS(P12).GT.PI2) STOP
      GAM=P12 + BET * PI / 180.0
      IF (GAM.GE.PI2) GAM=GAM-PI2
C
C
      IF(ABS(GAM).LT.0.1D-7) GO TO 202
      SG = SIN(GAM)
      CG = COS(GAM)
      T(1,1) =  CG**2
      T(1,2) =  SG**2
      T(1,3) =  CG* SG
      T(1,4) =  0.D0
      T(2,1) =  T(1,2)
      T(2,2) =  T(1,1)
      T(2,3) = -T(1,3)
      T(2,4) =  0.D0
      T(3,1) =  T(2,3)* 2.0
      T(3,2) = -T(3,1)
      T(3,3) =  T(1,1)- T(1,2)
      T(3,4) =  0.D0
      T(4,1) =  0.D0
      T(4,2) =  0.D0
      T(4,3) =  0.D0
      T(4,4) =  1.D0
  202 CONTINUE
C
C
      DUM = PROP(1)* PROP(2)* PROP(3)* PROP(7)
      IF (DUM.GT.0.1D-7) GO TO 25
      write(66,2010)
      STOP
   25 C(1,1) = 1.0/PROP(1)
      C(2,2) = 1.0/PROP(2)
      C(3,3) = 1.0/PROP(7)
      C(4,4) = 1.0/PROP(3)
      C(1,2) =-PROP(4)* C(2,2)
      C(1,4) =-PROP(5)* C(4,4)
      C(2,4) =-PROP(6)* C(4,4)
      C(1,3) = 0.D0
      C(2,3) = 0.D0
      C(3,4) = 0.D0
      DO 30 I=1,4
      DO 30 J=I,4
   30 C(J,I) = C(I,J)
C
C
      CALL POSINV (C,4,4)
C
C
      IF (ABS(GAM).LT.0.1D-7) GO TO 90
C
C
      DO 60 IR=1,4
      DO 60 IC=1,4
      D(IR,IC) = 0.D0
      DO 50 IN=1,4
   50 D(IR,IC) = D(IR,IC) + T(IN,IR)* C(IN,IC)
   60 CONTINUE
C
C
      DO 80 IR=1,4
      DO 80 IC=IR,4
      C(IR,IC) = 0.D0
      DO 70 IN=1,4
   70 C(IR,IC) = C(IR,IC) + D(IR,IN)* T(IN,IC)
   80 C(IC,IR)=C(IR,IC)
C
   90 IF (ITYP2D.LT.2) RETURN
C
C
   95 DO 110 I=1,3
      A=C(I,4)/C(4,4)
      DO 110 J=I,3
      C(I,J)=C(I,J) - C(4,J)*A
  110 C(J,I)=C(I,J)
      RETURN
C
C
C
 2000 FORMAT(10H0*** ERROR,/
     +       43H ZERO LENGTH BETWEEN NODES 1-2 IN ELEMENT (,I4,1H)  )