a25b.for

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

      RETURN
C
   75 IF (NINTRS.GT.4) GO TO 80
      NINTRS=4
      DO 77 I=1,NINTRS
      XGRS(I,1)=2.0*TRLW4(I,1) - 1.0
      TEMP=4.0/(1.0 - XGRS(I,1))
      XGRS(I,2)=TRLW4(I,2)*TEMP - 1.0
   77 WGTRS(I)=TRLW4(I,3)*TEMP
      RETURN
C
   80 IF (NINTRS.GT.9) GO TO 85
      NINTRS=7
      DO 82 I=1,NINTRS
      XGRS(I,1)=2.0*TRLW7(I,1) - 1.0
      TEMP=4.0/(1.0 - XGRS(I,1))
      XGRS(I,2)=TRLW7(I,2)*TEMP - 1.0
   82 WGTRS(I)=TRLW7(I,3)*TEMP
      RETURN
C
   85 NINTRS=13
      DO 87 I=1,NINTRS
      XGRS(I,1)=2.0*TRLWD(I,1) - 1.0
      TEMP=4.0/(1.0 - XGRS(I,1))
      XGRS(I,2)=TRLWD(I,2)*TEMP - 1.0
   87 WGTRS(I)=TRLWD(I,3)*TEMP
C
      RETURN
      END
      SUBROUTINE NWCOTT
C
C
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /NEWCOT/ XGNC(7,3),WGTNC(7,3)
C
      DO 10 I=1,7
      DO 10 J=1,3
      XGNC(I,J)=0.D0
   10 WGTNC(I,J)=0.D0
C
      XGNC(1,1)=-1.D0
      XGNC(3,1)=1.D0
      XGNC(1,2)=-1.D0
      XGNC(2,2)=-0.5D0
      XGNC(4,2)=0.5D0
      XGNC(5,2)=1.D0
      XGNC(1,3)=-1.D0
      XGNC(2,3)=-2.D0/3.D0
      XGNC(3,3)=-1.D0/3.D0
      XGNC(5,3)=-XGNC(3,3)
      XGNC(6,3)=-XGNC(2,3)
      XGNC(7,3)=1.D0
C
      WGTNC(1,1)=1.D0/3.D0
      WGTNC(2,1)=4.D0/3.D0
      WGTNC(3,1)= WGTNC(1,1)
      WGTNC(1,2)=7.D0/45.D0
      WGTNC(2,2)=32.D0/45.D0
      WGTNC(3,2)=12.D0/45.D0
      WGTNC(4,2)= WGTNC(2,2)
      WGTNC(5,2)= WGTNC(1,2)
      WGTNC(1,3)=41.D0/420.D0
      WGTNC(2,3)=216.D0/420.D0
      WGTNC(3,3)=27.D0/420.D0
      WGTNC(4,3)=272.D0/420.D0
      WGTNC(5,3)= WGTNC(3,3)
      WGTNC(6,3)= WGTNC(2,3)
      WGTNC(7,3)= WGTNC(1,3)
C
      RETURN
      END
      SUBROUTINE SHSTIF (ND,B,S,XYZ,PROP,RE,EDIS,WA,NDOPT,THICK,
     1                   BV,COSXY,VNI,VNT)
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            ,NDOFDM,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /SHELL1/ DISD(9),IELD,IELP,NPT,IDW,NDROT
      COMMON /SHELMT/ D(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
      COMMON /SHELL5/ ISHAPE
      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 /ELSTP/ TIME,IDTHF
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /TEGRSH/ NTEGRA,N6C
      COMMON /ORTHOT/ N112B
C
      COMMON A(1)
      REAL A
      DIMENSION B(*),S(*),XYZ(*),PROP(*),RE(*),EDIS(*),WA(*),NDOPT(*)
     1         ,C(6,6),TAU(6),DI(6,6),XXX(96),THICK(*),BV(*)
     2         ,COSXY(*),VNT(*),VNI(*),BVD(18)
C
      EQUIVALENCE (NPAR(3),INDNL),(NPAR(10),NINTR),(NPAR(11),NINTS)
     2           ,(NPAR(12),NINTT),(NPAR(15),MODEL)
C
C
C
C
      NDX=3*IELD
      DO 50 J=1,NDX
   50 XXX(J)=XYZ(J)
      IF (INDNL.LE.2) GO TO 55
      DO 52 J=1,NDX
   52 XXX(J)=XXX(J) + EDIS(J)
C
C
   55 CALL MAT1 (PROP,C)
C
C
      CALL SHBASE (NINTR,NINTS,NINTRS)
C
      IPT=0
      DO 100 LXY=1,NINTRS
      E1=XGRS(LXY,1)
      E2=XGRS(LXY,2)
      RSWGT=WGTRS(LXY)
      DO 100 LZ=1,NINTT
      IF (NEWGAU.EQ.0) GO TO 56
      E3=XG(LZ,NINTT)
      WT=RSWGT*WGT(LZ,NINTT)
      GO TO 58
   56 NCNC=NINTT/2
      E3=XGNC(LZ,NCNC)
      WT=RSWGT*WGTNC(LZ,NCNC)
   58 IPT=IPT + 1
      IF (INDNL.EQ.3) GO TO 310
C
C
C
C
C
      CALL SHDERV (XYZ,B,BV,DET,E1,E2,E3,NDOPT,COSXY,THICK,EDIS,VNI,VNT)
C
C
      IF (MODEL.NE.3) GO TO 69
      CALL MAT3 (PROP,C,E1,E2,E3,THICK,WA,A(N6C),A(N6B),NDOPT)
   69 IF (MODEL.EQ.2) CALL MAT2 (PROP,A(N112B),XYZ,NDOPT,C)
      CALL MATROT (C,D,1)
      IF (IND.GE.4) CALL STSTSH
C
      IF (INDNL.LE.1) GO TO 332
C
C
      FAC=WT*DET
      DO 170 I=1,6
 170  TAU(I)=STRESS(I)*FAC
      L=0
      K=0
      DO 180 KK=1,IELD
      K=K + 3
      I=K-2
      J=K-1
      M=L+6
      N=L+12
      DO 179 II=1,6
      RE(I)=RE(I) + BV(L+II)*TAU(II)
      RE(J)=RE(J) + BV(M+II)*TAU(II)
  179 RE(K)=RE(K) + BV(N+II)*TAU(II)
C
      IF (NDOPT(KK)) 175,175,180
  175 K=K + NDROT
      NL=N + 6
      DO 178 II=1,6
      RE(K-1)=RE(K-1) + BV(NL +II)*TAU(II)
  178 RE(K)=RE(K) + BV(NL+II+6)*TAU(II)
      L=L + 12
  180 L=L + 18
C
      IF (ICOUNT.GT.2) GO TO 100
      IF (IREF.NE.0) GO TO 100
C
C
C
      DO 201 I=1,6
      DO 201 J=1,6
  201 DI(I,J)=D(I,J)*FAC
C
      CALL SHBTDB (B,BV,DI,S,ND)
C
      GO TO 465
C
C
C
C
C
  310 CALL SHDERV (XXX,B,BV,DET,E1,E2,E3,NDOPT,COSXY,THICK,EDIS,VNI,VNT)
C
C
      IF (MODEL.NE.3) GO TO 315
      CALL MAT3 (PROP,C,E1,E2,E3,THICK,WA,A(N6C),A(N6B),NDOPT)
  315 IF (MODEL.EQ.2) CALL MAT2 (PROP,A(N112B),XYZ,NDOPT,C)
      CALL MATROT (C,D,1)
      CALL STSTSH
C
C
  332 FAC=WT*DET
      IF (IND.LT.4) GO TO 379
      DO 340 I=1,6
  340 TAU(I)=STRESS(I)*FAC
      KL=0
      K=0
      DO 350 KK=1,IELD
      K=K + 3
      I=K-2
      J=K-1
      KL=KL + 3
      JL=KL - 1
      IL=KL - 2
      RE(I)=RE(I) + B(IL)*TAU(1) + B(JL)*TAU(4) + B(KL)*TAU(5)
      RE(J)=RE(J) + B(JL)*TAU(2) + B(IL)*TAU(4) + B(KL)*TAU(6)
      RE(K)=RE(K) + B(KL)*TAU(3) + B(IL)*TAU(5) + B(JL)*TAU(6)
      IF (NDOPT(KK)) 320,320,350
  320 NL=6*K
      K=K + NDROT
      KL=KL + 9
      NLL=NL + 6
      DO 345 II=1,6
      RE(K-1)=RE(K-1) + BV(NL + II)*TAU(II)
  345 RE(K)=RE(K) + BV(NLL + II)*TAU(II)
  350 CONTINUE
C
  379 IF (ICOUNT.GT.2) GO TO 100
      IF (IREF.NE.0) GO TO 100
C
C
      DO 380 I=1,6
      DO 380 J=1,6
  380 DI(I,J)=D(I,J)*FAC
C
      CALL SHBTDB (B,BV,DI,S,ND)
C
C
C
C
C
  465 IF (INDNL.LE.1) GO TO 100
      IF (IELP) 480,480,500
  480 KL=1
      DO 491 J=1,ND,3
      DB1=TAU(1)*B(J) + TAU(4)*B(J+1) + TAU(5)*B(J+2)
      DB2=TAU(4)*B(J) + TAU(2)*B(J+1) + TAU(6)*B(J+2)
      DB3=TAU(5)*B(J) + TAU(6)*B(J+1) + TAU(3)*B(J+2)
      KS1=KL
      KS2=KS1+ND-J+1
      KS3=KS2+ND-J
      DO 490 I=J,ND,3
      DUM=B(I)*DB1 + B(I+1)*DB2 + B(I+2)*DB3
      S(KS1)=S(KS1) + DUM
      S(KS2)=S(KS2) + DUM
      S(KS3)=S(KS3) + DUM
      KS1=KS1+3
      KS2=KS2+3
  490 KS3=KS3+3
  491 KL=KL+3*ND-3*J
C
      GO TO 100
C
C
C     MID-SURFACE NODES
C
  500 KBV=0
      KB=0
      DO 510 I=1,IELD
      IF (NDOPT(I)) 515,510,510
  515 LL=0
      DO 517 K=4,6
      DO 517 L=1,3
      LL=LL + 1
      BV(KBV+LL)=B(KB+K)*B(KB+L+6)
  517 BV(KBV+LL+9)=B(KB+K)*B(KB+L+9)
      KBV=KBV + 18
      KB=KB + 9
  510 KB=KB + 3
C
C
      KL=1
      KBJ=-17
      KK=1
      JJ=1
      DO 550 J=1,IELD
      DB1=TAU(1)*B(KK) + TAU(4)*B(KK+1) + TAU(5)*B(KK+2)
      DB2=TAU(4)*B(KK) + TAU(2)*B(KK+1) + TAU(6)*B(KK+2)
      DB3=TAU(5)*B(KK) + TAU(6)*B(KK+1) + TAU(3)*B(KK+2)
      KS1=KL
      KS2=KS1 + ND - JJ + 1
      KS3=KS2 + ND - JJ
      KI=KK
      KB=KBJ + 18
      DO 560 I=J,IELD
      DUM=B(KI)*DB1 + B(KI+1)*DB2 + B(KI+2)*DB3
      S(KS1)=S(KS1) + DUM
      S(KS2)=S(KS2) + DUM
      S(KS3)=S(KS3) + DUM
      IF (NDOPT(I)) 558,558,562
  558 DO 565 L=1,2
      S(KS1+L+2)=S(KS1+L+2) + BV(KB)*DB1 + BV(KB+3)*DB2 + BV(KB+6)*DB3
      S(KS2+L+1)=S(KS2+L+1) + BV(KB+1)*DB1 + BV(KB+4)*DB2 + BV(KB+7)*DB3
      S(KS3+L)=S(KS3+L) + BV(KB+2)*DB1 + BV(KB+5)*DB2 + BV(KB+8)*DB3
  565 KB=KB + 9
      KS1=KS1 + NDROT
      KS2=KS2 + NDROT
      KS3=KS3 + NDROT
      KI=KI + 9
  562 KS1= KS1 + 3
      KS2=KS2 + 3
      KS3=KS3 + 3
  560 KI=KI + 3
      KL=KS3 - 2
      JJ=JJ + 3
C
      IF (NDOPT(J)) 580,580,550
C
  580 KBJ=KBJ + 18
      JJ=JJ + NDROT
      KS4=KL
C
      KB=KBJ - 1
      LL=0
      DO 570 L=1,2
      DO 571 I=1,3
      LL=LL + 1
      BVD(LL)=BV(KB+LL)*TAU(1) + BV(KB+LL+3)*TAU(4)
     +        + BV(KB+LL+6)*TAU(5)
      BVD(LL+3)=BV(KB+LL)*TAU(4) + BV(KB+LL+3)*TAU(2)
     +          + BV(KB+LL+6)*TAU(6)
      BVD(LL+6)=BV(KB+LL)*TAU(5) + BV(KB+LL+3)*TAU(6)
     +          + BV(KB+LL+6)*TAU(3)
  571 CONTINUE
  570 LL=9
C
      LL=0
      DO 572 I=1,2
      TEMP=0.D0
      DO 573 L=1,9
      LL=LL + 1
  573 TEMP=TEMP + BVD(L)*BV(KB+LL)
      S(KS4)=S(KS4) + TEMP
  572 KS4=KS4 + 1
C
      KS5=KS4 + ND - JJ + 1
      TEMP=0.D0
      DO 575 L=10,18
  575 TEMP=TEMP + BVD(L)*BV(KB+L)
      S(KS5)=S(KS5) + TEMP
C
      IF (IELD - J) 550,550,589
  589 KI=KK + 12
      IL=J + 1
      KB=KBJ
      DO 590 I=IL,IELD
      S(KS4)=S(KS4) + BVD(1)*B(KI) + BVD(4)*B(KI+1) + BVD(7)*B(KI+2)
      S(KS4+1)=S(KS4+1) + BVD(2)*B(KI) + BVD(5)*B(KI+1) + BVD(8)*B(KI+2)
      S(KS4+2)=S(KS4+2) + BVD(3)*B(KI) + BVD(6)*B(KI+1) + BVD(9)*B(KI+2)
      S(KS5+1)=S(KS5+1) + BVD(10)*B(KI) + BVD(13)*B(KI+1) +
     1         BVD(16)*B(KI+2)
      S(KS5+2)=S(KS5+2) + BVD(11)*B(KI) + BVD(14)*B(KI+1) +
     1         BVD(17)*B(KI+2)
      S(KS5+3)=S(KS5+3) + BVD(12)*B(KI) + BVD(15)*B(KI+1) +
     1         BVD(18)*B(KI+2)
C
      IF (NDOPT(I)) 591,591,598
C
  591 KB=KB + 18
      KS4=KS4 + NDROT
      KS5=KS5 + NDROT
C
      LL=-1
      DO 593 L=1,2
      TEMP=0.D0
      DO 594 JL=1,9
      LL=LL + 1
  594 TEMP=TEMP + BVD(JL)*BV(KB+LL)
  593 S(KS4+L)=S(KS4+L) + TEMP
C
      LL=-1
      DO 595 L=1,2
      TEMP=0.D0
      DO 596 JL=1,9
      LL=LL + 1
  596 TEMP=TEMP + BVD(JL+9)*BV(KB+LL)
  595 S(KS5+L+1)=S(KS5+L+1) + TEMP
      KI=KI + 9
C
  598 KS4=KS4 + 3
      KS5=KS5 + 3
  590 KI=KI + 3
C
      KL=KS5 + 1
      KK=KK + 9
  550 KK=KK + 3
C
  100 CONTINUE
C
      RETURN
C
      END
      SUBROUTINE SHSTCO (S,RE,BC,B,XYZ,EDIS,PROP,NDOPT,WA,THICK,COSXY
     1                   ,VNI,VNT)
C
C
C     MITC4
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 /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN
      COMMON /GT/ XJT(3,3),EPS(4)
      COMMON /SHELL4/ XJI(3,3),P(3,32),H(32)
      COMMON /SHELMT/ D(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
      COMMON /SHROT/ XJ(3,3),DCA(3,3)
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),R,Z,T
      COMMON /THINT/ NEWGAU
      COMMON /NEWCOT/ XGNC(7,3),WGTNC(7,3)
      COMMON /LOCCAR/ EH(9)
      COMMON /PLSTA/ IPELD
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /TEGRSH/ NTEGRA,N6C
      COMMON /ALED/ ALTH,TEMP,TREF
      COMMON /ORTHOT/ N112B
C
      COMMON A(1)
      REAL A
C
      DIMENSION S(*),RE(*),XYZ(*),EDIS(*),PROP(*),NDOPT(*),THICK(*),
     1          WA(*),COSXY(*),VNI(*),VNT(*)
      DIMENSION BC(*),B(*),GCOV(3,3),GCON(3,3)
      DIMENSION C(6,6),ST(6,6)
      DIMENSION GTDI(3),GSDI(3),GTBI(3),GSBI(3),
     1          GTAI(3),GRAI(3),GTCI(3),GRCI(3)
C
      EQUIVALENCE (NPAR(10),NINTR),(NPAR(11),NINTS),(NPAR(12),NINTT),
     1             (NPAR(15),MODEL)
C
      IF (MODEL.EQ.1) CALL MAT1 (PROP,C)
      IPT=0
C
      DO 500 IGR=1,NINTR
      DO 500 IGS=1,NINTS
      DO 500 IGT=1,NINTT
      IPT=IPT+1
      R=XG(IGR,NINTR)
      Z=XG(IGS,NINTS)
      IF (NEWGAU.EQ.0) GO TO 06
      T=XG(IGT,NINTT)
      WT=WGT(IGR,NINTR)*WGT(IGS,NINTS)*WGT(IGT,NINTT)
      GO TO 08
   06 NCNC=NINTT/2
      T=XGNC(IGT,NCNC)
      WT=WGT(IGR,NINTR)*WGT(IGS,NINTS)*WGTNC(IGT,NCNC)
   08 CONTINUE
C
C
C
C
      CALL SHFUNT (R,Z,T,NDOPT,DET,XYZ,VNI,THICK,0)
C
      IF (MODEL.EQ.3)
     1 CALL MAT3 (PROP,C,R,Z,T,THICK,WA,A(N6C),A(N6B),NDOPT)
      IF (MODEL.EQ.2) CALL MAT2 (PROP,A(N112B),XYZ,NDOPT,C)
      WTD=DET*WT
C
C
      DO 10 I=1,3
      DO 10 J=I,3
      GCOV(I,J)=0.D0
      DO 10 K=1,3
   10 GCOV(I,J)=GCOV(I,J) + XJ(I,K)*XJ(J,K)
      DO 12 I=2,3
      II=I-1
      DO 12 J=1,II
   12 GCOV(I,J)=GCOV(J,I)
      GDET=DET*DET
C
C
      GCON(1,1)=(GCOV(2,2)*GCOV(3,3) - GCOV(2,3)*GCOV(2,3))/GDET
      GCON(1,2)=(GCOV(2,3)*GCOV(1,3) - GCOV(1,2)*GCOV(3,3))/GDET
      GCON(1,3)=(GCOV(1,2)*GCOV(2,3) - GCOV(2,2)*GCOV(1,3))/GDET
      GCON(2,2)=(GCOV(1,1)*GCOV(3,3) - GCOV(1,3)*GCOV(1,3))/GDET
      GCON(2,3)=(GCOV(1,2)*GCOV(1,3) - GCOV(1,1)*GCOV(2,3))/GDET
      GCON(3,3)=(GCOV(1,1)*GCOV(2,2) - GCOV(1,2)*GC