a36.for

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

      JSA=ISEG(NSA)
      JSB=ISEG(NSB)
      IF (JSA.NE.1 .AND. JSA.NE.0) GO TO 310
      IF (JSB.NE.1 .AND. JSB.NE.0) GO TO 310
C
      IF (ITITLE.EQ.0) CALL PRINT2 (IPS,ISV,3)
      ITITLE=1
C
      TY = T(1,NSA) + T(1,NSB)
      TZ = T(2,NSA) + T(2,NSB)
      RESULT = SQRT(TY*TY + TZ*TZ)
      TY=TY/RESULT
      TZ=TZ/RESULT
      FN =-TZ*CPR(1,K) + TY*CPR(2,K)
      FT = TY*CPR(1,K) + TZ*CPR(2,K)
      AFN = MAX(FN,ZERO)
      AFT = ABS( FT )
      FSIGN = SIGN( ONE,FT )
C
      CHECK = FN
      IF (CHECK.GE.0.0D0) GO TO 330
      CFORCE(1,K)=0.0D0
      CFORCE(2,K)=0.0D0
      INODE(K) = -LONELY
      GO TO 320
C
  330 JTSEG=ITS(K)
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,2)
      SCF = FCOFF(1,IPAIR)
      DCF = FCOFF(2,IPAIR)
      PRCENT = FCOFF(3,IPAIR)
C
      PCF = SCF
      IF (KSNA.EQ.3) PCF=DCF
C
      RESIST = PCF*AFN
      CHECK = RESIST - AFT
      IF (KSNA.EQ.2 .AND. CHECK.GE.0.0D0) GO TO 340
      CHECK = ( 1.0 - PRCENT )*RESIST - AFT
      IF (KSNA.EQ.3 .AND. CHECK.GE.0.0D0) GO TO 340
C
      INODE(K)=3*LONELY
      FR = DCF*AFN
      FT = MIN(FR,AFT)
      FT = FSIGN*FT
      GO TO 350
  340 INODE(K) = 2*LONELY
      IF (LAST.EQ.1) INODE(K) = -2*LONELY
C
  350 FY = (TY*FT - TZ*FN)
      FZ = (TZ*FT + TY*FN)
      CFORCE(1,K) = FY
      CFORCE(2,K) = FZ
  320 IF (KCPRI.EQ.0) CALL PRINT2 (IPS,ISV,4)
      GO TO 300
C
  310 NOSOL = NOSOL + 1
C
  300 CONTINUE
      IF (NOSOL.LT.NTOTAL) GO TO 360
      IF (KCPRI.EQ.0) CALL PRINT2 (IPS,ISV,5)
C
  360 IF (KES.EQ.KLS) GO TO 370
C
      IF (INODE(KFS).EQ.1) GO TO 370
      JSA = ISEG(JFSD)
      JSB = ISEG(JFS)
      IF (JSA.NE.1 .AND. JSA.NE.0) GO TO 370
      IF (JSB.NE.1 .AND. JSB.NE.0) GO TO 370
      INODE(KLS)=INODE(KFS)
C
  370 IF (KPRI.EQ.0) GO TO 700
C
C
C
  600 JOLD=ISEG(JFSD)
      JOLDA=IABS(JOLD)
C
      DO 610 K=KFS,KLS
      JNEW=ISEG(K+IS)
      JNEWA=IABS(JNEW)
      IF (IABS(INODE(K)).LT.LONELY) GO TO 620
      INODE(K)=INODE(K)/LONELY
      GO TO 650
  620 INODE(K)=1
      IF (JOLD.EQ.-2 .OR. JNEW.EQ.-2) GO TO 630
      IF (JOLD.EQ.2 .OR. JNEW.EQ.2) GO TO 630
      IF (JOLDA.EQ.3 .OR. JNEWA.EQ.3) GO TO 640
      IF (JOLD.EQ.-1 .OR. JNEW.EQ.-1) INODE(K)=-1
      GO TO 650
  630 INODE(K)=2
      GO TO 650
  640 INODE(K)=3
  650 JOLD=JNEW
      JOLDA=JNEWA
  610 CONTINUE
  700 CONTINUE
C
      DO 660 J=1,NSEG
      IF (ISEG(J).EQ.-2) ISEG(J)=3
      IF (IABS(ISEG(J)).EQ.3) JSLIDE(J)=1
  660 CONTINUE
      DO 670 K=1,NTOUCH
      IF (INODE(K).EQ.-2) INODE(K)=3
      IF (IABS(INODE(K)).EQ.3) KSLIDE(K)=1
  670 CONTINUE
C
      DO 675 I=1,NSURF
      IS=I
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,1)
      IF (ISR.EQ.0) GO TO 675
      IF (KES.LT.KLS) KSLIDE(KLS) = KSLIDE(KFS)
  675 CONTINUE
C
C
C
      DO 680 K=1,NTOUCH
      CPR(1,K) = CFORCE(1,K)
      CPR(2,K) = CFORCE(2,K)
      CFORCE(1,K) = 0.0D0
      CFORCE(2,K) = 0.0D0
  680 CONTINUE
C
      IF (IDBA.LT.1) GO TO 685
      CALL EXAMIN (CPR,2,NTOUCH,IMDUMY,1,1,4,13)
      CALL EXAMIN (RMDUMY,1,1,ISEG,1,NSEG,4,-16)
      CALL EXAMIN (RMDUMY,1,1,INODE,1,NTOUCH,4,-17)
  685 IF (IDBA.LT.2) GO TO 690
      CALL EXAMIN (RMDUMY,1,1,JSLIDE,1,NSEG,4,-12)
      CALL EXAMIN (RMDUMY,1,1,KSLIDE,1,NTOUCH,4,-13)
C
C
C
  690 IF (KPRI.GT.0 .OR. ICPRNT.EQ.2) RETURN
C
      DO 710 K=1,NSNOD
      CFORCE(1,K)=0.0D0
  710 CFORCE(2,K)=0.0D0
C
      DO 720 I=1,NSURF
      IS=I
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,1)
      IF (ISR.EQ.0) GO TO 720
      DO 730 K=KFS,KES
      JTSEG=ITS(K)
      IF (JTSEG.EQ.0) GO TO 730
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,2)
      IT=ISURFP(1,IPAIR)
      NA=JTSEG-IT
      NB=NA+1
      BT=BETA(K)
      DO 740 L=1,2
      CPRCON = CPR(L,K)
      CFORCE(L,K) = CFORCE(L,K) + CPRCON
      CFORCE(L,NA) = CFORCE(L,NA) - (1.0-BT)*CPRCON
  740 CFORCE(L,NB) = CFORCE(L,NB) - BT*CPRCON
  730 CONTINUE
  720 CONTINUE
C
      DO 750 I=1,NSURF
      IS=I
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,1)
      IF (KES.EQ.KLS) GO TO 750
      DO 760 L=1,2
      CFORCE(L,KFS) = CFORCE(L,KFS) + CFORCE(L,KLS)
  760 CFORCE(L,KLS) = 0.0D0
  750 CONTINUE
C
      DO 770 I=1,NSURF
      IS=I
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,1)
      CALL PRINT2 (IPS,ISV,6)
      DO 770 K=KFS,KES
      NA=K
      NODEA=NODSF(NA)
      FY=CFORCE(1,NA)
      FZ=CFORCE(2,NA)
      CFORCE(1,NA)=0.0D0
      CFORCE(2,NA)=0.0D0
      CALL PRINT2 (IPS,ISV,7)
  770 CONTINUE
      RETURN
      END
      SUBROUTINE TRGET2 ( ISURFP,IFSN,ISEG,NSURF,NSURFP,IFLAG)
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /SRFACE/ IC,KFC,KLC,KEC,JFC,JLC,JFCD,JLCD
      COMMON /PAIR/ IT,KFT,KLT,KET,JFT,JLT,JFTD,JLTD
      COMMON /MATCH/ ISR,IPAIR,JFSEG,JLSEG,JTSEG
C
      DIMENSION ISURFP(2,NSURFP),IFSN(*),ISEG(*)
C
      GO TO ( 50,80,200 ) IFLAG
C
C
C
   50 ISR=0
      DO 60 L=1,NSURFP
      IF (ISURFP(2,L).NE.IC) GO TO 60
      ISR=IC
      GO TO 250
   60 CONTINUE
      GO TO 250
C
C
C
   80 DO 100 I=1,NSURF
      IF (I.EQ.IC) GO TO 100
      ITRGET=I
      JFSEG=IFSN(I)+I
      JLSEG=IFSN(I+1)+I-2
      IF (JFSEG.LE.JTSEG .AND. JTSEG.LE.JLSEG) GO TO 110
  100 CONTINUE
C
C
C
  110 DO 150 L=1,NSURFP
      IPAIR=L
      IF (ISURFP(1,L).EQ.ITRGET .AND. ISURFP(2,L).EQ.IC) GO TO 160
  150 CONTINUE
  160 RETURN
C
  200 KFT=IFSN(IT)
      KLT=IFSN(IT+1)-1
      JFT=KFT+IT
      JLT=KLT+IT-1
      JFTD=JFT-1
      JLTD=JLT+1
      KET=KLT
      IF (ISEG(JFTD).NE.0) KET=KLT-1
  250 KFC=IFSN(IC)
      KLC=IFSN(IC+1)-1
      JFC=KFC+IC
      JLC=KLC+IC-1
      JFCD=JFC-1
      JLCD=JLC+1
      KEC=KLC
      IF (ISEG(JFCD).NE.0) KEC=KLC-1
      RETURN
C
      END
      SUBROUTINE ADJOIN (ISURFP,IFSN,NODSF,
     1                   INODE,ISEG,ITS,IDBUG,
     1                   XYZ,
     1                   XYZS,XLN,T,
     1                   BETA,DELTA,
     1                   NSURF,NSURFP,NTOUCH,NSNOD,NSEG)
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 /CNTACT/ NEQI,LEADEQ,NCE2D,NCE3D
      COMMON /SRFACE/ IC,KFC,KLC,KEC,JFC,JLC,JFCD,JLCD
      COMMON /MATCH/ ISR,IPAIR,JFSEG,JLSEG,JTSEG
      COMMON /PAIR/ IT,KFT,KLT,KET,JFT,JLT,JFTD,JLTD
C
      DIMENSION ISURFP(2,NSURFP),IFSN(*),NODSF(*),
     1          INODE(*),ISEG(*),ITS(*),IDBUG(20)
      DIMENSION XYZ(2,NSNOD),
     1          XYZS(2,NSNOD),XLN(*),T(2,NSEG),
     1          BETA(*),DELTA(2,NSNOD)
      DIMENSION XNC(2),XNT(2)
      DIMENSION DIST(2)
      DIMENSION RMDUMY(1,1),IMDUMY(1,1)
      DATA SMALL / 1.0D-05 /
C
      IDBA=IDBUG(9)
      IDBB=IDBUG(10)
C
      IF (IDBA.LE.1) GO TO 10
      CALL EXAMIN (RMDUMY,1,1,INODE,1,NTOUCH,6,-17)
      CALL EXAMIN (RMDUMY,1,1,ISEG,1,NSEG,6,-16)
      CALL EXAMIN (RMDUMY,1,1,ITS,1,NTOUCH,6,-7)
      CALL EXAMIN (BETA,1,NTOUCH,IMDUMY,1,1,6,5)
   10 CONTINUE
C
      DO 20 K=1,NTOUCH
      ITS(K)=0
      BETA(K)=0.0D0
   20 CONTINUE
      DO 30 K=1,NSNOD
      DELTA(1,K)=0.0D0
      DELTA(2,K)=0.0D0
   30 CONTINUE
C
C
C
      ISTOP=0
      DO 400 I=1,NSURF
      IC=I
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,1)
      IF (ISR.EQ.0) GO TO 400
C
      DO 100 KC=KFC,KEC
C
      IRATIO=100
      DMIN=1.0D+20
      JSEG=0
      KSN=INODE(KC)
      KSNA=IABS(KSN)
      IF (KSN.EQ.-1) GO TO 100
C
C
C
      YC=XYZS(1,KC)
      ZC=XYZS(2,KC)
      NSA=KC+IC-1
      NSB=KC+IC
      XNC(1) = -T(2,NSA)-T(2,NSB)
      XNC(2) =  T(1,NSA)+T(1,NSB)
C
      DO 40 L=1,NSURFP
      IF (ISURFP(2,L).NE.IC) GO TO 40
      IT=ISURFP(1,L)
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,3)
C
      DO 110 K=KFT,KET
      YDIST=YC-XYZS(1,K)
      ZDIST=ZC-XYZS(2,K)
      DSQ = YDIST*YDIST + ZDIST*ZDIST
      IF (DSQ.GE.DMIN) GO TO 110
      DMIN=DSQ
      KNEAR=K
      ITNEAR=IT
  110 CONTINUE
   40 CONTINUE
C
      IT=ITNEAR
      CALL TRGET2 (ISURFP,IFSN,ISEG,NSURF,NSURFP,3)
      NSA=KNEAR+IT-1
      NSB=KNEAR+IT
      XNT(1) = -T(2,NSA)-T(2,NSB)
      XNT(2) =  T(1,NSA)+T(1,NSB)
C
C
C
      DIST(1) = YC - XYZS(1,KNEAR)
      DIST(2) = ZC - XYZS(2,KNEAR)
      DUMY = XNT(2)*DIST(1) - XNT(1)*DIST(2)
      XLSUM = XLN(NSA) + XLN(NSB)
      IF (ISEG(NSA).NE.0 .AND. ISEG(NSB).NE.0) GO TO 50
      IRATIO=INT( 100.0*DUMY/XLSUM )
      IF (IABS(IRATIO).GT.2) GO TO 50
C
C
C
      JSSA=IABS(ISEG(NSA))
      JSSB=IABS(ISEG(NSB))
      IF (JSSA.GT.0) JSEG=NSA
      IF (JSSB.GT.0) JSEG=NSB
      IF (JSEG.EQ.JFTD) JSEG=JLT
      IF (JSEG.EQ.JLTD) JSEG=JFT
      OVRLAP = XNT(1)*DIST(1) + XNT(2)*DIST(2)
      IF (KSNA.GT.1) GO TO 65
      IF (OVRLAP.LT.0.0D0) GO TO 100
C
   65 DUMY = XNC(1)*XNT(1) + XNC(2)*XNT(2)
      IF (DUMY.LT.0.0D0) GO TO 95
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,1000)
      write(66,1100) ISTOP,NODSF(KC),NODSF(KNEAR)
   95 IF (ISTOP.GT.0) GO TO 100
C
      IF (KSNA.EQ.1) INODE(KC)=2
      DELTA(1,KC) = OVRLAP*XNT(1)
      DELTA(2,KC) = OVRLAP*XNT(2)
      NA=JSEG-IT
      DIST(1) = (YC-XYZS(1,NA)) - DELTA(1,KC)
      DIST(2) = (ZC-XYZS(2,NA)) - DELTA(2,KC)
      DUMY = DIST(1)*T(1,JSEG) + DIST(2)*T(2,JSEG)
      GO TO 150
C
C
C
   50 JSEG=NSA
      IF (DUMY.GT.0.0D0) JSEG=NSB
   80 IF (ISEG(JSEG).EQ.0) GO TO 155
      IF (JSEG.EQ.JFTD) JSEG=JLT
      IF (JSEG.EQ.JLTD) JSEG=JFT
      NA=JSEG-IT
      NB=NA+1
      OVRLAP =-(YC-XYZS(1,NA))*T(2,JSEG)
     1        +(ZC-XYZS(2,NA))*T(1,JSEG)
C
      IF (INODE(KC).GT.1) GO TO 120
      IF (OVRLAP.LE.-1.0D-12) GO TO 100
C
  120 DUMY = XNC(1)*XNT(1) + XNC(2)*XNT(2)
      IF (DUMY.LT.0.0D0) GO TO 130
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,1000)
      write(66,1100) ISTOP,NODSF(KC),NODSF(KNEAR)
  130 IF (ISTOP.GT.0) GO TO 100
C
C
C
      IF (KSNA.EQ.1) INODE(KC)=2
      DELTA(1,KC)=-OVRLAP*T(2,JSEG)
      DELTA(2,KC)= OVRLAP*T(1,JSEG)
      DIST(1) = (YC-XYZS(1,NA)) - DELTA(1,KC)
      DIST(2) = (ZC-XYZS(2,NA)) - DELTA(2,KC)
      DUMY = DIST(1)*T(1,JSEG) + DIST(2)*T(2,JSEG)
C
C
C
      BDUMY1=-SMALL*XLN(JSEG)
      BDUMY2=(1.0+SMALL)*XLN(JSEG)
      IF (DUMY.LT.BDUMY1) GO TO 140
      IF (DUMY.LT.BDUMY2) GO TO 150
      DUMY=DUMY-XLN(JSEG)
      ISIDE=1
      GO TO 160
  140 ISIDE=-1
  160 LSEG=JSEG+ISIDE
      EXTRA = DUMY*( -T(1,JSEG)*T(2,LSEG) + T(2,JSEG)*T(1,LSEG) )
      DELTA(1,KC) = DELTA(1,KC) - EXTRA*T(2,LSEG)
      DELTA(2,KC) = DELTA(2,KC) + EXTRA*T(1,LSEG)
      JSEG=LSEG
      IF (ISEG(JSEG).EQ.0) GO TO 155
      IF (JSEG.EQ.JFTD) JSEG=JLT
      IF (JSEG.EQ.JLTD) JSEG=JFT
      NA=JSEG-IT
C
      DIST(1) = (YC-XYZS(1,NA)) - DELTA(1,KC)
      DIST(2) = (ZC-XYZS(2,NA)) - DELTA(2,KC)
      DUMY = DIST(1)*T(1,JSEG) + DIST(2)*T(2,JSEG)
C
C
C
  150 BT=DUMY/XLN(JSEG)
      BETA(KC)=BT
      ITS(KC)=JSEG
      GO TO 100
C
  155 INODE(KC)=1
C
  100 CONTINUE
C
C
      IF (ISTOP.EQ.0) GO TO 170
      GO TO 400
C
C
C
  170 IF (KEC.LT.KLC) INODE(KLC)=INODE(KFC)
      INODE(KFC)=IABS(INODE(KFC))
      DO 200 J=JFC,JLC
      NA=J-IC
      NB=NA+1
      KSNA=IABS(INODE(NA))
      KSNB=IABS(INODE(NB))
      INODE(NB)=KSNB
      IF (KSNA.EQ.1 .OR. KSNB.EQ.1) GO TO 210
      IF (ISEG(J).EQ.-1) GO TO 210
      IF (IABS(ISEG(J)).GT.1) GO TO 200
      ISEG(J)=2
      IF (KSNA.EQ.3 .OR. KSNB.EQ.3) ISEG(J)=3
      GO TO 200
  210 ISEG(J)=1
  200 CONTINUE