a37.for

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

      SUBROUTINE THDCON
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /DIM/ N0,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14,N15
      COMMON /DIMETC/ N01,N02,N03,N04,N05,N06,N07,N08,N09
      COMMON /FREQIF/ ISTOH,N1A,N1B,N1C,N1S
      COMMON /LOA/ NLOAD,NPR2,NPR3,NPBM,NP3DB,NPPL,NPSH,NODE3,IDGRAV,
     1             NPDIS,NTEMP,IDCFGL,ISMASS
      COMMON /SOL/ NUMNP,NEQ,NWK,NWM,NWC,NUMEST,MIDEST,MAXEST,NSTE,MA
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1           ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NGALL,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /EQIT/ METHOD,ILNSCH,NLSTPD,NLSTEP,ITEDIV,IPEQIT
      COMMON /SEQNCE/ ITA(18),IT19,ITB(10)
      COMMON /ADDB/ NEQL,NEQR,MLA,NBLOCK
      COMMON /DPR/ ITWO
      COMMON /ELGLTH/ NFIRST,NLAST,NBCEL
      COMMON /ELSTP/ TIME,IDTHF
      COMMON /PRCON/ IDATWR,IPRIC,NPB,IDC,IVC,IAC,IPC,IPNODE(3,15)
      COMMON /SKEW/ NSKEWS
      COMMON /PORT/ INPORT,JNPORT,NPUTSV,IBPORT,IFPORT,JDC,JVC,JAC
      COMMON /CNTACT/ NEQI,LEADEQ,NCE2D,NCE3D
      COMMON /ITRATE/ RCENRM,RCNORM,RCONSM,RCTOL
      COMMON /GROUPS/ NTGNL,NSGNL
      COMMON /ISET3/ ISET,NCEQN
      COMMON /ISSAC/ NEWREF
      COMMON /TAPES/ IIN,IOUT
      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 A(1)
      REAL A
      DIMENSION IA(1)
      DIMENSION IMDUMY(1,1),RMDUMY(1,1)
      EQUIVALENCE ( A(1),IA(1) )
C
      EQUIVALENCE ( NPAR(1),NPAR1  ),( NPAR(2),NSURF  ),
     1            ( NPAR(4),NSURFP ),( NPAR(6),NEGSKS ),
     1            ( NPAR(7),NECON  ),( NPAR(8),NEALL  ),
     1            ( NPAR(9),NTOUCH ),( NPAR(10),NSNOD ),
     1            ( NPAR(11),NJOIN ),( NPAR(13),ICPRNT),
     1            ( NPAR(15),MODEL )
C
      IF (IND.GE.1 .AND. IND.LE.3) RETURN
      IF (IND.EQ.4) GO TO 500
      NG = NGALL-NEGL-NEGNL
C
C
C
      ISTOP=0
      IF (NEWREF.EQ.1) GO TO 30
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,2550) ISTOP
C
   30 IF (ILNSCH.EQ.0) GO TO 40
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,2560) ISTOP
C
C
   40 IF ( NSURF.GE.2 ) GO TO 50
      ISTOP=ISTOP+1
      IF ( ISTOP.EQ.1 ) write(66,2500) NG
      write(66,2600) ISTOP,NSURF
C
   50 IF ( NSURFP.GE.1 ) GO TO 100
      ISTOP=ISTOP+1
      IF ( ISTOP.EQ.1 ) write(66,2500) NG
      write(66,2700) ISTOP,NSURFP
C
  100 IF (NEGSKS.EQ.0) GO TO 150
      IF (NSKEWS.GT.0) GO TO 150
      ISTOP=ISTOP+1
      IF ( ISTOP.EQ.1 ) write(66,2500) NG
      write(66,2800) ISTOP,NSKEWS,NEGSKS
C
  150 IF ( NECON.GE.1 ) GO TO 200
      ISTOP=ISTOP+1
      IF ( ISTOP.EQ.1 ) write(66,2500) NG
      write(66,2900) ISTOP,NECON
C
  200 IF ( NEALL.GE.2 ) GO TO 310
      ISTOP=ISTOP+1
      IF ( ISTOP.EQ.1 ) write(66,2500) NG
      write(66,3000) ISTOP,NEALL
C
  310 IF ( NTOUCH.GE.3 ) GO TO 320
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,3210) ISTOP,NTOUCH
C
  320 IF ( NSNOD.GE.6 ) GO TO 330
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,3220) ISTOP,NSNOD
C
  330 IF ( NJOIN.EQ.0 ) NJOIN = 4
      IF ( NJOIN.GE.1 ) GO TO 350
      ISTOP=ISTOP+1
      IF ( ISTOP.EQ.1 ) write(66,2500) NG
      write(66,3250) ISTOP,NJOIN
C
  350 ICTEMP=ICPRNT
      IF (ICTEMP.EQ.1) ICPRNT=2
      IF (ICTEMP.EQ.2) ICPRNT=1
      IF (ICTEMP.GE.0 .AND. ICTEMP.LE.2) GO TO 360
      ISTOP=ISTOP+1
      IF ( ISTOP.EQ.1 ) write(66,2500) NG
      write(66,3300) ISTOP,ICTEMP
C
  360 MTEMP = MODEL
      MODEL = MODEL + 1
      IF (MTEMP.EQ.0 .OR. MTEMP.EQ.1) GO TO 400
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,3310) ISTOP,MTEMP
C
  400 IF ( ISTOP.EQ.0 ) GO TO 450
      write(66,4000) ISTOP
      STOP
C
  450 IF (IDATWR.GT.1) GO TO 460
C
C
      write(66,2000) NPAR1
      write(66,2050) NSURF,NSURFP
      write(66,2100) NEGSKS
      write(66,2150) NECON,NEALL
      write(66,2200) NTOUCH,NSNOD
      write(66,2250) NJOIN
      write(66,2300) ICTEMP
      write(66,2310) MTEMP
C
C
  460 IF (JNPORT.EQ.0 .OR. NPUTSV.EQ.0) GO TO 500
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'TYPE-17 ',NSUB,NGALL,(NPAR(I),I=1,20)
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9000) 'TYPE-17 ',NSUB,NGALL,(NPAR(I),I=1,20)
C
 9000 FORMAT ( A,/,(8I10) )
C
C
C
C
C
  500 NFIRST=N6
      IF ( IND.EQ.4 ) NFIRST=N10
C
      NXALL = 3*NSNOD
      NXCON = 3*NTOUCH
      NCEQN = 1
      IF (MODEL.EQ.2) NCEQN = 3
      NCADD = NCEQN*NTOUCH
C
      IF ( IND.NE.4 ) LNWK = 1
      IF ( IND.EQ.4 ) LNWK = NPAR(20)
C
      N110 = NFIRST + 20
      N120 = N110 + 2*NSURFP
      N130 = N120 + NSURF + 1
      N135 = N130 + NSURF + 1
      N136 = N135 + NSURF
      N140 = N136 + NSURF
      N150 = N140 + NSNOD
      N160 = N150 + NXALL
      N170 = N160 + NSNOD
      N180 = N170 + 2
      N190 = N180 + 4*NEALL
      N200 = N190 + NTOUCH + 1
      N210 = N200
      N220 = N210 + NTOUCH
      N230 = N220 + NECON
      N240 = N230 + NTOUCH
      N250 = N240 + NTOUCH
      N260 = N250 + NJOIN*NSNOD
      N270 = N260
      N280 = N270
      N290 = N280 + NCADD
      N300 = N290 + 16
      N500 = N300
      NILAST = N500 - 1
C
      N510 = N500 + 3*NSURFP*ITWO
      N520 = N510 + NXALL*ITWO
      N530 = N520 + NXALL*ITWO
      N540 = N530 + NXCON*ITWO
      N550 = N540 + NXCON*ITWO
      N560 = N550 + NXCON*ITWO
      N570 = N560 + 4*NECON*ITWO
      N580 = N570 + 3*NEALL*ITWO
      N700 = N580
      NRLAST = N700 - 1
C
      N710 = N700 + LNWK*ITWO
      NLAST = N710 - 1
C
C
C         BLANK COMMON
C
      IF (IND.NE.0) GO TO 520
      CALL SIZE (NLAST)
      LINTE =   NILAST - NFIRST + 1
      LREAL = ( NRLAST - N500   + 1 )/ITWO
      CALL LIGHT (RMDUMY,1,1,IA(NFIRST),1,LINTE,0,-100)
      CALL LIGHT (A(N500),1,LREAL,IMDUMY,1,1,0,100)
  520 CALL LIGHT (A(N700),1,LNWK ,IMDUMY,1,1,0,100)
C
      IF ( IND.EQ.4 ) GO TO 570
      M2 = N2
      M3 = N3
      M4 = N4
      M5 = N5
      GO TO 600
C
  570 M2 = N2
      M3 = N7
      M4 = N8
      M5 = N3
C
      IF ( IEQUIT.EQ.1 ) GO TO 600
      M2 = N6
      M5 = N5
C
  600 CALL TOUCH3 ( A(N01),A(N02),A(N04),
     1              A(N1),A(M2),A(M3),A(M4),A(M5),
     1              A(N06),A(N1A),
     1              A(N110),A(N120),A(N130),A(N135),A(N136),A(N140),
     1              A(N150),A(N160),A(N170),A(N180),
     1              A(N190),A(N210),A(N220),
     1              A(N230),A(N240),A(N250),
     1              A(N280),A(N290),
     1              A(N500),A(N510),A(N520),A(N530),
     1              A(N540),A(N550),A(N560),A(N570),
     1              A(N700),
     1              NDOF,NIDM,NSURFP,NECON,NEALL,
     1              NTOUCH,NSNOD,NJOIN)
C
C
C
      IF ( IND.NE.0 ) RETURN
      J = NFIRST - 1
      DO 650 L=1,20
      J=J+1
  650 IA(J) = NPAR(L)
C
      NCEQ = LEADEQ - NEQI
      LWK  = LNWK*ITWO
      LEXTRA = MAX0( NCEQ,LWK )
      MIDEST = NLAST - NFIRST + 1
      NLAST  = NLAST + LEXTRA
      LENGTH = NLAST - NFIRST + 1
      IF (IDATWR.LE.1 ) write(66,2350) LENGTH
      CALL SIZE (NLAST)
C
C
C
      NOLD  = NCEQ  - NCADD
      N280A = N280  - 1
      N700A = N700  - 1
      N700B = N700A + NOLD
      DO 660 L=1,NCADD
  660 IA(N700B+L) = IA(N280A+L)
      IF (NOLD.EQ.0) GO TO 670
      REWIND IT19
      READ  (IT19) ( IA(N700A+L),L=1,NOLD )
  670 REWIND IT19
      WRITE (IT19) ( IA(N700A+L),L=1,NCEQ )
      RETURN
C
C
C
 2000 FORMAT(36H S U R F A C E   D E F I N I T I O N  ///,
     1       14H SURFACE TYPE ,13(2H .),16H( NPAR(1) ). . =,I5/,
     D       32H     EQ.16, 2-D CONTACT SURFACES/,
     E       32H     EQ.17, 3-D CONTACT SURFACES//)
 2050 FORMAT(28H NUMBER OF CONTACT SURFACES ,6(2H .),11H( NPAR(2) ),
     1         5H. . =,I5,
     2      //,24H NUMBER OF SURFACE PAIRS,8(2H .),16H( NPAR(4) ). . = ,
     3         I5 )
 2100 FORMAT(//23H SKEW COORDINATE SYSTEM,
     1      /,5X,19HREFERENCE INDICATOR,8(2H .),16H( NPAR(6) ). . =,I5,
     2      /,5X,23HEQ.0  ALL SURFACE NODES ,
     2      /,11X,37HUSE THE GLOBAL COORDINATE SYSTEM ONLY,
     3      /,5X,25HEQ.1  SURFACE NODES REFER ,
     4      /,11X,25HTO SKEW COORDINATE SYSTEM )
 2150 FORMAT(//40H NUMBER OF CONTACTOR SURFACE SEGMENTS . ,
     1          16H( NPAR(7) ). . =,I5,
     1        //34H TOTAL NUMBER OF SURFACE SEGMENTS ,3(2H .),
     1          16H( NPAR(8) ). . =,I5 )
 2200 FORMAT(//26H NUMBER OF CONTACTOR NODES,7(2H .),13H( NPAR(9) ). ,
     1          3H. =,I5,
     2       //30H TOTAL NUMBER OF SURFACE NODES,5(2H .),
     2          16H( NPAR(10) ) . =,I5)
 2250 FORMAT(//36H MAXIMUM NUMBER OF SURFACE SEGMENTS ,2(2H .),
     1          16H( NPAR(11) ) . =,I5,/,
     1          38H SHARING A COMMON CONTACT SURFACE NODE )
 2300 FORMAT(/40H PRINT CONTROL FOR CONTACT SURFACES . . ,
     1       16H( NPAR(13) ) . =,I5/
     2       5X,43HEQ.0  PRINT CONSISTENT CONTACT NODAL FORCES,/,
     1    5X,47HEQ.1  PRINT AVERAGE CONTACTOR SEGMENT TRACTIONS ,/,
     1    5X,43HEQ.2  PRINT CONSISTENT CONTACT NODAL FORCES,/,
     1   11X,39HAND AVERAGE CONTACTOR SEGMENT TRACTIONS )
 2310 FORMAT(//16H FRICTION MODEL ,12(2H .),16H( NPAR(15) ) . =,I5,/,
     1    5X,26HEQ.0, FRICTIONLESS CONTACT,/,
     1    5X,23HEQ.1, INFINITE FRICTION ,/)
 2350 FORMAT(//43H LENGTH OF ARRAY NEEDED FOR STORING SURFACE/
     1       43H DATA. . . . . . . . . . . . .(MIDEST). . =,I5/)
C
C
C
 2500 FORMAT(1H1,38HERROR IN CONTACT SURFACE CONTROL CARDS/
     1       16H SURFACE GROUP =,I5/)
 2550 FORMAT(I5,31H. FULL NEWTON ITERATION MUST BE  ,/,
     1       7X,34HEMPLOYED WHEN CONTACT SURFACES ARE  ,/,
     1       7X,40HPRESENT.  (SEE MASTER CONTROL CARD NO.7)   )
 2560 FORMAT(I5,35H. LINE SEARCH SOLUTION MUST NOT BE  ,/,
     1       7X,34HEMPLOYED WHEN CONTACT SURFACES ARE   ,/,
     1       7X,40HPRESENT.  (SEE MASTER CONTROL CARD NO.7)  )
 2600 FORMAT(I5,41H. AT LEAST TWO SURFACES MUST BE SPECIFIED/
     1          30X,7HNSURF =,I2)
 2700 FORMAT(I5,38H. NPAR(4) MUST BE .GE. 1 ... NPAR(4) =,I2)
 2800 FORMAT(I5,10H. NSKEWS =,I5,14H AND NPAR(6) =,I5,
     1       7X,19H ARE NOT COMPATIBLE )
 2900 FORMAT(I5,48H. AT LEAST ONE CONTACTOR SEGMENT MUST BE PRESENT/
     1           30X,7HNECON =,I5,/)
 3000 FORMAT(I5,53H. AT LEAST TWO CONTACT SURFACE SEGMENTS MUST BE INPUT
     1        /,30X,7HNEALL =,I5,/)
 3210 FORMAT(I5,23H. NPAR(9) MUST BE GE. 3 ,/,
     1       7X,14HINPUT NPAR(9)=,I5 )
 3220 FORMAT (I5,23H. NPAR(10) MUST BE GE.6 ,/,
     1       7X,15HINPUT NPAR(10)=,I5 )
 3250 FORMAT(I5,48H. AT LEAST ONE SEGMENT MUST JOIN AT EACH CONTACT,/,
     1       8X,12HSURFACE NODE,/,
     1      30X,10HNPAR(11) =,I5,/)
 3300 FORMAT(I5,40H. NPAR(13) MUST BE .LE. 2 ... NPAR(13) =,I2)
 3310 FORMAT(I5,31H. NPAR(15) MUST BE EQ.0 OR EQ.1 ,/,
     1        7X,15HINPUT NPAR(15)=,I5 )
 4000 FORMAT(//25H TOTAL NUMBER OF ERRORS =,I5////
     1       34H S T O P  (ERRORS IN SURFACE DATA) )
C
C
      END
      SUBROUTINE TOUCH3 (NID,IDI,NOD,
     1                   ID,X,Y,Z,HT,
     1                   RSDCOS,NODSYS,
     1                   ISURFP,IFSN,IFSE,IPS,ISV,NODSF,
     1                   LMS,ISKEW,INUMEQ,NCA,
     1                   LMAXA,INODE,ISECT,
     1                   ITS,ITSP,JOIN,
     1                   LMHT,IDBUG,
     1                   FCOFF,XYZ,XYZS,CPROLD,
     1                   CPR,DELTA,GUSTAV,VN,
     1                   T,
     1                   NDOF,NIDM,NSURFP,NECON,NEALL,
     1                   NTOUCH,NSNOD,NJOIN)
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 /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1           ,NDFD,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /DIM/ N0,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14,N15
      COMMON /LOA/ NLOAD,NPR2,NPR3,NPBM,NP3DB,NPPL,NPSH,NODE3,IDGRAV,
     1             NPDIS,NTEMP,IDCFGL,ISMASS
      COMMON /ELSTP/ TIME,IDTHF
      COMMON /VAR/ NGALL,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 /FREQIF/ ISTOH,N1A,N1B,N1C,N1S
      COMMON /RANDI/ N0A,N1D,IELCPL
      COMMON /ADDB/ NEQL,NEQR,MLA,NBLOCK
      COMMON /MDFRDM/ IDOF(12)
      COMMON /SKEW/ NSKEWS
      COMMON /DISCON/ NDISCE,MIDM
      COMMON /CNTACT/ NEQI,LEADEQ,NCE2D,NCE3D
      COMMON /ITRATE/ RCENRM,RCNORM,RCONSM,RCTOL
      COMMON /SURF3/ IS,KFS,KLS,JFS,JLS
      COMMON /TRGT3/ IT,KFT,KLT,JFT,JLT
      COMMON /MATCH3/ ISR,IPAIR,JTSECT
      COMMON /ISET3/ ISET,NCEQN
      COMMON /VALUES/ C1,C2,C3,C4,C5
      COMMON /TAPES/ IIN,IOUT
      COMMON /DPR/ ITWO
      COMMON A(1)
      REAL A
      DIMENSION IA(1)
C
      DIMENSION NID(*),IDI(NIDM,*),NOD(*),
     1          ID(NDOF,*),X(*),Y(*),Z(*),HT(*),
     1          RSDCOS(9,*),NODSYS(*)
      DIMENSION ISURFP(2,NSURFP),IFSN(*),IFSE(*),IPS(*),ISV(*),NODSF(*),
     1          LMS(3,NSNOD),ISKEW(*),INUMEQ(2),NCA(4,NEALL),
     1          LMAXA(*),INODE(*),ISECT(*),
     1          ITS(*),ITSP(*),JOIN(NJOIN,NSNOD),
     1          LMHT(*),IDBUG(16)
      DIMENSION FCOFF(3,NSURFP),XYZ(3,NSNOD),
     1          XYZS(3,NSNOD),CPROLD(3,NTOUCH),
     1          CPR(3,NTOUCH),DELTA(3,NTOUCH),GUSTAV(4,NECON),
     1          VN(3,NEALL),T(*)
C
      DIMENSION IMDUMY(1,1),RMDUMY(1,1)
C
      EQUIVALENCE ( A(1),IA(1) )
      EQUIVALENCE (  NPAR(2),NSURF  ),( NPAR(6),NEGSKS ),
     1            ( NPAR(13),ICPRNT ),( NPAR(15),MODEL ),
     1            ( NPAR(16),IDBPR  ),( NPAR(20),LNWK  )
C
      C1 = 0.3868945341D0
      C2 = 0.1036680780D0
      C3 = 0.0277777778D0
      C4 = 0.0074430331D0
      C5 = 0.0019943547D0
C
      IF (IND.GT.0) GO TO 600
      IF (IDATWR.GT.1) GO TO 10
      write(66,2000)
      write(66,2050)
C
C
C
   10 ISTOP=0
      IFSN(1)=1
      IFSE(1)=1
      NKSUM=0
      NJSUM=0
      NG   =NGALL-NEGL-NEGNL
C
      DO 20 I=1,NSURF
C
      READ (IIN,1000) ISURF,NSEG,NNODE,IPRINT,IPLOT
C
      NKSUM = NKSUM + NNODE
      NJSUM = NJSUM + NSEG
      IFSN(I+1) = IFSN(I) + NNODE
      IFSE(I+1) = IFSE(I) + NSEG
C
      IF (ISURF.EQ.I) GO TO 30
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,2550) ISTOP,I,ISURF
C
   30 IF (NSEG.GE.1) GO TO 40
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,2600) ISTOP,ISURF,NSEG
C
   40 IF (NNODE.GE.3) GO TO 50
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,2650) ISTOP,ISURF,NNODE
C
   50 IPS(I)=IPRINT
      IF (IPRINT.EQ.0 .OR. IPRINT.EQ.1) GO TO 55
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG
      write(66,2700) ISTOP,ISURF,IPRINT
C
   55 ISV(I)=IPLOT
      IF (IPLOT.EQ.0 .OR. IPLOT.EQ.1) GO TO 60
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,2500) NG