a36.for

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

C
      DIMENSION X(*),R(*),
     1          RSDCOS(9,*)
      DIMENSION ISURFP(2,NSURFP),IFSN(*),LMS(2,NSNOD),
     1          INODE(*),ISEG(*),IDBUG(20),
     1          ISV(*),
     1          INUMEQ(2),ISKEW(*)
      DIMENSION XYZ(2,NSNOD),CPROLD(2,NTOUCH),
     1          XYZS(2,NSNOD),XLN(*),T(2,NSEG),
     1          CPR(2,NTOUCH),GAUSS(2,NSNOD),CLOSED(*)
      DIMENSION RA(4)
      DIMENSION RMDUMY(1,1),IMDUMY(1,1)
      EQUIVALENCE (NPAR(5),IAXIS),(NPAR(6),NEGSKS)
C
      IDBA=IDBUG(5)
      IDBB=IDBUG(6)
      IF (IDBA.LT.3) GO TO 50
      CALL EXAMIN (XYZ,2,NSNOD,IMDUMY,1,1,3,2)
      CALL EXAMIN (X,1,NEQ,IMDUMY,1,1,3,10)
      CALL EXAMIN (RMDUMY,1,1,INODE,1,NTOUCH,3,-5)
      CALL EXAMIN (RMDUMY,1,1,ISEG,1,NSEG,3,-6)
C
C
C
   50 DO 60 K=1,NSNOD
      DO 60 L=1,2
   60 XYZS(L,K)=0.0D0
      DO 70 K=1,NTOUCH
      CPR(1,K)=0.0D0
   70 CPR(2,K)=0.0D0
      DO 80 J=1,NSEG
      XLN(J)=0.0D0
      T(1,J)=0.0D0
   80 T(2,J)=0.0D0
C
      NLOWER = INUMEQ(1) + 1
      NUPPER = INUMEQ(2)
      DO 90 L=NLOWER,NUPPER
   90 X(L)=0.0D0
C
C
C
C
      DO 110 K=1,NSNOD
      DO 110 L=1,2
      KDOF=LMS(L,K)
      IF (KDOF.GT.NEQ) KDOF = KDOF - NEQ
      IF (KDOF.LT. 0 ) KDOF = NEQ  - KDOF
      XMOVE=0.0D0
      IF (KDOF.GT.0) XMOVE=X(KDOF)
      XYZS(L,K)=XMOVE
  110 CONTINUE
      IF (NEGSKS.GT.0) CALL DIRCOS (RSDCOS,XYZS,ISKEW,NSNOD,2,1)
C
      DO 130 K=1,NSNOD
      DO 130 L=1,2
      XYZS(L,K) = XYZ(L,K) + XYZS(L,K)
  130 CONTINUE
C
      DO 150 I=1,NSURF
      IS=I
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,1)
C
      DO 170 J=JFS,JLS
      IF (ISEG(J).EQ.0) GO TO 170
      NA=J-IS
      NB=NA+1
      YDIST=XYZS(1,NB) - XYZS(1,NA)
      ZDIST=XYZS(2,NB) - XYZS(2,NA)
      XLEN=SQRT( YDIST*YDIST + ZDIST*ZDIST )
      IF (XLEN.LE.0.0D0) GO TO 180
      XLN(J)=XLEN
      T(1,J)=YDIST/XLEN
      T(2,J)=ZDIST/XLEN
  170 CONTINUE
C
      IF (KES.EQ.KLS) GO TO 200
C
      XLN(JFSD)=XLN(JLS)
      T(1,JFSD)=T(1,JLS)
      T(2,JFSD)=T(2,JLS)
C
      XLN(JLSD)=XLN(JFS)
      T(1,JLSD)=T(1,JFS)
      T(2,JLSD)=T(2,JFS)
  200 CONTINUE
  150 CONTINUE
      GO TO 190
  180 NAL = NA - KFS + 1
      NBL = NAL + 1
      write(66,2000) IS,NAL,NBL,XLEN
      STOP
C
  190 IF (IDBA.GE.1) CALL EXAMIN (XYZS,2,NSNOD,IMDUMY,1,1,3,3)
      IF (IDBA.GE.2) CALL EXAMIN (XLN,1,NSEG,IMDUMY,1,1,3,6)
      IF (IDBA.GE.2) CALL EXAMIN (T,2,NSEG,IMDUMY,1,1,3,7)
      IF (IDBB.GE.3) CALL EXAMIN (R,1,NEQ,IMDUMY,1,1,3,11)
C
C
C
C
      IF (KPRI.GT.0) GO TO 250
      DO 210 K=1,NTOUCH
      CPR(1,K) = CPROLD(1,K)
  210 CPR(2,K) = CPROLD(2,K)
      IF (NEGSKS.GT.0) CALL DIRCOS (RSDCOS,CPR,ISKEW,NTOUCH,2,1)
      GO TO 600
C
  250 CSMALL = 0.0D0
      CTOTAL = 0.0D0
C
      DO 260 K=1,NTOUCH
      FY=0.0D0
      FZ=0.0D0
      IF (INODE(K).EQ.1) GO TO 270
      IYDOF=LMS(1,K)
      IZDOF=LMS(2,K)
      IF (IYDOF.LE.NEQI .AND. IYDOF.GT.0) FY = -R(IYDOF)
      IF (IZDOF.LE.NEQI .AND. IZDOF.GT.0) FZ = -R(IZDOF)
      CPR(1,K)= FY
      CPR(2,K)= FZ
C
  270 CTOTAL = CTOTAL + FY*FY + FZ*FZ
      DIFY = CPROLD(1,K) - FY
      DIFZ = CPROLD(2,K) - FZ
      CSMALL = CSMALL + DIFY*DIFY + DIFZ*DIFZ
      CPROLD(1,K) = FY
      CPROLD(2,K) = FZ
  260 CONTINUE
C
      IF (NEGSKS.GT.0) CALL DIRCOS (RSDCOS,CPR,ISKEW,NTOUCH,2,1)
C
  280 RCN = RCNORM*RCNORM + CTOTAL
      RCE = RCENRM*RCENRM + CSMALL
      IF (RCN.GT.0.0D0) RCNORM = SQRT(RCN)
      IF (RCE.GT.0.0D0) RCENRM = SQRT(RCE)
C
      IF (IDBB.GE.1) CALL EXAMIN (CPR,2,NTOUCH,IMDUMY,1,1,3,8)
      IF (ISET.LE.1) RETURN
C
C
C .                                                                   .
C .                                                                   .
C
  600 DO 550 I=1,NSURF
C
      IS=I
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,1)
      IF (ISR.EQ.0) GO TO 550
C
      DO 630 M=1,2
      DO 640 K=KFS,KES
      CLOSED(K)=0.0D0
      DO 640 L=1,2
  640 GAUSS(L,K)=0.0D0
C
      GEC=0.0D0
      R1=XYZS(1,KES)
      DO 500 K=KFS,KES
      NSA=K+IS-1
      NSB=K+IS
      JSA=ISEG(NSA)
      JSB=ISEG(NSB)
      LA=1
      LB=1
      FIXED=1.0D0
      IF (JSA.EQ.0 .OR. JSA.EQ.1) LA=0
      IF (JSB.EQ.0 .OR. JSB.EQ.1) LB=0
      IF (LA.EQ.0 .AND. LB.EQ.0) GO TO 510
C
C
C
      DO 570 L=1,4
  570 RA(L)=1.0D0
      IF (IAXIS.GT.0) GO TO 580
      IF (K.EQ.KFS) GO TO 590
      R1=XYZS(1,K-1)
  590 R2=XYZS(1,K)
      RA(1) = (R1+R2)/2.0
      RA(2) = (3.0*R2+R1)/4.0
      IF (K.EQ.KLS) GO TO 580
      R1=XYZS(1,K)
      R2=XYZS(1,K+1)
      RA(3) = (3.0*R1+R2)/4.0
      RA(4) = (R1+R2)/2.0
  580 CONTINUE
C
C
C
      DSA=XLN(NSA)*FLOAT(LA)/6.0
      DSB=XLN(NSB)*FLOAT(LB)/6.0
      GTEMP = 2.0*( RA(2)*DSA + RA(3)*DSB ) - RA(1)*DSA*GEC
      GAUSS(1,K) = GTEMP*FIXED
      GEC = RA(4)*DSB/GAUSS(1,K)
      GO TO 520
  510 GAUSS(1,K)=FIXED
      GEC=0.0D0
  520 GAUSS(2,K)=GEC
  500 CONTINUE
C
      IF (KES.EQ.KLS) GO TO 610
      CA = -RA(4)*DSB
      CB = 1.0D0
      KESDO=KES-1
      DO 620 K=KFS,KESDO
      CA = -CA*CB
      CB = GAUSS(2,K)
      CLOSED(K) = CA/GAUSS(1,K)
      GAUSS(1,KES) = GAUSS(1,KES) - CA*CLOSED(K)
  620 CONTINUE
C
      GAUSS(1,KES) = GAUSS(1,KES) - CA*GAUSS(2,KES-1)
     1                            - RA(1)*DSA*CLOSED(KES-1)
C
  610 DO 530 K=KFS,KES
      IF (K.GT.KFS) CPR(M,K) = CPR(M,K) - GAUSS(2,K-1)*CPR(M,K-1)
      IF (KES.EQ.KLS) GO TO 530
      IF (K.LT.KES) CPR(M,KES) = CPR(M,KES) - CLOSED(K)*CPR(M,K)
  530 CONTINUE
C
C
C
      CPR(M,KES)=CPR(M,KES)/GAUSS(1,KES)
C
      KDUMY=KES-1
      KTEMP=KDUMY+KFS
      DO 540 KK=KFS,KDUMY
      K=KTEMP-KK
      CPR(M,K) = CPR(M,K)/GAUSS(1,K) - GAUSS(2,K)*CPR(M,K+1)
      IF (KES.LT.KLS) CPR(M,K) = CPR(M,K) - CLOSED(K)*CPR(M,KES)
  540 CONTINUE
C
      IF (KES.LT.KLS) CPR(M,KLS)=CPR(M,KFS)
  630 CONTINUE
  550 CONTINUE
      IF (IDBB.GE.2) CALL EXAMIN (CPR,2,NTOUCH,IMDUMY,1,1,3,12)
C
 2000 FORMAT(//,5X,19HFOR SURFACE NUMBER=,I5,/,
     1          5X,30HSEGMENT BOUNDED BY LOCAL NODES,I5,3X,3HAND,I5,/,
     1          5X,11HHAS LENGTH=,E15.4,/,
     1          5X,16HSTOPPED IN YZNEW ,/)
C
      END
      SUBROUTINE FRICT2 (ISURFP,IFSN,NODSF,
     1                   INODE,ISEG,ITS,IDBUG,
     1                   IPS,ISV,
     1                   JSLIDE,KSLIDE,
     1                   FCOFF,
     1                   XYZS,XLN,T,
     1                   CPR,BETA,CFORCE,
     1                   NSURF,NSURFP,NTOUCH,NSNOD,NSEG)
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            ,NDFD,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /CNTACT/ NEQI,LEADEQ,NCE2D,NCE3D
      COMMON /SRFACE/ IS,KFS,KLS,KES,JFS,JLS,JFSD,JLSD
      COMMON /MATCH/ ISR,IPAIR,JFSEG,JLSEG,JTSEG
      COMMON /COSINE/ NODEA,NODEB,T1,T2,XLEN,JSS
      COMMON /SPRESS/ PNA,PNB,PTA,PTB,APTAV,PTAV,PNAV
      COMMON /SFORCE/ RTALL,RNALL,FY,FZ
C
      DIMENSION ISURFP(2,NSURFP),IFSN(*),NODSF(*),
     1          INODE(*),ISEG(*),ITS(*),IDBUG(20),
     1          IPS(*),ISV(*),
     1          JSLIDE(*),KSLIDE(*)
      DIMENSION FCOFF(3,NSURFP),
     1          XYZS(2,NSNOD),XLN(*),T(2,NSEG),
     1          CPR(2,NTOUCH),BETA(*),CFORCE(2,NSNOD)
      DIMENSION RA(8)
      DIMENSION RMDUMY(1,1),IMDUMY(1,1)
      EQUIVALENCE (NPAR(5),IAXIS),(NPAR(13),ICPRNT),
     1            (NPAR(15),MODEL)
C
      DATA LONELY / 10 /
C
C
C
C
C
      ZERO = 0.D0
      ONE  = 1.D0
      IDBA=IDBUG(7)
      IDBB=IDBUG(8)
      KCPRI=1
      IF (IDBA.LE.2) GO TO 50
      CALL EXAMIN (RMDUMY,1,1,INODE,1,NTOUCH,4,-5)
      CALL EXAMIN (RMDUMY,1,1,ISEG,1,NSEG,4,-6)
      CALL EXAMIN (CPR,2,NTOUCH,IMDUMY,1,1,4,12)
   50 IF (KPRI.EQ.0 .OR. IDBB.GE.1) KCPRI=0
      IF (KPRI.EQ.0) CALL PRINT2 (IPS,ISV,0)
C
      DO 60 K=1,NSNOD
      CFORCE(1,K)=0.0D0
   60 CFORCE(2,K)=0.0D0
C
      DO 700 I=1,NSURF
      IS=I
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,1)
      IF (ISR.EQ.0) GO TO 700
      IF (KCPRI.EQ.0) CALL PRINT2 (IPS,ISV,1)
C
      DO 250 J=JFS,JLS
C
      JSS = ISEG(J)
      JSSA = IABS( JSS )
      LAST = JSLIDE(J)
      NA=J-IS
      NB=NA+1
      NODEA=NODSF(NA)
      NODEB=NODSF(NB)
      T1=T(1,J)
      T2=T(2,J)
      XLEN=XLN(J)
C
      PNA=0.0D0
      PNB=0.0D0
      PTA=0.0D0
      PTB=0.0D0
      PNAV=0.0D0
      PTAV=0.0D0
      APTAV=0.0D0
      RNALL=0.0D0
      RTALL=0.0D0
      IF (JSSA.LT.2) GO TO 200
C
C
C
      IF ( IAXIS.EQ.0 ) GO TO 110
      R1=1.0D0
      R2=1.0D0
      RA(1)=1.0D0
      RA(2)=1.0D0
      RA(3)=1.0D0
      RA(4)=0.0D0
      RA(5)=2.0D0
      RA(6)=9.0D0
      RA(7)=3.0D0
      RA(8)=3.0D0
      GO TO 130
  110 R1=XYZS(1,NA)
      R2=XYZS(1,NB)
      RA(1) = ( 3.0*R1 + R2 )/4.0
      RA(2) = ( R1 + R2 )/2.0
      RA(3) = ( 3.0*R2 + R1 )/4.0
      RA(4)=R1-R2
      RA(5)=R1+R2
      RA(6)=7.0*R1+2.0*R2
      RA(7)=2.0*R1+R2
      RA(8)=2.0*R2+R1
  130 CONTINUE
C
C
C
      PNA = -T2*CPR(1,NA) + T1*CPR(2,NA)
      PNB = -T2*CPR(1,NB) + T1*CPR(2,NB)
      PTA =  T1*CPR(1,NA) + T2*CPR(2,NA)
      PTB =  T1*CPR(1,NB) + T2*CPR(2,NB)
C
      PNAV = ( RA(7)*PNA + RA(8)*PNB )/( 3.0*RA(5) )
      PTAV = ( RA(7)*PTA + RA(8)*PTB )/( 3.0*RA(5) )
      APNAV = MAX(PNAV,ZERO)
      APTAV = ABS( PTAV )
      PSIGN = SIGN( ONE,PTAV )
C
      RNA = XLEN*( 2.0*RA(1)*PNA + RA(2)*PNB )/6.0
      RNB = XLEN*( 2.0*RA(3)*PNB + RA(2)*PNA )/6.0
      RTA = XLEN*( 2.0*RA(1)*PTA + RA(2)*PTB )/6.0
      RTB = XLEN*( 2.0*RA(3)*PTB + RA(2)*PTA )/6.0
C
      RNALL = RNA + RNB
      RTALL = RTA + RTB
C
C
C
      CHECK = PNAV
      IF (CHECK.GE.0.0D0) GO TO 150
      PNAV=0.0D0
      PTAV=0.0D0
      RNALL=0.0D0
      RTALL=0.0D0
      ISEG(J)=-1
      GO TO 200
C
C
C
  150 JTSEG=ITS(NA)
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,2)
      SCF=FCOFF(1,IPAIR)
      DCF=FCOFF(2,IPAIR)
      PRCENT=FCOFF(3,IPAIR)
C
      PCF = SCF
      IF (JSSA.EQ.3) PCF=DCF
C
      PNAV = APNAV
      PFAV = PCF*PNAV
      CHECK = PFAV - APTAV
      IF (JSSA.EQ.2 .AND. CHECK.GE.0.0D0) GO TO 160
      CHECK = ( 1.0 - PRCENT )*PFAV - APTAV
      IF (JSSA.EQ.3 .AND. CHECK.GE.0.0D0) GO TO 160
C
C
C
      ISEG(J)=3
      PFAV = DCF*APNAV
      PTAV = MIN(PFAV,APTAV)
      PTAV = PSIGN*PTAV
      DUMY = XLEN*PTAV/6.0
      RTA = RA(7)*DUMY
      RTB = RA(8)*DUMY
      RTALL = RTA + RTB
      GO TO 170
C
  160 ISEG(J)=2
      IF (LAST.EQ.1) ISEG(J)=-2
C
C
C
  170 CFORCE(1,NA) = CFORCE(1,NA) + RTA*T1 - RNA*T2
      CFORCE(2,NA) = CFORCE(2,NA) + RTA*T2 + RNA*T1
      CFORCE(1,NB) = CFORCE(1,NB) + RTB*T1 - RNB*T2
      CFORCE(2,NB) = CFORCE(2,NB) + RTB*T2 + RNB*T1
C
  200 JSS = ISEG(J)
      IF (KCPRI.EQ.0) CALL PRINT2 (IPS,ISV,2)
C
  250 CONTINUE
C
C
C
      IF (KES.EQ.KLS) GO TO 260
      DO 270 L=1,2
      CFORCE(L,KFS) = CFORCE(L,KFS) + CFORCE(L,KLS)
  270 CFORCE(L,KLS)=0.0D0
      ISEG(JFSD)=ISEG(JLS)
      ISEG(JLSD)=ISEG(JFS)
      JSLIDE(JFSD)=JSLIDE(JLS)
      JSLIDE(JLSD)=JSLIDE(JFS)
  260 CONTINUE
C
C
C
C
C
      NOSOL = 0
      ITITLE = 1
      IF (KCPRI.EQ.0) ITITLE = 0
      NTOTAL = KES - KFS + 1
C
      DO 300 K=KFS,KES
      FY=0.0D0
      FZ=0.0D0
      NODEA=NODSF(K)
      KSN = INODE(K)
      KSNA = IABS(KSN)
      LAST = KSLIDE(K)
      IF (KSN.EQ.1) GO TO 310
      NSA=K+IS-1
      NSB=K+IS