a25b.for

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

      SUBROUTINE SHSTRS (RSDCOS,NODSYS,X,DISPI,
     1                  LM,XYZ,IELTD,IELTP,IPST,MATP,
     2                  NDOPT,ETIMV,EDISB,ROTB,PROP,
     2                  BET,WA,ITABLE,THICK,
     3                  ISKEW,NTHT,VNI,VNT,V1,IGLOB,ISHAP,COSXY,
     4                  B,XM,RE,S,EDIS,BV,ETIMV2,ISV,
     5                  NTAB,NCON,IDWA,NDM6,NDM36,
     6                  MXTNOD,MXMNOD,NDMV,NDMX,NEQT)
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            ,NDOFDM,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /DIM/ N0,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14,N15
      COMMON/ELSTP/TIME,IDTHF
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /PRCON/ IDATWR,IPRIC,NPB,IDC,IVC,IAC,IPC,IPNODE(3,15)
      COMMON /PORT/ INPORT,JNPORT,NPUTSV,IBPORT,IFPORT,JDC,JVC,JAC
      COMMON /SHELL1/ DISD(9),IELD,IELP,NPT,IDW,NDROT
      COMMON /SHELL2/ NOD(32),NODM(32),NDOPTM(32)
      COMMON /SHELMT/ D(6,6),STRESS(6),STRAIN(6),IPT,N,IPS,ISVE
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      COMMON /TRANG/ TRLW4(4,3),TRLW7(7,3),TRLWD(13,3),
     1               XGRS(36,2),WGTRS(36)
      COMMON /THINT/ NEWGAU
      COMMON /NEWCOT/ XGNC(7,3),WGTNC(7,3)
      COMMON /MDFRDM/ IDOF(12)
      COMMON /RANDI/ N0A,N1D,IELCPL
      COMMON /ADDB/ NEQL,NEQR,MLA,NBLOCK
      COMMON /SKEW  / NSKEWS
      COMMON /SHELL4/ XJI(3,3),P(3,32),H(32)
      COMMON /SHELL5/ ISHAPE
      COMMON /SHROT/ XJ(3,3),DCA(3,3)
      COMMON /SOL/ NUMNP,NEQ,NWK,NWM,NWC,NUMEST,MIDEST,MAXEST,NSTE,MA
      COMMON /XATKA/ LMID(32)
      COMMON /ISUBST/ ISUBC,NSUBST,NSUB,NTUSE,NEGLS,NEGNLS,NUMNPS,
     1                NODCON,NODRET,IDOFS(12),NDOFS,NEQS,NWKS,MAXESC,
     2                MAS,NSTAPE,ILOA(13),KRSIZM,NEQC
      COMMON /PRSHAP/ KSHAPE
      COMMON /DPR/ ITWO
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /ALED/ ALTH,THETA,TREF
      COMMON /TEGRSH/ NTEGRA,N6C
      COMMON /RFORCE/ IREFOR,NEQSUM,NREQ,NBLANK
      COMMON /DISCON/ NDISCE,NIDM
C
      COMMON A(1)
      REAL A
C
      DIMENSION X(*),LM(NDM6,*),XYZ(NDMX,*),
     1          IELTD(*),IELTP(*),IPST(*),MATP(*),
     2          PROP(NCON,*),WA(IDWA,*),S(*),XM(*),B(*),RE(*),
     3          EDIS(*),ETIMV(*),THICK(MXMNOD,*),NDOPT(MXTNOD,*),
     4          ITABLE(NTAB,*),EDISB(NDMX,*),XXX(96),
     5          RSDCOS(9,*),NODSYS(*),ISKEW(MXTNOD,*),NTHT(*),
     6          C(6,6),VNI(NDMV,*),BV(*),
     7          VNT(NDMV,*),V1(NDMV,*),IGLOB(MXMNOD,*),ANG(2),COSXY(*)
      DIMENSION BET(*)
      DIMENSION ETIMV2(*),ISV(*),ROTB(NDMX,*)
      DIMENSION DISPI (*)
      DIMENSION ISHAP(*)
      DIMENSION ILSK(32)
      DIMENSION REFOR(6)
C
      EQUIVALENCE  (NPAR(2),NUME),(NPAR(3),INDNL),(NPAR(4),IDEATH),
     1             (NPAR(6),NEGSKS),(NPAR(9),IFUNCT),
     2             (NPAR(10),NINTR),(NPAR(11),NINTS),(NPAR(12),NINTT),
     3             (NPAR(15),MODEL),(NPAR(16),NUMMAT),(NPAR(5),ISTRES)
      EQUIVALENCE (NPAR(13),NTABLE)
      EQUIVALENCE (NPAR(19),IPRDV)
C
C
C
C
C
  800 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 811
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'NEWSTEP7',NSUB,NG,NEGL,KSTEP,
     2                       (NPAR(I),I=1,20),TIME
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9050) 'NEWSTEP7',NSUB,NG,NEGL,KSTEP,
     2                       (NPAR(I),I=1,20),TIME
C
 9050 FORMAT ( A,/,3(8I10,/),E20.13 )
C
C
  811 IFORCE=0
      ISIGMA=0
      IPRNT=0
      IF (NTABLE.GE.0) ISIGMA=1
      IF (NTABLE.EQ.-1) IFORCE=1
      IF (ISUBC.EQ.0 .AND. IREFOR.GT.0) IFORCE=1
C
      DO 850 IGOV=1,2
      IGOVE=IGOV
C
      IF (IGOVE.EQ.1 .AND. ISIGMA.EQ.0) GO TO 850
      IF (IGOVE.EQ.2 .AND. IFORCE.EQ.0) GO TO 850
C
      CALL GOVELR (IGOVE)
C
      IDISB=0
      IF (IDEATH.EQ.1 .OR. IDEATH.EQ.3) IDISB=1
      EETA=0.D0
C
      DO 840 N=1,NUME
      IF (IDEATH.EQ.0) GO TO 790
      ETIM=ABS(ETIMV(N))
      IF (IDEATH.EQ.2) GO TO 792
      IF (TIME.LT.ETIM) GO TO 899
      IF (IDEATH.NE.3) GO TO 790
      IF (TIME.GT.ETIMV2(N)) GO TO 899
      GO TO 790
  792 IF (TIME.GT.ETIM) GO TO 899
  790 MTYPE = MATP(N)
      IELD = IELTD(N)
      IELP = IELTP(N)
      NTH = NTHT(N)
      IF (NTH.EQ.0) NTH=1
      ISHAPE = ISHAP(N)
      ND = 3*IELD + NDROT*IELP
      NDX = 3*IELD
      NDM3 = ND*(ND+1)/2
      ND6=3*IELD
      NNMS=0
      DO 799 I=1,IELD
      IF (NDOPT(I,N).GE.0) GO TO 799
      ND6=ND6+2
      NNMS=NNMS + 1
      IF (IGLOB(NNMS,N).EQ.0) ND6=ND6 + 1
  799 CONTINUE
      NDM36 = ND6*(ND6+1)/2
      IPS = IPST(N)
      ISVE=ISV(N)
      IF (IPRI.GT.0) IPS=0
      IF (KPLOTE.GT.0) ISVE=0
      IF (IGOVE.EQ.2) GO TO 852
C
      JREFOR=0
      IF (NTABLE.GE.0) GO TO 854
      IPS=0
      ISVE=0
  854 IF (IPS.EQ.0 .AND. ISVE.EQ.0) GO TO 840
      IF (IPS.EQ.0) GO TO 802
      IPRNT=IPRNT + 1
      IF (IPRNT.NE.1) GO TO 802
      write(66,2020) NG
      IF (ISTRES.EQ.0) write(66,2021)
      IF (ISTRES.EQ.1) write(66,2022)
      GO TO (791,791,793,802),MODEL
  791 IF (ISTRES.EQ.0) write(66,2030)
      IF (ISTRES.EQ.1) write(66,2032)
      GO TO 802
  793 IF (ISTRES.EQ.0) write(66,2031)
      IF (ISTRES.EQ.1) write(66,2033)
      GO TO 802
C
  852 IF (NTABLE.EQ.-1) GO TO 856
      IPS = 0
      ISVE = 0
  856 JREFOR = IREFOR
      IF (ISUBC.EQ.1) JREFOR=0
      IF (IPRI.GT.0 ) JREFOR=0
      IF (IPS.GT.0 .OR. ISVE.GT.0) GO TO 858
C
      IF (JREFOR.EQ.0) GO TO 840
      DO 860 L=1,ND6
  860 IF (LM(L,N).GT.NEQT) GO TO 802
      GO TO 840
C
  858 IF (IPS.EQ.0) GO TO 802
      IPRNT=IPRNT+1
      IF (IPRNT.NE.1) GO TO 802
      write(66,2100) NG
      write(66,2110)
  802 CONTINUE
C
C
C
      I=0
      K=0
      NNMS=0
      DO 805 J=1,IELD
      DO 803 L=1,3
      K=K + 1
      I=I + 1
      EDIS(I)=0.D0
      II=LM(K,N)
      IF (II.EQ.0 .OR. II.GT.NEQT) GO TO 803
      IF (II.LT.0) II=NEQ - II
      EDIS(I)=X(II)
  803 CONTINUE
      IF (NDOPT(J,N).GE.0) GO TO 805
      NNMS=NNMS + 1
      K=K + 2
      IF (IGLOB(NNMS,N).EQ.0) K=K+1
  805 CONTINUE
C
      LL=0
      I=0
      IF(INDNL-2) 821,826,826
  821 IF (IDISB.EQ.1)
     1CALL ROTBTH (ROTB(1,N),X,LM(1,N),NDOPT(1,N),IGLOB(1,N),IELD,NEQ,
     2             NEQT,1)
      DO 823 K=1,IELD
      I=I + 3
      IF (NDOPT(K,N)) 824,823,823
  824 LANG=6*LL + 1
      LVN=3*LL + 1
      NNMS=LL+1
      IF (IGLOB(NNMS,N).EQ.0) GO TO 782
      DO 773 IR=1,2
      ANG(IR)=0.D0
      JJ=LM((I+IR),N)
      IF (JJ.EQ.0 .OR. JJ.GT.NEQT) GO TO 773
      IF (JJ.LT.0) JJ=NEQ-JJ
      ANG(IR)=X(JJ)
  773 CONTINUE
      GO TO 787
  782 CONTINUE
      CALL ANCAL (VNI(LVN,N),V1(LVN,1),ANG,RSDCOS,X,LM(I,N),NEGSKS,
     1            ISKEW(K,N),1)
      I=I + 1
  787 CONTINUE
      IVCOD=1
      CALL RSTNOD (COSXY(LANG),VNI(LVN,N),VNT(LVN,1),V1(LVN,1),ANG,
     1             IVCOD)
      LL=LL + 1
      I=I+2
  823 CONTINUE
      IF (IDISB.EQ.1)
     1CALL ROTBTH (ROTB(1,N),X,LM(1,N),NDOPT(1,N),IGLOB(1,N),IELD,NEQ,
     2             NEQT,2)
      GO TO 829
C
  826 CALL DIRECV (VNI(1,N),VNT(1,N),V1(1,N),IGLOB(1,N),LM(1,N),
     1             DISPI,RSDCOS,IELD,IELP,NDOPT(1,N),NEGSKS,
     2             ISKEW(1,N),1)
      DO 827 K=1,IELD
      IF (NDOPT(K,N)) 828,827,827
  828 LANG=6*LL + 1
      LVN=3*LL + 1
      LL=LL+1
      IVCOD=2
      CALL RSTNOD (COSXY(LANG),VNI(LVN,N),VNT(LVN,N),V1(LVN,N),ANG,
     1             IVCOD)
  827 CONTINUE
  829 IF (NEGSKS.EQ.0) GO TO 845
      IF (ISKEW(1,N).LT.0) GO TO 845
      CALL DIRCOS (RSDCOS,EDIS,ISKEW(1,N),IELD,3,1)
  845 CONTINUE
C
      IF (IDISB.NE.1) GO TO 801
      DO 812 I=1,NDX
  812 EDIS(I) =EDIS(I) -EDISB(I,N)
  801 DO 806 I=1,NDX
  806 XXX(I)=XYZ(I,N)
      IF (IDISB.NE.1) GO TO 809
      DO 804 I=1,NDX
  804 XXX(I)=XXX(I) + EDISB(I,N)
  809 CONTINUE
C
C
C
C
C
      IF (IGOVE.EQ.2) GO TO 862
C
      IF (INDNL.NE.3) GO TO 847
      DO 808 I=1,NDX
  808 XXX(I) = XXX(I) + EDIS(I)
C
  847 CALL MAT1 (PROP(1,MTYPE),C)
      IF (MODEL.NE.1 .AND. MODEL.NE.2) GO TO 831
C
      IF (IPRI.GT.0 .OR. IPS.EQ.0) GO TO 814
      IF (ISHAPE.EQ.0) write(66,2035) N
      IF (ISHAPE.EQ.1) write(66,2036) N
  814 CONTINUE
C
C
      IF (NTABLE.EQ.0) GO TO 831
      IPRSV=IPS
      IF (IPS.EQ.0) IPRSV=ISVE
      DO 830 II=1,16
      M=ITABLE(IPRSV,II)
      IF (M.EQ.0) GO TO 840
      IF (IELD.EQ.4) GO TO 870
      IF (INDNL.LE.1)
     1CALL SHDERV (XXX,B,BV,DET,EVAL3(M,1),EVAL3(M,2),EVAL3(M,3),
     2             NDOPT(1,N),COSXY,THICK(1,NTH),EDIS,VNI(1,N),VNT)
      IF (INDNL.GT.1)
     1CALL SHDERV (XXX,B,BV,DET,EVAL3(M,1),EVAL3(M,2),EVAL3(M,3),
     2             NDOPT(1,N),COSXY,THICK(1,NTH),EDIS,VNI(1,N),VNT(1,N))
C
C
C     TO GLOBAL AXES
C
      IF (MODEL.EQ.2) CALL MAT2 (PROP(1,MTYPE),BET,XYZ(1,N),NDOPT(1,N),
     1                           C)
      CALL MATROT (C,D,1)
      CALL STSTSH
C
C
C
C
      IF (ISTRES.EQ.0) GO TO 820
      CALL SIGROT (STRESS,1,1)
      GO TO 822
  820 IF (INDNL.NE.2) GO TO 822
C
      CALL CAUSHL
      GO TO 822
C
  870 IF (INDNL.LE.1)
     1CALL STRPR (XXX,EVAL3(M,1),EVAL3(M,2),EVAL3(M,3),NDOPT(1,N),
     2            THICK(1,NTH),EDIS,VNI(1,N),VNT,C,WA(1,N),
     3            PROP(1,MTYPE))
      IF (INDNL.GT.1)
     1CALL STRPR (XXX,EVAL3(M,1),EVAL3(M,2),EVAL3(M,3),NDOPT(1,N),
     2            THICK(1,NTH),EDIS,VNI(1,N),VNT(1,N),C,WA(1,N),
     3            PROP(1,MTYPE))
C
C
C
  822 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 830
      IF (ISVE.EQ.0) GO TO 830
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-7',N,M,(STRESS(I),I=1,6),
     2                       (STRAIN(I),I=1,6)
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9060) 'OUTPUT-7',N,M,(STRESS(I),I=1,6),
     2                       (STRAIN(I),I=1,6)
C
 9060 FORMAT ( A,/,2I10,/,(4E20.13) )
C
C
  830 IF (IPRI.EQ.0 .AND. IPS.NE.0) write(66,2040) M,STRESS
      GO TO 840
C
C
C
  831 IPT=0
      IF (MODEL.EQ.1) THETA=0.D0
      IF (MODEL.NE.3) GO TO 832
      IF (IPRI.GT.0) GO TO 832
      IF (ISHAPE.EQ.0) write(66,2035) N
      IF (ISHAPE.EQ.1) write(66,2036) N
  832 CONTINUE
C
      CALL SHBASE (NINTR,NINTS,NINTRS)
C
      DO 939 LXY=1,NINTRS
      E1=XGRS(LXY,1)
      E2=XGRS(LXY,2)
      DO 939 LZ=1,NINTT
      IF (NEWGAU.EQ.0) GO TO 940
      E3=XG(LZ,NINTT)
      GO TO 941
  940 NCNC=NINTT/2
      E3=XGNC(LZ,NCNC)
  941 IPT=IPT + 1
C
      IF (IELD.EQ.4) GO TO 990
      IF (INDNL.LE.1)
     1CALL SHDERV (XXX,B,BV,DET,E1,E2,E3,NDOPT(1,N),COSXY,THICK(1,NTH),
     2             EDIS,VNI(1,N),VNT)
      IF (INDNL.GT.1)
     1CALL SHDERV (XXX,B,BV,DET,E1,E2,E3,NDOPT(1,N),COSXY,THICK(1,NTH),
     2             EDIS,VNI(1,N),VNT(1,N))
C
C
C
      IF (MODEL.NE.3) GO TO 980
      CALL MAT3 (PROP(1,MTYPE),C,E1,E2,E3,THICK(1,NTH),WA(1,N),A(N6C),
     1           A(N6B),NDOPT(1,N))
  980 IF (MODEL.EQ.2) CALL MAT2 (PROP(1,MTYPE),BET,XYZ(1,N),NDOPT(1,N),
     1                           C)
      CALL MATROT (C,D,1)
      CALL STSTSH
C
C
C
C
      IF (MODEL.EQ.4) GO TO 939
      IF (ISTRES.EQ.0) GO TO 920
      CALL SIGROT (STRESS,1,1)
      GO TO 925
  920 IF (INDNL.NE.2) GO TO 925
C
      CALL CAUSHL
C
  925 IF (MODEL.EQ.1 .OR. MODEL.EQ.2) GO TO 930
      IF (IPRI.EQ.0 .AND. IPS.NE.0) write(66,2041) IPT,STRESS,THETA
      GO TO 938
  990 IF (INDNL.LE.1)
     1CALL STRPR (XXX,E1,E2,E3,NDOPT(1,N),THICK(1,NTH),EDIS,VNI(1,N),
     2           VNT,C,WA(1,N),PROP(1,MTYPE))
      IF (INDNL.GT.1)
     1CALL STRPR (XXX,E1,E2,E3,NDOPT(1,N),THICK(1,NTH),EDIS,VNI(1,N),
     2            VNT(1,N),C,WA(1,N),PROP(1,MTYPE))
      IF (MODEL.NE.4) GO TO 925
      GO TO 939
  930 IF (IPRI.EQ.0 .AND. IPS.NE.0) write(66,2040) IPT,STRESS
C
C
  938 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 939
      IF (ISVE.EQ.0) GO TO 939
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-7',N,IPT,(STRESS(I),I=1,6),
     2                       (STRAIN(I),I=1,6),THETA
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9070) 'OUTPUT-7',N,IPT,(STRESS(I),I=1,6),
     2                       (STRAIN(I),I=1,6),THETA
C
 9070 FORMAT ( A,/,2I20,/,(4E20.13) )
C
C
  939 CONTINUE
      GO TO 840
C
C
C
  862 IF (INDNL.NE.0) GO TO 864
      DO 866 I=1,NDM36
  866 S(I)=0.D0
C
  864 DO 865 I=1,ND6
      XM(I)=0.D0
  865 RE(I)=0.D0
C
      IF (IELD.EQ.4) GO TO 816
      IF (INDNL.LE.1)
     1CALL SHSTIF (ND,B,S,XXX,PROP(1,MTYPE),RE,EDIS,WA(1,N),
     1             NDOPT(1,N),THICK(1,NTH),BV,COSXY,VNI(1,N),VNT)
      IF (INDNL.GT.1)
     1CALL SHSTIF (ND,B,S,XXX,PROP(1,MTYPE),RE,EDIS,WA(1,N),
     1             NDOPT(1,N),THICK(1,NTH),BV,COSXY,VNI(1,N),VNT(1,N))
      GO TO 818
C
  816 IF (INDNL.LE.1)
     1CALL SHSTCO (S,RE,B,BV,XXX,EDIS,PROP(1,MTYPE),NDOPT(1,N),
     2             WA(1,N),THICK(1,NTH),COSXY,VNI(1,N),VNT)
      IF (INDNL.GT.1)
     1CALL SHSTCO (S,RE,B,BV,XXX,EDIS,PROP(1,MTYPE),NDOPT(1,N),
     2             WA(1,N),THICK(1,NTH),COSXY,VNI(1,N),VNT(1,N))
  818 CONTINUE
C
      NDAUX = ND
      IF (INDNL.EQ.0)
     1CALL GLOROT (VNI(1,N),V1,S,NDOPT(1,N),ND,IELP,IELD,
     2             IGLOB(1,N),IVCOD,EETA,1)
      IF (INDNL.EQ.1)
     1CALL GLOROT (VNI(1,N),V1,RE,NDOPT(1,N),ND,IELP,IELD,
     2             IGLOB(1,N),IVCOD,EETA,3)
      IF (INDNL.EQ.2)
     1CALL GLOROT (VNT(1,N),V1(1,N),RE,NDOPT(1,N),ND,IELP,IELD,
     2             IGLOB(1,N),IVCOD,EETA,3)
C
      IELDE=IELD
      IF (NEGSKS.EQ.0) GO TO 875
      IF (ISKEW(1,N).LT.0) GO TO 875
      NK=0
      NNMS=0
      DO 844 I=1,IELD
      NK=NK+1
      ILSK(NK)=ISKEW(I,N)
      LMID(NK)=NDOPT(I,N)
      IF (NDOPT(I,N).GT.0) GO TO 844
      NNMS=NNMS + 1
      IF (IGLOB(NNMS,N).EQ.1) GO TO 844
      LMID(NK)=-LMID(NK)
      IELDE=IELDE+1
      NK=NK+1
      ILSK(NK)=ILSK(NK-1)
      LMID(NK)=LMID(NK-1)
  844 CONTINUE
C
  875 IF (INDNL.NE.0) GO TO 868
C
      NNMS=0
      K=0
      DO 883 J=1,IELD
      DO 884 L=1,3
      K=K+1
      II=LM(K,N)
      IF (II.EQ.0 .OR. II.GT.NEQT) GO TO 884
      IF (II.LT.0) II=NEQ-II
      XM(K)=X(II)
  884 CONTINUE
      IF (NDOPT(J,N).GE.0) GO TO 883
      NNMS=NNMS+1
      NRON=2
      IF (IGLOB(NNMS,N).EQ.0) NRON=3
      DO 885 L=1,NRON
      K=K+1