a25b.for

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

      V22=V22*DUMI
      V23=V23*DUMI
C
      VN1R=VN1 - DANG1*V21 + DANG2*V11
      VN2R=VN2 - DANG1*V22 + DANG2*V12
      VN3R=VN3 - DANG1*V23 + DANG2*V13
      DUM=SQRT(VN1R*VN1R + VN2R*VN2R + VN3R*VN3R)
      DUMI=1./DUM
      VNT(KK)=VN1R*DUMI
      VNT(KK+1)=VN2R*DUMI
      VNT(KK+2)=VN3R*DUMI
C
      V11R=V11 - DANG2*VN1
      V12R=V12 - DANG2*VN2
      V13R=V13 - DANG2*VN3
      DUM=SQRT(V11R*V11R + V12R*V12R + V13R*V13R)
      DUMI=1./DUM
      V1(KK)=V11R*DUMI
      V1(KK+1)=V12R*DUMI
   50 V1(KK+2)=V13R*DUMI
  190 CONTINUE
C
      RETURN
C
      END
      SUBROUTINE GLOROT (VN,V1,S,NDOPT,ND,IELP,IELD,IGLOB,IVCOD,ETA,IGR)
C
C
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
C
      DIMENSION VN(*),V1(*),S(*),NDOPT(*),IGLOB(*)
      DIMENSION T(6,5),AUX(5),W(5,5)
C
      IF (IGR-2) 100,200,300
C
C     I G R  =  1
C
C
  100 NDA=ND
      MS=0
      N1=0
      NNPT=0
      DO 50 IN=1,IELD
      II=NDOPT(IN)
      IF (II) 2,50,1
    1 N1=N1+3
      GO TO 50
    2 MS=MS + 1
      IF (IGLOB(MS).EQ.0) GO TO 3
      N1=N1 + 5
      GO TO 50
    3 N2=3*IELD + 2*IELP + NNPT - N1 - 5
      NOREL=NDA*(NDA+1)/2
      NADEL=N1 + 6 + N2
      NFIEL=NOREL + NADEL
      IPOLD=NOREL+1
      IPNEW=NFIEL+1
      IF (N2.EQ.0) GO TO 40
      NELLT=N2*(N2+1)/2
      DO 10 I=1,NELLT
      IPOLD=IPOLD-1
      IPNEW=IPNEW-1
   10 S(IPNEW)=S(IPOLD)
      NZEROS=N2 + 1
      DO 15 I=1,NZEROS
      IPNEW=IPNEW-1
   15 S(IPNEW)=0.D0
      NEL=NOREL-NELLT
      INZE=N2
      II=0
      DO 25 I=1,NEL
      II=II+1
      IPOLD=IPOLD-1
      IPNEW=IPNEW-1
      S(IPNEW)=S(IPOLD)
      IF (II.LT.INZE) GO TO 25
      IPNEW=IPNEW-1
      S(IPNEW)=0.D0
      INZE=INZE+1
      II=0
   25 CONTINUE
   26 N1=N1+6
      NNPT=NNPT+1
      NDA=NDA+1
      GO TO 50
   40 IPNEW=IPNEW-1
      S(IPNEW)=0.D0
      N15=N1+5
      DO 42 I=1,N15
      IPNEW=IPNEW-1
      S(IPNEW)=0.D0
      DO 42 J=1,I
      IPNEW=IPNEW-1
      IPOLD=IPOLD-1
   42 S(IPNEW)=S(IPOLD)
      GO TO 26
   50 CONTINUE
      IF (NNPT.EQ.0) GO TO 198
C
C
      N1=0
      MS=0
      KK=0
      DO 199 IN=1,IELD
      II=NDOPT(IN)
      IF (II) 102,199,101
  101 N1=N1+3
      GO TO 199
  102 MS=MS + 1
      IF (IGLOB(MS).EQ.0) GO TO 103
      N1=N1 + 5
      KK=KK + 3
      GO TO 199
  103 VNX=VN(KK+1)
      VNY=VN(KK+2)
      VNZ=VN(KK+3)
      CALL V1CAL (VNX,VNY,VNZ,V1(KK+1),V1X,V1Y,V1Z,IVCOD)
      V2X=VNY*V1Z - V1Y*VNZ
      V2Y=VNZ*V1X - V1Z*VNX
      V2Z=VNX*V1Y - V1X*VNY
      DUM=SQRT(V2X*V2X+V2Y*V2Y+V2Z*V2Z)
      DUMI=1./DUM
      V2X=V2X*DUMI
      V2Y=V2Y*DUMI
      V2Z=V2Z*DUMI
      KK=KK+3
      N2=NDA - N1 - 6
C
C
      DO 110 I=1,6
      DO 110 J=1,5
  110 T(I,J)=0.D0
      DO 115 I=1,3
  115 T(I,I)=1.D0
      T(4,4)=V1X
      T(4,5)=V2X
      T(5,4)=V1Y
      T(5,5)=V2Y
      T(6,4)=V1Z
      T(6,5)=V2Z
C
C
      IP=N1
      IF (IP.EQ.0) GO TO 142
      INC=N1
  117 DO 120 I=1,5
  120 AUX(I)=S(IP+I)
      DO 130 IC=1,6
      S(IP+IC)=0.D0
      DO 130 I=1,5
  130 S(IP+IC)=S(IP+IC) + AUX(I)*T(IC,I)
      INC=INC-1
      IF (INC.LE.0) GO TO 140
      IP=IP + N2 + INC + 6
      GO TO 117
  140 IP=IP + 5
  142 IF (N1.EQ.0) IP=-(N2+1)
  141 IPP=IP
      DO 145 I=1,5
      IP=IP + N2 + 1
      DO 145 J=I,5
      IP=IP+1
  145 W(I,J)=S(IP)
      DO 150 I=2,5
      JK=I-1
      DO 150 J=1,JK
  150 W(I,J)=W(J,I)
      IP=IPP + 1
      DO 160 I=1,6
      IP=IP + N2
      DO 160 J=I,6
      IP=IP + 1
      S(IP)=0.D0
      DO 160 L=1,5
      DO 160 K=1,5
  160 S(IP)=S(IP) + T(I,L)*W(L,K)*T(J,K)
C
      IPP=IPP + N2 + 6
      IF (N2.EQ.0) GO TO 181
      DO 180 I=1,N2
      IP=IPP + I + 1
      IP=IP - N2 - 6
      INDA=IP
      DO 165 J=1,5
      INDA=INDA + N2 + (7-J)
  165 AUX(J)=S(INDA)
      DO 170 J=1,6
      IP=IP + N2 + (7-J)
      S(IP)=0.D0
      DO 170 K=1,5
  170 S(IP)=S(IP) + T(J,K)*AUX(K)
  180 CONTINUE
  181 N1=N1+6
      ND=ND + 1
  199 CONTINUE
  198 RETURN
C
C     I G R  =  2
C
  200 N1=0
      MS=0
      FM=2.*ETA/3.
      DO 299 IN=1,IELD
      II=NDOPT(IN)
      IF (II) 230,299,220
  220 N1=N1 + 3
      GO TO 299
  230 N1=N1 + 5
      MS=MS + 1
      IF (IGLOB(MS).EQ.1) GO TO 299
      IOLD=ND + 1
      ND=ND + 1
      INEW=ND + 1
      N1P1=N1 + 1
      NDM1=ND - 1
      IF (N1P1.GT.NDM1) GO TO 250
      K=0
      DO 240 I=N1P1,NDM1
      K=K + 1
  240 S(INEW-K)=S(IOLD-K)
  250 S(N1+1)=S(N1)
      N1=N1 + 1
      S(N1  )=FM*S(N1  )
      S(N1-1)=FM*S(N1-1)
      S(N1-2)=FM*S(N1-2)
  299 CONTINUE
      RETURN
C
C     I G R  =  3
C     TRANSFORM A VECTOR
C
C
  300 NDA=ND
      MS=0
      N1=0
      DO 305 IN=1,IELD
      II=NDOPT(IN)
      IF (II) 302,305,301
  301 N1=N1 + 3
      GO TO 305
  302 MS=MS + 1
      IF (IGLOB(MS).NE.1) GO TO 303
      N1=N1 + 5
      GO TO 305
  303 IOLD=NDA + 1
      NDA=NDA + 1
      INEW=NDA + 1
      N1P6=N1 + 6
      IF (N1P6.EQ.NDA) GO TO 306
      II=N1P6 + 1
      DO 304 I=II,NDA
      IOLD=IOLD-1
      INEW=INEW-1
  304 S(INEW)=S(IOLD)
  306 S(N1P6)=0.D0
      N1=N1 + 6
  305 CONTINUE
      IF (NDA.EQ.ND) GO TO 398
C
C
      N1=0
      MS=0
      KK=0
      DO 399 IN=1,IELD
      II=NDOPT(IN)
      IF (II) 312,399,311
  311 N1=N1 + 3
      GO TO 399
  312 MS=MS + 1
      IF (IGLOB(MS).NE.1) GO TO 313
      N1=N1 + 5
      KK=KK + 3
      GO TO 399
  313 VNX=VN(KK+1)
      VNY=VN(KK+2)
      VNZ=VN(KK+3)
      CALL V1CAL (VNX,VNY,VNZ,V1(KK+1),V1X,V1Y,V1Z,IVCOD)
      V2X=VNY*V1Z - V1Y*VNZ
      V2Y=VNZ*V1X - V1Z*VNX
      V2Z=VNX*V1Y - V1X*VNY
      DUM=SQRT(V2X*V2X+V2Y*V2Y+V2Z*V2Z)
      DUMI=1./DUM
      V2X=V2X*DUMI
      V2Y=V2Y*DUMI
      V2Z=V2Z*DUMI
      KK=KK + 3
C
C
      DO 310 I=1,6
      DO 310 J=1,5
  310 T(I,J)=0.D0
      DO 315 I=1,3
  315 T(I,I)=1.D0
      T(4,4)=V1X
      T(4,5)=V2X
      T(5,4)=V1Y
      T(5,5)=V2Y
      T(6,4)=V1Z
      T(6,5)=V2Z
C
C
      DO 320 I=1,5
  320 AUX(I)=S(N1+I)
      DO 330 I=1,6
      S(N1+I)=0.D0
      DO 330 J=1,5
  330 S(N1+I)=S(N1+I) + T(I,J)*AUX(J)
      N1=N1 + 6
      ND=ND + 1
  399 CONTINUE
  398 RETURN
      END
      SUBROUTINE MATRIT (N,DEN,ETA,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,IREF,
     1             IEQUIT,IPRI,KPLOTN,KPLOTE
      DIMENSION PROP(*)
      EQUIVALENCE (NPAR(15),MODEL),(NPAR(16),NUMMAT),(NPAR(17),NCON)
C
C
      IF (IDATWR.GT.1) GO TO 5
      IF (IMASS.NE.1) write(66,2100) N,DEN
      IF (IMASS.EQ.1) write(66,2102) N,DEN,ETA
    5 CONTINUE
C
      GO TO (1,4,3,2),MODEL
C
C
C
    1 IF (IDATWR.LE.1) write(66,2101) (PROP(I),I=1,NCON)
      IBUG=0
      IF (PROP(1).LE.0.D0) IBUG=1
      IF (PROP(2).GE.0.5D0 .OR. PROP(2).LE.-1.D0) IBUG=1
      IF (IBUG.EQ.0) GO TO 11
C
      write(66,3000) NG,N
      write(66,3050)
      write(66,3100)
      IF (MODEX.EQ.0) GO TO 11
      STOP
   11 CONTINUE
C
      RETURN
C
C
C
    2 IF (NCON.GT.4) GO TO 200
C
C
      IBUG=0
      IF (PROP(1).LE.0.D0) IBUG=1
      IF (PROP(2).GE.0.5D0 .OR. PROP(2).LE.-1.D0) IBUG=1
C
      IF (PROP(3).GT.0.D0) GO TO 150
      IBUG=1
  150 IF (PROP(4).LT.PROP(1) .AND. PROP(4).GE.0.D0) GO TO 152
      IBUG=1
  152 CONTINUE
      IF (IDATWR.LE.1) write(66,2106) (PROP(I),I=1,NCON)
      IF (IBUG.EQ.0) RETURN
      write(66,3000) NG,N
      write(66,3401)
      write(66,3100)
      IF (MODEX.EQ.0) RETURN
      STOP
C
  200 IF (IDATWR.GT.1) GO TO 160
      write(66,2111) (PROP(I),I=1,3)
      write(66,2112) PROP(3),PROP(4)
C
  160 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 162
      IF (PROP(3).GT.0.0D0) GO TO 161
      IBUG=1
  162 write(66,3000) NG,N
      write(66,3401)
      write(66,3100)
  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,3000) NG,N
      write(66,3404) ICP,ICP2
      write(66,3100)
  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
      IF (ET.LT.0.D0) GO TO 175
      IF (ET.GE.PROP(1)) GO TO 175
      IF (ET.LE.ETOLD) GO TO 176
  175 write(66,3000) NG,N
      write(66,3405)
      write(66,3100)
      IF (MODEX.EQ.0) RETURN
      STOP
  176 ETOLD=ET
C
  210 CONTINUE
      RETURN
C
C
C
C
    3 NPTS=INT(PROP(65))
      IF (NPTS.GT.0) GO TO 10
      PROP(65)=16.D0
      NPTS=16
   10 write(66,2110)
      DO 15 KED=1,16
      IP1=KED + 16
      IP2=KED + 32
      IP3=KED + 48
   15 write(66,2120) PROP(KED),PROP(IP1),PROP(IP2),PROP(IP3)
      write(66,2125) PROP(65),PROP(66)
C
      IBUG=0
      IF (NPTS.GE.2 .AND. NPTS.LE.16) GO TO 20
      IBUG=1
      write(66,3000) NG,N
      write(66,2130)
      write(66,3100)
      GO TO 40
C
   20 DO 30 JE=2,NPTS
      JJ=JE-1
      IF (PROP(JE).GE.PROP(JJ)) GO TO 30
      IBUG=1
      write(66,3000) NG,N
      write(66,2140)
      write(66,3100)
      GO TO 40
   30 CONTINUE
C
      DO 22 I=1,NPTS
      IK1=I + 16
      IK2=I + 32
      IF (PROP(IK1).LE.0.D0) IBUG=1
   22 IF (PROP(IK2).LE.-1.D0 .OR. PROP(IK2).GE.0.5D0) IBUG=1
C
      IF (IBUG.EQ.0) RETURN
      write(66,3000) NG,N
      write(66,3150)
      write(66,3050)
      write(66,3100)
C
   40 IF (MODEX.EQ.0) RETURN
      STOP
C
C
    4 IF (IDATWR.LE.1) write(66,2400) (PROP(I),I=1,NCON)
      IBUG=0
      IF (PROP(1).LE.0.D0 .OR. PROP(3).LE.0.D0) IBUG=IBUG+1
      IF (PROP(4).LE.0.D0 .OR. PROP(5).LE.0.D0) IBUG=IBUG+1
      R1=PROP(1)*PROP(3)
      R2=PROP(2)*PROP(2)
      IF (R2.GE.R1) IBUG=IBUG+1
      IF (IBUG.EQ.0) RETURN
      write(66,3450) NG,N
      IF (MODEX.EQ.0) RETURN
      STOP
C
C
 2100 FORMAT (30H MATERIAL CONSTANTS SET NUMBER,6H .... ,I5//,
     1        1H ,4X,29HDEN ..........( DENSITY ).. =, E14.6/)
 2101 FORMAT (1H ,4X,29HE ............( PROP(1) ).. =, E14.6/,
     1        1H ,4X,29HVNU ..........( PROP(2) ).. =, E14.6/,
     2        1H ,4X,29HRKAPA ........( PROP(3) ).. =, E14.6///)
 2102 FORMAT (//31H MATERIAL CONSTANTS SET NUMBER ,16(1H.),2H =,I5,//,
     1     4X,45H DEN .......................... (DENSITY).. =,E14.6/,
     2     4X,45H ETA  (LUMPED ROTATIONAL MASS MULTIPLIER).. =,E14.6/)
 2106 FORMAT (1H ,4X,29HE ............( PROP(1) ).. =,E14.6/
     1        1H ,4X,29HVNU ..........( PROP(2) ).. =,E14.6/
     2        1H ,4X,29HYIELD ........( PROP(3) ).. =,E14.6/
     3        1H ,4X,29HE (HARDEN) ...( PROP(4) ).. =,E14.6///)
 2110 FORMAT (1H ,4X,17HTEMP.(PROP(1-16)),5X,15HE (PROP(17-32)),5X,
     1        17HVNU (PROP(33-48)),4X,19HALPHA (PROP(49-64)),/)