a25b.for

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

      II=LM(K,N)
      IF (II.EQ.0 .OR. II.GT.NEQT) GO TO 885
      IF (II.LT.0) II=NEQ-II
      XM(K)=X(II)
  885 CONTINUE
  883 CONTINUE
C
      IF (NEGSKS.EQ.0) GO TO 871
      IF (ISKEW(1,N).LT.0) GO TO 871
      CALL ATKA (RSDCOS,S,ILSK,IELDE,3)
C
  871 IJ=0
      NU=ND6-1
      DO 886 I=1,NU
      IJ=IJ+1
      JL=I+1
      DO 886 J=JL,ND6
      IJ=IJ+1
      RE(I)=RE(I) + S(IJ)*XM(J)
  886 CONTINUE
C
      IJ=0
      DO 872 I=1,ND6
      DO 872 J=I,ND6
      IJ=IJ+1
      RE(J)=RE(J) + S(IJ)*XM(I)
  872 CONTINUE
      GO TO 869
C
  868 IF (NEGSKS.EQ.0) GO TO 869
      IF (ISKEW(1,N).LT.0) GO TO 869
      CALL DIRCOS (RSDCOS,RE,ILSK,IELDE,3,2)
C
  869 IF (JREFOR.EQ.1) CALL READD (A(N3),RE,LM(1,N),ND6)
      IF (IPS.EQ.0) GO TO 874
C
C
      IJ=0
      write(66,2120) N
      NNMS=0
      DO 876 I=1,IELD
      NIELP=NDOPT(I,N)
      NABS=IABS (NIELP)
      DO 877 J=1,6
  877 REFOR(J)=0.D0
      IF (NIELP.LT.0) NNMS=NNMS + 1
      NDOFN=6
      IF (NIELP.LT.0 .AND. IGLOB(NNMS,N).EQ.1) NDOFN=5
      DO 878 J=1,NDOFN
      IF (J.GT.3 .AND. NIELP.GT.0) GO TO 878
      IJ=IJ+1
      REFOR(J) = RE(IJ)
  878 CONTINUE
      write(66,2130) NABS,(REFOR(K),K=1,6)
  876 CONTINUE
C
C
  874 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 840
      IF (ISVE.EQ.0) GO TO 840
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-7',N,ND6,(RE(I),I=1,ND6)
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9080) 'OUTPUT-7',N,ND6,(RE(I),I=1,ND6)
C
 9080 FORMAT ( A,/,2I10,/,(4E20.13) )
C
C
      GO TO 840
  899 IPS=IPST(N)
      IF (IPRI.GT.0) IPS=0
      IF (IPS.EQ.0) GO TO 840
      IPRNT=IPRNT + 1
      IF (IGOVE.EQ.2) GO TO 897
      IF (NTABLE.LT.0) GO TO 840
      IF (IPRNT.NE.1) GO TO 898
      write(66,2020) NG
      IF (ISTRES.EQ.0) write(66,2021)
      IF (ISTRES.EQ.1) write(66,2022)
      IF (MODEL.EQ.1 .AND. ISTRES.EQ.0) write(66,2030)
      IF (MODEL.EQ.1 .AND. ISTRES.EQ.1) write(66,2032)
      IF (MODEL.EQ.3 .AND. ISTRES.EQ.0) write(66,2031)
      IF (MODEL.EQ.3 .AND. ISTRES.EQ.1) write(66,2033)
      IF (MODEL.EQ.2 .AND. ISTRES.EQ.0) write(66,2030)
      IF (MODEL.EQ.2 .AND. ISTRES.EQ.1) write(66,2032)
  898 write(66,2700) N
      GO TO 840
  897 IF (NTABLE.GE.0) GO TO 840
      IF (IPRNT.NE.1) GO TO 896
      write(66,2100) NG
      write(66,2110)
  896 write(66,2700) N
  840 CONTINUE
  850 CONTINUE
      IF (INDNL.NE.2) GO TO 960
      IF (IPRDV.NE.2) GO TO 960
      IF (ISIGMA.EQ.0 .AND. IFORCE.EQ.0) GO TO 960
      IF (IPRI.GT.0 .AND. KPLOTE.GT.0) GO TO 960
      IF (IPRNT.NE.0) write(66,2093) NG
      DO 959 N=1,NUME
      IF (IDEATH.EQ.0) GO TO 950
      ETIM=ABS(ETIMV(N))
      IF (IDEATH.EQ.2) GO TO 952
      IF (TIME.LT.ETIM) GO TO 959
      IF (IDEATH.NE.3) GO TO 950
      IF (TIME.GT.ETIMV2(N)) GO TO 959
      GO TO 950
  952 IF (TIME.GT.ETIM) GO TO 959
  950 IPS=IPST(N)
      ISVE=ISV(N)
      IF (IPRI.EQ.0 .AND. IPS.NE.0) write(66,2092) N
      IELD=IELTD(N)
      IELP=IELTP(N)
      I1=1
      DO 954 IN=1,IELD
      IF (NDOPT(IN,N).GE.0) GO TO 954
      NN=-NDOPT(IN,N)
      I2=I1 + 1
      I3=I1 + 2
      IF (IPRI.EQ.0 .AND. IPS.NE.0)
     1write(66,2094) NN,VNT(I1,N),VNT(I2,N),VNT(I3,N)
C
C
      IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 955
      IF (ISVE.EQ.0) GO TO 955
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'DIREVECT',N,NN,VNT(I1,N),VNT(I2,N),
     2                       VNT(I3,N)
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9090) 'DIREVECT',N,NN,VNT(I1,N),VNT(I2,N),
     2                       VNT(I3,N)
C
 9090 FORMAT ( A,/,2I10,/,3E20.13 )
C
C
  955 I1=I1 + 3
  954 CONTINUE
  959 CONTINUE
  960 IGOVE=3
      CALL GOVELR (IGOVE)
C
      RETURN
C
C
 2020 FORMAT (1H1,45HS T R E S S   C A L C U L A T I O N S   F O R, 3X,
     1        25HE L E M E N T   G R O U P ,I5,3X,13H( 3/D SHELL ) /)
 2021 FORMAT (/ 56H STRESSES ARE CALCULATED IN THE GLOBAL COORDINATE SYS
     1TEM //)
 2022 FORMAT (/ 55H STRESSES ARE CALCULATED IN THE LOCAL COORDINATE SYST
     1EM  //)
 2030 FORMAT (8H ELEMENT,4X,6HOUTPUT,/ 2X,6HNUMBER,2X,8HLOCATION,6X,
     2        9HSTRESS-XX,6X,9HSTRESS-YY,6X,9HSTRESS-ZZ,6X,9HSTRESS-XY,
     3        6X,9HSTRESS-XZ,6X,9HSTRESS-YZ / 1X)
 2031 FORMAT (5X,7HELEMENT,6X,9HSTRESS-XX,6X,
     2        9HSTRESS-YY,6X,9HSTRESS-ZZ,6X,9HSTRESS-XY,6X,9HSTRESS-XZ,
     3        6X,9HSTRESS-YZ,4X,11HTEMPERATURE / 5X,7HNUM/IPT /)
 2032 FORMAT (8H ELEMENT,4X,6HOUTPUT,/ 2X,6HNUMBER,2X,8HLOCATION,6X,
     2        9HSTRESS-RR,6X,9HSTRESS-SS,6X,9HSTRESS-TT,6X,9HSTRESS-RS,
     3        6X,9HSTRESS-RT,6X,9HSTRESS-ST /)
 2033 FORMAT (5X,7HELEMENT,6X,9HSTRESS-RR,6X,
     2        9HSTRESS-SS,6X,9HSTRESS-TT,6X,9HSTRESS-RS,6X,9HSTRESS-RT,
     3        6X,9HSTRESS-ST,4X,11HTEMPERATURE / 5X,7HNUM/IPT /)
 2035 FORMAT (I8)
 2036 FORMAT (I8,3X,10H(TRIANGLE) )
 2040 FORMAT (13X,I5,6E15.4)
 2041 FORMAT (7X,I5,7E15.4)
 2093 FORMAT (1H1,43HN O D A L   D I R E C T O R   V E C T O R S,3X,
     1        36HF O R   E L E M E N T   G R O U P   ,I5,//,
     2        5X,7HELEMENT,5X,4HNODE,5X,5HVNT-X,5X,5HVNT-Y,5X,5HVNT-Z,/)
 2092 FORMAT (5X,I7)
 2094 FORMAT (17X,I4,3F10.5)
 2100 FORMAT(1H1,49H N O D A L   F O R C E   C A L C U L A T I O N     ,
     1           35HF O R   E L E M E N T   G R O U P  ,3X,I2,
     1           3X,8H(SHELLS) ,/)
 2110 FORMAT(8H ELEMENT,5X,5HLOCAL,/,
     1       7H NUMBER,6X,4HNODE,11X,7HFORCE-X,10X,7HFORCE-Y,10X,
     1       7HFORCE-Z,9X,8HMOMENT-X,9X,8HMOMENT-Y,9X,8HMOMENT-Z,/)
 2120 FORMAT(/,I5)
 2130 FORMAT (11X,I5,3X,6(5X,E12.4))
 2700 FORMAT (/,17(2H .),7HELEMENT,I5,3X,28HNOT ACTIVE AT THIS TIME STEP
     1        ,/)
C
C
      END
      SUBROUTINE VNIDET (IELD,IELP,NDOPTM,NODM,NMIDSS,IDV,DV,IGLS,
     1                   VCOMP,IGLOB,
     1                   XYZ,THICK,VNI,N,NG)
C
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /SHELL4/ XJI(3,3),P(3,32),H(32)
      COMMON /SHROT/ XX(3,3),DCA(3,3)
      COMMON /SHELL5/ ISHAPE
C
      DIMENSION NDOPTM(*),NODM(*),IDV(*),DV(3,*),XYZ(*),THICK(*),
     1          IGLS(*),VCOMP(3,*),IGLOB(*),VNI(*)
      DIMENSION RV(16),SV(16),GR(3),GS(3)
C
C
C
      PI=3.141592654D0
      DO 01 I=1,16
      RV(I)=0.D0
   01 SV(I)=0.D0
      RV(1)=1.D0
      SV(1)=1.D0
      RV(2)=-1.D0
      SV(2)=1.D0
      RV(3)=-1.D0
      SV(3)=-1.D0
      RV(4)=1.D0
      SV(4)=-1.D0
      RV(5)=0.D0
      SV(5)=1.D0
      RV(6)=-1.D0
      SV(6)=0.D0
      RV(7)=0.D0
      SV(7)=-1.D0
      RV(8)=1.D0
      SV(8)=0.D0
      IF (IELD.EQ.4) GO TO 30
      I13=0
      DRS=1.D0/6.D0
      DO 29 I=5,IELD
      IF (NDOPTM(I).GT.0) GO TO 29
      JJ=-NDOPTM(I)
      IF (JJ-13) 27,29,28
   27 IF (JJ.LE.8) GO TO 29
      JJI=JJ - 8
      GO TO (9,10,11,12) JJI
    9 RV(JJ)=-DRS
      SV(JJ)=1.D0
      RV(5)=DRS
      GO TO 29
   10 RV(JJ)=-1.D0
      SV(JJ)=-DRS
      SV(6)=DRS
      GO TO 29
   11 RV(JJ)=DRS
      SV(JJ)=-1.D0
      RV(7)=-DRS
      GO TO 29
   12 RV(JJ)=1.D0
      SV(JJ)=DRS
      SV(8)=-DRS
      GO TO 29
   28 I13=I13 + 1
   29 CONTINUE
      IF (I13.EQ.0) GO TO 30
      RV(13)=DRS
      SV(13)=DRS
      RV(14)=-DRS
      SV(14)=DRS
      RV(15)=-DRS
      SV(15)=-DRS
      RV(16)=DRS
      SV(16)=-DRS
C
C
   30 NMS=0
      DO 59 I=1,IELD
      JJ=NDOPTM(I)
      IF (JJ.GT.0) GO TO 59
      JJ=IABS(JJ)
      NMS=NMS + 1
      IVX=3*(NMS-1)+1
      IVY=IVX+1
      IVZ=IVY+1
      NN=NODM(JJ)
      IGLOB(NMS)=IGLS(NN)
      IF (NMIDSS.EQ.0) GO TO 40
      NCOSDV=IDV(NN)
      IF (NCOSDV.EQ.0) GO TO 40
C
C
      VNI(IVX)=DV(1,NCOSDV)
      VNI(IVY)=DV(2,NCOSDV)
      VNI(IVZ)=DV(3,NCOSDV)
      GO TO 50
C
C
   40 R=RV(JJ)
      S=SV(JJ)
      T=0.D0
      IF (ISHAPE.EQ.0) GO TO 41
      IF (JJ.EQ.1 .OR. JJ.EQ.4) GO TO 70
   41 CALL SHFUNT (R,S,T,NDOPTM,DET,XYZ,VNI,THICK,1)
      DO 45 JI=1,3
      GR(JI)=0.D0
   45 GS(JI)=0.D0
C
C
      DO 46 J=1,3
      IXYZ=J
      DO 46 L=1,IELD
      GR(J)=GR(J) + P(1,L)*XYZ(IXYZ)
      GS(J)=GS(J) + P(2,L)*XYZ(IXYZ)
   46 IXYZ=IXYZ + 3
C
C
   60 VNX=GR(2)*GS(3) - GR(3)*GS(2)
      VNY=GR(3)*GS(1) - GR(1)*GS(3)
      VNZ=GR(1)*GS(2) - GR(2)*GS(1)
C
      VN=SQRT (VNX*VNX + VNY*VNY + VNZ*VNZ)
      IF (VN.GT.1.D-10) GO TO 47
      write(66,3000) N,NG
      STOP
   47 VNI(IVX)=VNX/VN
      VNI(IVY)=VNY/VN
      VNI(IVZ)=VNZ/VN
C
C
   50 SUM=0.D0
      DO 52 J=1,3
   52 SUM=SUM + VCOMP(J,NN)
      IF (SUM.NE.0.0D0) GO TO 53
      VCOMP(1,NN)=VNI(IVX)
      VCOMP(2,NN)=VNI(IVY)
      VCOMP(3,NN)=VNI(IVZ)
      GO TO 59
   53 IF (IGLS(NN).NE.1) GO TO 59
      PROD=VCOMP(1,NN)*VNI(IVX) + VCOMP(2,NN)*VNI(IVY) +
     1     VCOMP(3,NN)*VNI(IVZ)
      AA=VCOMP(1,NN)*VCOMP(1,NN) + VCOMP(2,NN)*VCOMP(2,NN) +
     1   VCOMP(3,NN)*VCOMP(3,NN)
      AA=SQRT(AA)
      BB=VNI(IVX)*VNI(IVX) + VNI(IVY)*VNI(IVY) + VNI(IVZ)*VNI(IVZ)
      BB=SQRT(BB)
      COSAL=PROD/(AA*BB)
      COSRE=COS(PI/180.)
      IF (COSAL.GE.COSRE) GO TO 58
      IGLS(NN)=99
      IGLOB(NMS)=IGLS(NN)
      GO TO 59
   58 VNI(IVX)=VCOMP(1,NN)
      VNI(IVY)=VCOMP(2,NN)
      VNI(IVZ)=VCOMP(3,NN)
      GO TO 59
C
C
   70 IF (JJ.EQ.4) GO TO 100
      DO 72 IT=1,3
      GR(IT)=0.D0
   72 GS(IT)=0.D0
      NNR=1
      NNS=1
      DO 75 IT=1,IELD
      JJT=NDOPTM(IT)
      IF (JJT.GT.0) GO TO 75
      JJT=IABS(JJT)
      IF (JJT.EQ.5 .OR. JJT.EQ.9) NNR=NNR + 1
      IF (JJT.EQ.7 .OR. JJT.EQ.11) NNS=NNS + 1
   75 CONTINUE
      GO TO (76,77,78) NNR
   80 GO TO (86,87,88) NNS
   76 GR(1)=(XYZ(1) - XYZ(4))/2.
      GR(2)=(XYZ(2) - XYZ(5))/2.
      GR(3)=(XYZ(3) - XYZ(6))/2.
      GO TO 80
   77 GR(1)=0.5*XYZ(4) + 1.5*XYZ(1) - 2.*XYZ(13)
      GR(2)=0.5*XYZ(5) + 1.5*XYZ(2) - 2.*XYZ(14)
      GR(3)=0.5*XYZ(6) + 1.5*XYZ(3) - 2.*XYZ(15)
      GO TO 80
   78 DO 79 ITI=6,IELD
      JJI=IABS(NDOPTM(ITI))
      IF (JJI.EQ.9) IN9=3*(ITI-1) + 1
   79 CONTINUE
      GR(1)=-7.*XYZ(4)/8.+7.*XYZ(1)/4.+9.*XYZ(IN9  )/4.-9.*XYZ(13)/2.
      GR(2)=-7.*XYZ(5)/8.+7.*XYZ(2)/4.+9.*XYZ(IN9+1)/4.-9.*XYZ(14)/2.
      GR(3)=-7.*XYZ(6)/8.+7.*XYZ(3)/4.+9.*XYZ(IN9+2)/4.-9.*XYZ(15)/2.
      GO TO 80
   86 GS(1)=(XYZ(1)-XYZ(7))/2.
      GS(2)=(XYZ(2)-XYZ(8))/2.
      GS(3)=(XYZ(3)-XYZ(9))/2.
      GO TO 60
   87 DO 89 ITI=5,IELD
      JJI=IABS(NDOPTM(ITI))
      IF (JJI.EQ.7) IN7=3*(ITI-1)+1
   89 CONTINUE
      GS(1)=0.5*XYZ(7) + 1.5*XYZ(1) - 2.*XYZ(IN7)
      GS(2)=0.5*XYZ(8) + 1.5*XYZ(2) - 2.*XYZ(IN7+1)
      GS(3)=0.5*XYZ(9) + 1.5*XYZ(3) - 2.*XYZ(IN7+2)
      GO TO 60
   88 DO 91 ITI=5,IELD
      JJI=IABS(NDOPTM(ITI))
      IF (JJI.EQ.7) IN7=3*(ITI-1)+1
      IF (JJI.EQ.11) IN11=3*(ITI-1)+1
   91 CONTINUE
      GS(1)=-7.*XYZ(7)/8.+7.*XYZ(1)/4.+9.*XYZ(IN7  )/4.-9.*XYZ(IN11  )
     1/2.
      GS(2)=-7.*XYZ(8)/8.+7.*XYZ(2)/4.+9.*XYZ(IN7+1)/4.-9.*XYZ(IN11+1)
     1/2.
      GS(3)=-7.*XYZ(9)/8.+7.*XYZ(3)/4.+9.*XYZ(IN7+2)/4.-9.*XYZ(IN11+2)
     1/2.
      GO TO 60
  100 VNI(IVX)=VNI(1)
      VNI(IVY)=VNI(2)
      VNI(IVZ)=VNI(3)
      GO TO 50
C
   59 CONTINUE
C
      RETURN
 3000 FORMAT (///,5X,7HELEMENT,I5,13HELEMENT GROUP,I5,/,
     1        10X,20HINCORRECT NODAL DATA,/,10X,4HSTOP)
C
      END
      SUBROUTINE DIRECV (VNI,VNT,V1,IGLOB,LM,X,RSDCOS,IELD,IELP,
     1                   NDOPT,NEGSKS,ISKEW,KKK)
C
C      THIS SUBROUTINE
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /SOL/ NUMNP,NEQ,NWK,NWM,NWC,NUMEST,MIDEST,MAXEST,NSTE,MA
      COMMON /DISCON/ NDISCE,NIDM
C
      DIMENSION VNI(*),VNT(*),V1(*),LM(*),X(*),NDOPT(*),ISKEW(*),
     1          RSDCOS(9,*),IGLOB(*)
      DIMENSION ANG (2)
C
      NEQT=NEQ + NDISCE
      IF (KKK.EQ.1) GO TO 100
C
C     KKK=0
C
      IELP3=3*IELP
      DO 10 I=1,IELP3
   10 VNT(I)=VNI(I)
      DO 30 I=1,IELP
      IVX=3*(I-1) + 1
      IVY=IVX + 1
      IVZ=IVY + 1
      VNX=VNT(IVX)
      VNY=VNT(IVY)
      VNZ=VNT(IVZ)
C
      CALL V1CAL (VNX,VNY,VNZ,V1,V1X,V1Y,V1Z,1)
C
      V1(IVX)=V1X
      V1(IVY)=V1Y
      V1(IVZ)=V1Z
C
   30 CONTINUE
      RETURN
C
C     KKK=1
C
  100 NNMS=0
      II=0
      DO 190 I=1,IELD
      II=II + 3
      IF (NDOPT(I).GE.0) GO TO 190
      NNMS=NNMS + 1
C
C
      KK=3*(NNMS-1) + 1
      IF (IGLOB(NNMS).EQ.0) GO TO 110
      DO 105 IR=1,2
      ANG(IR)=0.D0
      JJ=LM(II+IR)
      IF (JJ.EQ.0 .OR. JJ.GT.NEQT) GO TO 105
      IF (JJ.LT.0) JJ=NEQ-JJ
      ANG(IR)=X(JJ)
  105 CONTINUE
      II=II + 2
      GO TO 115
  110 CONTINUE
      CALL ANCAL (VNT(KK),V1(KK),ANG,RSDCOS,X,LM(II),NEGSKS,ISKEW(I),2)
      II=II + 3
  115 CONTINUE
C
C
      REFANG=0.01D0
      INTER=ABS(ANG(1))/REFANG
      INT2 =ABS(ANG(2))/REFANG
      IF (INT2.GT.INTER) INTER=INT2
      IF (INTER.GT.20) INTER=20
      IF (INTER.LT.1) INTER=1
      XINTER=INTER
      DANG1=ANG(1)/XINTER
      DANG2=ANG(2)/XINTER
      DO 50 ITRO=1,INTER
      VN1=VNT(KK)
      VN2=VNT(KK+1)
      VN3=VNT(KK+2)
C
      V11=V1(KK)
      V12=V1(KK+1)
      V13=V1(KK+2)
C
      V21=VN2*V13 - VN3*V12
      V22=VN3*V11 - VN1*V13
      V23=VN1*V12 - VN2*V11
      DUM=SQRT(V21*V21 + V22*V22 + V23*V23)
      DUMI=1./DUM
      V21=V21*DUMI