a30.for

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

C
  280 NPROP(N)=MPRP
      IPS(N)=IPST
      ISV(N)=ISVPH
      IELD(N)=IELDT
C
      IF (MODEL.EQ.2) GO TO 290
      NDIR(1,N)=IDIR(1)
      NDIR(2,N)=IDIR(2)
      IF (NSTRES.EQ.1) CS(1,N)=CSS
  290 CONTINUE
C
      ND=6*IELDT
C
C
      K=0
C
      DO 310 I = 1,IELDT
      L=NODEM(I,N)
      LL=0
      DO 310 J = 1,6
      K=K + 1
      LM(K,N)=0
      IF (IDOF(J).EQ.1) GO TO 310
      LL=LL + 1
      LM(K,N)=ID(LL,L)
  310 CONTINUE
C
C
      CALL COLHT (HT,ND,LM(1,N))
C
C
      IF (IDATWR.GT.1) GO TO 330
      IF (MODEL.EQ.1) GO TO 320
      NNPRI=9
      IF (IELDT.LE.9) NNPRI=IELDT
      write(66,2400) N,IELDT,IPST,ISVPH,MPRP,KKK,
     1               (NODEM(I,N),I=1,NNPRI)
      IF (IELDT.GT.9) write(66,2405) (NODEM(I,N),I=10,IELDT)
      GO TO 330
C
  320 write(66,2410) N,IELDT,IPST,ISVPH,MPRP,KKK,
     1               (NODEM(I,N),IDIR(I),I=1,IELDT)
  330 CONTINUE
C
C
      IF (JNPORT.EQ.0 .OR. NPUTSV.EQ.0) GO TO 340
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'ELEMENT9',N,IELDT,IPST,MPRP,ISVPH,
     2                       (NODEM(I,N),I=1,IELDT)
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9000) 'ELEMENT9',N,IELDT,IPST,MPRP,ISVPH,
     2                       (NODEM(I,N),I=1,IELDT)
C
 9000 FORMAT ( A,/,(8I10) )
C
C
  340 CONTINUE
C
      IF (N.EQ.NUME) GO TO 360
      N=N + 1
      DO 350 LL = 1,IELDT
      NODEM(LL,N)=NODEM(LL,N-1) + KKK
  350 CONTINUE
C
      IF (N - M) 230,210,150
C
  360 IF (NEGSKS.EQ.0) GO TO 380
C
      DO 370 N = 1,NUME
      IF (ISKEW(1,N).GE.0) GO TO 380
  370 CONTINUE
C
      write(66,2500) NG,NEGSKS
C
  380 CONTINUE
C
      RETURN
C
C
  400 GO TO (410,500,500,600,550),IND
C
C
C
C
C
  410 DO 490 N = 1,NUME
      IELDT=IELD(N)
      MPRP=NPROP(N)
      IF ( IMS(1,MPRP).NE.1 ) GO TO 490
      ND=6*IELDT
      CALL ECHECK (LM(1,N),ND,ICODE,IUPDT)
      IF (ICODE.EQ.1) GO TO 490
C
C
      CALL GLESTF (STF(1,MPRP),ST,AS,ND,MODEL,NDIR(1,N),
     1             IELDT,NODEM(1,N))
C
C
      IF (NEGSKS.EQ.0) GO TO 430
      IF (ISKEW(1,N).LT.0) GO TO 430
      J=-1
      DO 420 I = 1,IELDT
      J=J + 2
      ILSK(J)=ISKEW(I,N)
      ILSK(J+1)=ILSK(J)
  420 CONTINUE
      ITELD=2*IELDT
      CALL ATKA (RSDCOS,ST,ILSK,ITELD,3)
C
  430 CONTINUE
C
C
      CALL ADDBAN (A(N2),A(N1),ST,RE,LM(1,N),ND,1)
  490 CONTINUE
C
      RETURN
C
C
C
C
  500 DO 540 N=1,NUME
      IELDT=IELD(N)
      MPRP=NPROP(N)
      IF ( IMS(2,MPRP).NE.1 ) GO TO 540
      ND=6*IELDT
C
      IF (IMASS.EQ.1) GO TO 505
      CALL ECHECK (LM(1,N),ND,ICODE,IUPDT)
      IF (ICODE.EQ.1) GO TO 540
C
C
  505 CONTINUE
      CALL GLEMAS (XM(1,MPRP),ST,AS,ND,MODEL,IMASS,
     1             NDIR(1,N),IELDT,NODEM(1,N))
C
C
      IF (NEGSKS.EQ.0) GO TO 520
      IF (ISKEW(1,N).LT.0) GO TO 520
      J=-1
      DO 510 I = 1,IELDT
      J=J + 2
      ILSK(J)=ISKEW(I,N)
      ILSK(J+1)=ILSK(J)
  510 CONTINUE
      ITELD=2*IELDT
      CALL ATKA (RSDCOS,ST,ILSK,ITELD,3)
C
  520 CONTINUE
      IF (IMASS.EQ.2) GO TO 530
C
C     LUMPED MASS MATRIX
C
      CALL ADDMA (A(N4),ST,LM(1,N),ND)
      GO TO 540
C
C
  530 CALL ADDBAN (A(N2),A(N1),ST,RE,LM(1,N),ND,1)
  540 CONTINUE
C
      RETURN
C
C
C
C
  550 DO 590 N=1,NUME
      IELDT=IELD(N)
      MPRP=NPROP(N)
      IF ( IMS(3,MPRP).NE.1 ) GO TO 590
      ND=6*IELDT
C
      IF (IDAMP.EQ.1) GO TO 555
      CALL ECHECK (LM(1,N),ND,ICODE,IUPDT)
      IF (ICODE.EQ.1) GO TO 590
C
C
  555 CONTINUE
      CALL GLEDMP (XC(1,MPRP),ST,AS,ND,MODEL,IDAMP,NDIR(1,N),IELDT,
     A             NODEM(1,N) )
C
C
      IF (NEGSKS.EQ.0) GO TO 570
      IF (ISKEW(1,N).LT.0) GO TO 570
      J=-1
      DO 560 I=1,IELDT
      J=J+2
      ILSK(J)=ISKEW(I,N)
      ILSK(J+1)=ILSK(J)
  560 CONTINUE
      ITELD=2*IELDT
      CALL ATKA (RSDCOS,ST,ILSK,ITELD,3)
C
  570 CONTINUE
      IF (IDAMP.EQ.2) GO TO 580
C
C
      CALL ADDMA (A(N5),ST,LM(1,N),ND)
      GO TO 590
C
C
  580 CALL ADDBAN (A(N2),A(N1),ST,RE,LM(1,N),ND,1)
  590 CONTINUE
C
      RETURN
C
C
C
C
  600 ISTIF=0
      IF (ICOUNT.EQ.3) GO TO 610
      IF (IREF.EQ.0) ISTIF=1
  610 CONTINUE
      MADR=N3
      IF (IEQUIT.EQ.0) MADR=N5
C
      DO 700 N = 1,NUME
      IELDT=IELD(N)
      MPRP=NPROP(N)
      IF ( IMS(1,MPRP).NE.1 ) GO TO 700
      ND=6*IELDT
C
      CALL ECHECK (LM(1,N),ND,ICODE,IUPDT)
      IF (ICODE.EQ.1) IELCPL=IELCPL + 1
      IF (ICODE.EQ.1) GO TO 700
C
C
      ISKEL=0
      IF (NEGSKS.EQ.0) GO TO 615
      IF (ISKEW(1,N).LT.0) GO TO 615
      ISKEL=1
      ITELD=2*IELDT
C
      J=-1
      DO 612 I = 1,IELDT
      J=J + 2
      ILSK(J)=ISKEW(I,N)
      ILSK(J+1)=ILSK(J)
  612 CONTINUE
  615 CONTINUE
C
C
      DO 640 L = 1,ND
      DISP(L)=0.D0
      II=LM(L,N)
      IF (II.GT.NEQT) GO TO 640
      IF (II) 620,640,630
  620 II=NEQ-II
  630 DISP(L)=X(II)
  640 CONTINUE
C
C
      IF (ISKEL.EQ.1) CALL DIRCOS (RSDCOS,DISP,ILSK,ITELD,3,1)
C
C
      CALL GLESTF (STF(1,MPRP),ST,AS,ND,MODEL,NDIR(1,N),
     1             IELDT,NODEM(1,N))
C
C
      DO 670 I = 1,ND
      RE(I)=0.D0
  670 CONTINUE
C
      CALL FORCAL (ST,RE,AS,DISP,ND)
C
C
      IF (ISKEL.EQ.0) GO TO 680
      CALL DIRCOS (RSDCOS,RE,ILSK,ITELD,3,2)
C
  680 CALL ADDBAN (A(MADR),A(N1),ST,RE,LM(1,N),ND,2)
C
      IF (ISTIF.EQ.0) GO TO 700
C
C
C
      IF (ISKEL.EQ.0) GO TO 690
      CALL ATKA (RSDCOS,ST,ILSK,ITELD,3)
C
  690 CALL ADDBAN (A(N4),A(N1),ST,RE,LM(1,N),ND,1)
  700 CONTINUE
C
      IF (IELCPL.EQ.NUME) IELCPL=-1
C
      RETURN
C
C
C
C
  800 CONTINUE
C
C
      IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 810
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'NEWSTEP9',NSUB,NG,NEGL,KSTEP,
     2                       (NPAR(I),I=1,20),TIME
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9010) 'NEWSTEP9',NSUB,NG,NEGL,KSTEP,
     2                       (NPAR(I),I=1,20),TIME
C
 9010 FORMAT ( A,/,3(8I10,/),E20.13 )
C
C
  810 CONTINUE
      IPRNT=0
C
      DO 950 N = 1,NUME
      IPST=IPS(N)
      ISVE=ISV(N)
      IF (IPST.EQ.0 .AND. ISVE.EQ.0 .AND. IREFOR.EQ.0) GO TO 950
      MPRP=NPROP(N)
      IF ( IMS(1,MPRP).NE.1 ) GO TO 950
C
      IF (IPRI.NE.0) GO TO 820
      IF (IPST.EQ.0) GO TO 820
      IPRNT=IPRNT + 1
      IF (IPRNT.NE.1) GO TO 820
      IF (NSTRES.LE.0) write(66,2510) NG
      IF (NSTRES.GT.0) write(66,2520) NG
      IF (NSTRES.LE.0) write(66,2530)
      IF (NSTRES.GT.0) write(66,2540)
C
  820 IELDT=IELD(N)
      ND=6*IELDT
C
C
      ISKEL=0
      IF (NEGSKS.EQ.0) GO TO 825
      IF (ISKEW(1,N).LT.0) GO TO 825
      ISKEL=1
      ITELD=2*IELDT
C
      J=-1
      DO 823 I = 1,IELDT
      J=J + 2
      ILSK(J)=ISKEW(I,N)
      ILSK(J+1)=ILSK(J)
  823 CONTINUE
  825 CONTINUE
C
C
      DO 850 I = 1,ND
      DISP(I)=0.D0
      II=LM(I,N)
      IF (II.GT.NEQT) GO TO 850
      IF (II) 830,850,840
  830 II=NEQ - II
  840 DISP(I)=X(II)
  850 CONTINUE
C
C
      IF (ISKEL.EQ.1) CALL DIRCOS (RSDCOS,DISP,ILSK,ITELD,3,1)
C
C
C
      DO 865 I = 1,ND
      RE(I)=0.D0
  865 CONTINUE
C
      CALL GLESTF (STF(1,MPRP),ST,AS,ND,MODEL,NDIR(1,N),
     1             IELDT,NODEM(1,N))
C
      CALL FORCAL (ST,RE,AS,DISP,ND)
C
C
C
      IF (ISKEL.EQ.1) CALL DIRCOS (RSDCOS,RE,ILSK,ITELD,3,2)
C
C     OBTAIN REACTION
C
      IF ((ISUBC.EQ.0).AND.(IREF.EQ.1)) CALL READD (A(N3),RE,LM(1,N),
     1                                              ND)
      IF (NSTRES.EQ.0) GO TO 880
C
C
      CALL STRCAL (CS(1,MPRP),AS,RE,DISP,NDIR(1,N),NSTRES,
     1             MODEL,ND,IELDT,NODEM(1,N))
      GO TO 930
  880 CONTINUE
C
C
C
      IF (IPRI.NE.0 .OR. IPST.EQ.0) GO TO 910
      write(66,2610) N
C
      DO 920 I = 1,IELDT
      II=NODEM(I,N)
      L=0
      IF (ISKEL.EQ.0) GO TO 890
      IF (ISKEW(I,N).NE.0) L=2
  890 CONTINUE
C
      DO 900 J = 1,6
      K=(I - 1)*6 + J
      FORCE(J)=RE(K)
  900 CONTINUE
C
      write(66,2600) II,AIND(L+1),AIND(L+2),
     1               (FORCE(J),J = 1,6)
C
  920 CONTINUE
C
C
  910 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 950
      IF (ISVE.EQ.0) GO TO 950
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-9',N,ND,(RE(I),I=1,ND)
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9020) 'OUTPUT-9',N,ND,(RE(I),I=1,ND)
C
 9020 FORMAT ( A,/,2I10,/,(4E20.13) )
C
C
      GO TO 950
C
C     PRINT OUT STRESSES
C
  930 IF (IPRI.NE.0 .OR. IPST.EQ.0) GO TO 935
      write(66,2610) N
      write(66,2700) (I,RE(I),I=1,NSTRES)
C
C
  935 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 950
      IF (ISVE.EQ.0) GO TO 950
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-9',N,NSTRES,(RE(I),I=1,NSTRES)
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9030) 'OUTPUT-9',N,NSTRES,(RE(I),I=1,NSTRES)
C
 9030 FORMAT ( A,/,2I10,/,(4E20.13) )
C
C
  950 CONTINUE
C
      RETURN
C
C
 1000 FORMAT (4I5)
 1010 FORMAT (4F10.0)
 1020 FORMAT (8F10.0)
 1100 FORMAT (6I5)
 1200 FORMAT (16I5)
C
 2000 FORMAT (//52H S T I F F N E S S / M A S S / D A M P I N G    M A ,
     A 60HT R I C E S    S E T S                                      /)
C
 2010 FORMAT (//,20H MATRICES SET NUMBER,6H .... ,I5,/)
 2015 FORMAT (1H ,2X,18HSTIFFNESS ...... =,E14.6)
 2016 FORMAT (1H ,2X,18HMASS ........... =,E14.6)
 2017 FORMAT (1H ,2X,18HSTRESS TRANS ... =,E14.6)
 2018 FORMAT (1H ,2X,18HDAMPING ........ =,E14.6)
C
 2300 FORMAT (///21H *** INPUT ERROR *** /,
     1        1H ,4X,18HSUBSTRUCTURE NO. =,I5/,
     2        1H ,4X,18HELEMENT GROUP NO.=,I5/,
     3        1H ,4X,18HELEMENT NUMBER   =,I5/,
     4        1H ,4X,28HNUMBER OF NODES IN ELEMENT (,I5,2H) ,
     5        42HEXCEEDS MAXIMUM ALLOWED NODES IN ELEMENT (,I5,2H).,//,
     6        1H ,10X,8H S T O P )
C
 2310 FORMAT (///21H *** INPUT ERROR *** /,
     1        1H ,4X,18HSUBSTRUCTURE  NO.=,I5/,
     2        1H ,4X,18HELEMENT GROUP NO.=,I5/,
     3        1H ,4X,30HFIRST ELEMENT NUMBER MUST BE 1 ,//,
     4        1H ,10X,8H S T O P )
C
 2320 FORMAT (///21H *** INPUT ERROR *** /,
     1        1H ,4X,18HSUBSTRUCTURE  NO.=,I5/,
     2        1H ,4X,18HELEMENT GROUP NO.=,I5/,
     3        1H ,4X,18HELEMENT NUMBER   =,I5/,
     4        1H ,4X,21HPROPERTY SET NUMBER (,I3,
     5        19H) IS GT. NPAR(16)   ,//,
     6        1H ,10X,8H S T O P )
C
 2330 FORMAT (///21H *** INPUT ERROR *** /,
     1        1H ,4X,18HSUBSTRUCTURE  NO.=,I5/,
     2        1H ,4X,18HELEMENT GROUP NO.=,I5/,
     3        1H ,4X,18HELEMENT NUMBER   =,I5/,
     4        1H ,4X,12H NODE NUMBER ,I5,11H IS LARGER ,
     5        31HTHAN THE TOTAL NUMBER OF NODES ,I5,//,
     6        1H ,10X,8H S T O P )
C
 2340 FORMAT (///21H *** INPUT ERROR *** /,
     1        1H ,4X,18HSUBSTRUCTURE  NO.=,I5/,
     2        1H ,4X,18HELEMENT GROUP NO.=,I5/,
     3        1H ,4X,30HGLOBAL DIRECTION MUST BE GE.1 ,
     4        8HAND LE.6/,
     5        1H ,4X,43HFOR TWO DEGREES OF FREEDOM GENERAL ELEMENTS ,
     6        //1H ,10X,8H S T O P )
C
 2350 FORMAT (///40H E L E M E N T   I N F O R M A T I O N    )
C
 2360 FORMAT (///5X,1HN,3X,4HIELD,3X,3HIPS,2X,5HISVPH,2X,4HPROP,4X,2HKG,
     1        15X,A4,I1,3X,A4,I1,3X,A4,I1,3X,A4,I1)
C
 2370 FORMAT (///5X,1HN,3X,4HIELD,3X,3HIPS,2X,5HISVPH,2X,4HPROP,4X,2HKG,
     1        15X,9(A4,I1,3X))
C
 2375 FORMAT (54X,9(A4,I2,2X)/)
C
 2380 FORMAT (///21H *** INPUT ERROR *** /,
     1        1H ,4X,18HSUBSTRUCTURE  NO.=,I5/,
     2        1H ,4X,18HELEMENT GROUP NO.=,I5/,
     3        1H ,4X,18HELEMENT NUMBER   =,I5/,
     4        1H ,4X,42HSINCE NODES OF THIS ELEMENT REFER TO SKEW ,
     5        47HCOORDINATE SYSTEM(S) NPAR(6) MUST BE SET TO 1   ,
     6        //1H ,10X,8H S T O P )
C
 2400 FORMAT (/I6,3I6,I7,2X,I6,11X,9(3X,I5))
 2405 FORMAT (1H ,49X,9(3X,I5))
C
 2410 FORMAT (/,4I6,I8,I6,15X,I5,2X,I5,4X,I5,2X,I5)
C
 2500 FORMAT (///16H ELEMENT GROUP =,I3,
     1        35H ( GENERAL MASS/STIFFNESS ELEMENT ) /,
     2        19H ALTHOUGH NPAR(6) =,I2,22H, NONE OF THE ELEMENTS/
     3        49H IN THIS GROUP REFERS TO SKEW COORDINATE SYSTEMS./
     4        50H IT IS ADVISED THAT NPAR(6) BE SET TO ZERO TO SAVE,
     5        15H STORAGE SPACE.//
     6        39H THIS IS ONLY AN INFORMATIONAL MESSAGE.///)
C
 2510 FORMAT (1H1,43HF O R C E   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,3X,I2,3X,
     2        20H( GENERAL ELEMENTS ),/)
C
 2520 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,3X,I2,3X,
     2        20H( GENERAL ELEMENTS ),/)
C
 2530 FORMAT (8H ELEMENT,4X,4HNODE,4X,4HREF.,/ 2X,6HNUMBER,2X,
     1        6HNUMBER,2X,6HSYSTEM,7X,8H FORCE-1,7X,8H FORCE-2,
     2        7X,8H FORCE-3,7X,8HMOMENT-1,7X,8HMOMENT-2,7X,
     3        8HMOMENT-3 / 1X)
C
 2540 FORMAT (8H ELEMENT,4X,6HSTRESS,/ 2X,6HNUMBER,4X,6HNUMBER,7X,
     1        8H  STRESS / 1X)
C
 2600 FORMAT (10X,I5,2X,A4,A4,1X,6(1X,E14.6))
 2610 FORMAT (I6)
 2700 FORMAT (10X,I6,6X,E14.6/)
 3001 FORMAT (28H *** I N P U T    E R R O R ,/,
     A 61H DAMPING ELEMENTS ARE NOT ALLOWED UNLESS IDAMP.EQ.1 OR 3     ,
     B //,13H *** S T O P  )
 3801 FORMAT (///,19H SUBSTRUCTURE NO  =,I5)
 3802 FORMAT (///,19H ELEMENT GROUP NO =,I5,
     A          /,19H PROPERTY SET NO  =,I5)
 3901 FORMAT (//37H * * *     STOP OF SOLUTION     * * *,//,
     1          28H INPUT ERROR IN ELEMENT DATA,/,
     2          19H ELEMENT GROUP NO =,I5,/,
     3          13H ELEMENT NO =,I5)
 3902 FORMAT (//37H * * *     STOP OF SOLUTION     * * *,//,
     1          28H INPUT ERROR IN ELEMENT DATA,/,
     *          18H SUBSTRUCTURE NO =,I5,/,
     2          19H ELEMENT GROUP NO =,I5,/,
     3          13H ELEMENT NO =,I5)
 3905 FORMAT (  41H   NPROP MUST BE .GE. 1 AND .LE. NPAR(16),/,
     1          10H   NPROP =,I5)
 3920 FORMAT (/ 40H SEE USER*S MANUAL FOR INPUT DESCRIPTION,//,
     1          37H * * *   END OF ERROR MESSAGE   * * *,//)
C
      END
      SUBROUTINE GLESTF (STF,ST,AS,ND,MODEL,NDIR,IELDT,NODEM)
C
C
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /MDFRDM/ IDOF(12)
C
      DIMENSION STF(*),ST(*),AS(ND,*),NDIR(*),NODEM(*)
C