a37.for

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

      JOINT = 0
      IF (JFS.LE.KJOIN .AND. KJOIN.LE.JLS) JOINT = KJOIN
      RETURN
      END
      SUBROUTINE TPOLGY (ISURFP,IFSN,IFSE,NODSF,
     1                   NCA,
     1                   INODE,ISECT,
     1                   ITS,ITSP,JOIN,IDBUG,
     1                   XYZ,XYZS,
     1                   DELTA,VN,
     1                   NSURFP,NEALL,NTOUCH,NSNOD,NJOIN)
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           ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /SURF3/ IC,KFC,KLC,JFC,JLC
      COMMON /TRGT3/ IT,KFT,KLT,JFT,JLT
      COMMON /MATCH3/ ISR,IPAIR,JTSECT
      EQUIVALENCE ( NPAR(2),NSURF  ),( NPAR(7),NECON )
C
      DIMENSION ISURFP(2,NSURFP),IFSN(*),IFSE(*),NODSF(*),
     1          NCA(4,NEALL),
     1          INODE(*),ISECT(*),
     1          ITS(*),ITSP(*),JOIN(NJOIN,NSNOD),IDBUG(16)
      DIMENSION XYZ(3,NSNOD),XYZS(3,NSNOD),
     1          DELTA(3,NTOUCH),VN(3,NEALL)
      DIMENSION NCYCLE(8),VCYCLE(8)
      DIMENSION VC(3),VB(3),VT(3),XO(3)
      DIMENSION IMDUMY(1,1),RMDUMY(1,1)
      DATA SMALL / 1.0D-02 /
C
      ISTOP=0
      IDBA = IDBUG(9)
      IDBB = IDBUG(10)
C
      DO 40 K=1,NTOUCH
      ITS(K)=0
      ITSP(K)=0
      DO 40 L=1,3
   40 DELTA(L,K)=0.0D0
C
      IF (IDBB.GT.2) CALL LIGHT (XYZS,3,NSNOD,IMDUMY,1,1,5,3)
C
C
C
      DO 100 I=1,NSURF
C
      IC=I
      CALL TRGET3 (ISURFP,IFSN,IFSE,NSURF,NSURFP,1)
      IF (ISR.EQ.0) GO TO 100
C
      DO 120 KCON=KFC,KLC
C
      KSN  = INODE(KCON)
      KSNA = IABS(KSN)
      IF (KSN.EQ.-10) GO TO 120
      XKCON = XYZS(1,KCON)
      YKCON = XYZS(2,KCON)
      ZKCON = XYZS(3,KCON)
C
C
C
      DMIN = 1.0D+20
      DO 140 L=1,NSURFP
      IF (ISURFP(2,L).NE.IC) GO TO 140
      IT = ISURFP(1,L)
      CALL TRGET3 (ISURFP,IFSN,IFSE,NSURF,NSURFP,3)
C
      DO 160 K=KFT,KLT
      XDIST = XKCON - XYZS(1,K)
      YDIST = YKCON - XYZS(2,K)
      ZDIST = ZKCON - XYZS(3,K)
      DSQ = XDIST*XDIST + YDIST*YDIST + ZDIST*ZDIST
      IF (DSQ.GE.DMIN) GO TO 160
      DMIN = DSQ
      KNEAR=K
      ITNEAR=IT
  160 CONTINUE
C
  140 CONTINUE
      IF (IDBB.GT.2) write(66,2000) KCON,KNEAR
C
C
C
      DO 180 L=1,3
      VC(L)=0.0D0
      VB(L)=0.0D0
      DO 180 N=1,NJOIN
      JNC = JOIN(N,KCON)
      JNB = JOIN(N,KNEAR)
      CALL JOINTS (IFSN,IFSE,KCON ,JNC)
      CALL JOINTS (IFSN,IFSE,KNEAR,JNB)
      IF (JNC.GT.0) VC(L) = VC(L) + VN(L,JNC)
  180 IF (JNB.GT.0) VB(L) = VB(L) + VN(L,JNB)
C
      VCMAG = SQRT( VC(1)*VC(1) + VC(2)*VC(2) + VC(3)*VC(3) )
      VBMAG = SQRT( VB(1)*VB(1) + VB(2)*VB(2) + VB(3)*VB(3) )
C
      DO 200 L=1,3
      VC(L)=VC(L)/VCMAG
  200 VB(L)=VB(L)/VBMAG
C
C
C
      DO 220 N=1,NJOIN
      JN=JOIN(N,KNEAR)
      CALL JOINTS (IFSN,IFSE,KNEAR,JN)
      IF (JN.EQ.0) GO TO 220
C
      DO 240 L=1,4
      KNODE = NCA(L,JN)
      NCYCLE(L)   = KNODE
  240 NCYCLE(L+4) = KNODE
C
C
C
C
C
      DO 260 L=3,6
      LOCAL=L
  260 IF (NCYCLE(L).EQ.KNEAR) GO TO 280
C
  280 KB = KNEAR
      KA = NCYCLE(LOCAL-1)
      KC = NCYCLE(LOCAL+1)
      IF ( KA.EQ.KB ) KA = NCYCLE(LOCAL-2)
      IF ( KC.EQ.KB ) KC = NCYCLE(LOCAL+2)
C
      XB = XYZS(1,KB)
      YB = XYZS(2,KB)
      ZB = XYZS(3,KB)
C
      XA = XB - XYZS(1,KA)
      YA = YB - XYZS(2,KA)
      ZA = ZB - XYZS(3,KA)
C
      ASQ = XA*XA + YA*YA + ZA*ZA
      ASL = SQRT( ASQ )
      XA  = XA/ASL
      YA  = YA/ASL
      ZA =  ZA/ASL
C
      XC = XYZS(1,KC) - XB
      YC = XYZS(2,KC) - YB
      ZC = XYZS(3,KC) - ZB
C
      CSQ = XC*XC + YC*YC + ZC*ZC
      CSL = SQRT( CSQ )
      XC  = XC/CSL
      YC  = YC/CSL
      ZC  = ZC/CSL
C
      REFER = SMALL*SQRT( ASQ+CSQ )
C
      XKB = XKCON - XB
      YKB = YKCON - YB
      ZKB = ZKCON - ZB
C
C
      IF (IDBB.LE.2) GO TO 3000
      write(66,2240) KB,KA,KC
      write(66,2200) (VC(LC),LC=1,3),(VB(LB),LB=1,3)
      write(66,2210) XB,YB,ZB
      write(66,2210) XA,YA,ZA
      write(66,2210) XC,YC,ZC
      write(66,2210) XKB,YKB,ZKB
 3000 CONTINUE
C
      VOLA  =   XKB*(  YA*VB(3)  -  VB(2)*ZA  )
     1       +  YKB*(  ZA*VB(1)  -  VB(3)*XA  )
     1       +  ZKB*(  XA*VB(2)  -  VB(1)*YA  )
C
      VOLC  =   XKB*(  YC*VB(3)  -  VB(2)*ZC  )
     1       +  YKB*(  ZC*VB(1)  -  VB(3)*XC  )
     1       +  ZKB*(  XC*VB(2)  -  VB(1)*YC  )
C
      VOLA =   VOLA  +  REFER
      VOLC =   VOLC  +  REFER
      IF (VOLA.GE.0.0D0 .AND. VOLC.GE.0.0D0) GO TO 300
C
      IF (IDBB.LE.2) GO TO 220
      write(66,2010) KCON,KNEAR,JN,LOCAL,(NCYCLE(LN),LN=1,4),VOLA,VOLC
  220 CONTINUE
C
      INODE(KCON)=10
      IF (IDBB.GT.2) write(66,2020) KCON,KNEAR
      GO TO 120
C
C
C
  300 LOPP = LOCAL+2
      KOPP = NCYCLE(LOPP)
      IF (IDBB.GT.2) write(66,2030) KCON,KNEAR,LOPP,KOPP
C
      DO 320 L=1,3
      VT(L)=0.0D0
      DO 320 N=1,NJOIN
      JNT = JOIN(N,KOPP)
      CALL JOINTS (IFSN,IFSE,KOPP,JNT)
  320 IF (JNT.NE.0) VT(L) = VT(L) + VN(L,JNT)
C
      VTMAG = SQRT( VT(1)*VT(1) + VT(2)*VT(2) + VT(3)*VT(3) )
      VT(1) = VT(1)/VTMAG
      VT(2) = VT(2)/VTMAG
      VT(3) = VT(3)/VTMAG
C
C
C
C    (TRAILING NODE)
C
      DO 340 N=1,NJOIN
      JN = JOIN(N,KOPP)
      CALL JOINTS (IFSN,IFSE,KOPP,JN)
      IF (JN.EQ.0) GO TO 340
C
      DO 360 L=1,4
      KNODE = NCA(L,JN)
      NCYCLE(L)   = KNODE
  360 NCYCLE(L+4) = KNODE
C
      DO 380 L=3,6
      LOPP=L
  380 IF (NCYCLE(L).EQ.KOPP) GO TO 400
C
  400 KD = KOPP
      KC = NCYCLE(LOPP-1)
      KA = NCYCLE(LOPP+1)
      IF (KC.EQ.KD) KC = NCYCLE(LOPP-2)
      IF (KA.EQ.KD) KA = NCYCLE(LOPP+2)
C
      XD = XYZS(1,KD)
      YD = XYZS(2,KD)
      ZD = XYZS(3,KD)
C
      XC = XD - XYZS(1,KC)
      YC = YD - XYZS(2,KC)
      ZC = ZD - XYZS(3,KC)
C
      CSQ = XC*XC + YC*YC + ZC*ZC
      CSL = SQRT( CSQ )
      XC  = XC/CSL
      YC  = YC/CSL
      ZC  = ZC/CSL
C
      XA = XYZS(1,KA) - XD
      YA = XYZS(2,KA) - YD
      ZA = XYZS(3,KA) - ZD
C
      ASQ = XA*XA + YA*YA + ZA*ZA
      ASL = SQRT( ASQ )
      XA  = XA/ASL
      YA  = YA/ASL
      ZA  = ZA/ASL
C
      REFER = SMALL*SQRT( CSQ+ASQ )
C
      XKD = XKCON - XD
      YKD = YKCON - YD
      ZKD = ZKCON - ZD
C
C
      IF (IDBB.LE.2) GO TO 3010
      write(66,2240) KD,KC,KA
      write(66,2200) (VC(LC),LC=1,3),(VT(LT),LT=1,3)
      write(66,2210) XD,YD,ZD
      write(66,2210) XC,YC,ZC
      write(66,2210) XA,YA,ZA
      write(66,2210) XKD,YKD,ZKD
 3010 CONTINUE
      VOLC =   XKD*(  YC*VT(3)  -  VT(2)*ZC  )
     1      +  YKD*(  ZC*VT(1)  -  VT(3)*XC  )
     1      +  ZKD*(  XC*VT(2)  -  VT(1)*YC  )
C
      VOLA =   XKD*(  YA*VT(3)  -  VT(2)*ZA  )
     1      +  YKD*(  ZA*VT(1)  -  VT(3)*XA  )
     1      +  ZKD*(  XA*VT(2)  -  VT(1)*YA  )
C
      VOLC =   VOLC  +  REFER
      VOLA =   VOLA  +  REFER
      IF (VOLC.GE.0.0D0 .AND. VOLA.GE.0.0D0) GO TO 420
C
      IF (IDBB.LE.2) GO TO 340
      write(66,2040) KCON,KOPP,JN,LOPP,(NCYCLE(LN),LN=1,4),VOLC,VOLA
  340 CONTINUE
C
      INODE(KCON)=10
      IF (IDBB.GT.2) write(66,2050) KCON,KOPP
      GO TO 120
C
C
C
  420 JSAVE = JN
      VX = VN(1,JSAVE)
      VY = VN(2,JSAVE)
      VZ = VN(3,JSAVE)
C
      DO 440 L=1,3
      XSUM=0.0D0
      DO 460 M=1,4
      KNODE=NCYCLE(M)
  460 XSUM = XSUM + XYZS(L,KNODE)
  440 XO(L)=XSUM/4.0
C
      XKO = XKCON - XO(1)
      YKO = YKCON - XO(2)
      ZKO = ZKCON - XO(3)
C
      IF (IDBB.GT.2) write(66,2060) KCON,JSAVE,VX,VY,VZ,(XO(LX),LX=1,3)
C
C
C
C
C
      DO 480 L=1,4
      KNODE = NCYCLE(L)
      XR = XYZS(1,KNODE) - XO(1)
      YR = XYZS(2,KNODE) - XO(2)
      ZR = XYZS(3,KNODE) - XO(3)
C
      RSQ = XR*XR + YR*YR + ZR*ZR
      RSL = SQRT( RSQ )
      XR  = XR/RSL
      YR  = YR/RSL
      ZR  = ZR/RSL
C
      VOLUME =  XKO*(  YR*VZ  -  VY*ZR  )
     1        + YKO*(  ZR*VX  -  VZ*XR  )
     1        + ZKO*(  XR*VY  -  VX*YR  )
C
      VCYCLE(L  ) = VOLUME
      VCYCLE(L+4) = VOLUME
  480 CONTINUE
      IF (IDBB.GT.2) write(66,2070) KCON,JSAVE,(VCYCLE(NV),NV=1,4)
C
      DO 500 L=1,4
      LA = L
      LB = L+1
      KA = NCYCLE(LA)
      KB = NCYCLE(LB)
      IF (KA.EQ.KB) GO TO 500
      VOLB = -VCYCLE(LB) + SMALL
      VOLA =  VCYCLE(LA) + SMALL
      IF (VOLA.GE.0.0D0 .AND. VOLB.GE.0.0D0) GO TO 520
  500 CONTINUE
C
  520 LSAVE = LA
      XB = XO(1) - XYZS(1,KB)
      YB = XO(2) - XYZS(2,KB)
      ZB = XO(3) - XYZS(3,KB)
C
      XA = XYZS(1,KA) - XO(1)
      YA = XYZS(2,KA) - XO(2)
      ZA = XYZS(3,KA) - XO(3)
C
      VOX = YA*ZB - YB*ZA
      VOY = ZA*XB - ZB*XA
      VOZ = XA*YB - XB*YA
      VOMAG = SQRT( VOX*VOX + VOY*VOY + VOZ*VOZ )
C
      VOX = VOX/VOMAG
      VOY = VOY/VOMAG
      VOZ = VOZ/VOMAG
C
      COSINE =    VX*VOX +  VY*VOY +  VZ*VOZ
      OVRLAP = ( XKO*VOX + YKO*VOY + ZKO*VOZ )/COSINE
C
      IF (KSN.GT.10) GO TO 540
      IF (OVRLAP.GT.-1.0D-12) GO TO 540
      INODE(KCON)=10
      GO TO 120
C
  540 DOT = VX*VC(1) + VY*VC(2) + VZ*VC(3)
      IF (DOT.LE.0.0D0) GO TO 580
      ISTOP=ISTOP+1
      IF (ISTOP.EQ.1) write(66,1000)
      write(66,1100) ISTOP,NODSF(KCON),NODSF(KNEAR)
  580 IF (ISTOP.GT.0) GO TO 120
C
      ITS(KCON) = JSAVE
      ITSP(KCON)= LSAVE
      IF (KSNA.EQ.10) INODE(KCON) = 20
C
      DELTA(1,KCON) = OVRLAP*VX
      DELTA(2,KCON) = OVRLAP*VY
      DELTA(3,KCON) = OVRLAP*VZ
C
  120 CONTINUE
  100 CONTINUE
C
      IF (IDBA.EQ.0) GO TO 390
      CALL LIGHT (RMDUMY,1,1,INODE,1,NTOUCH,5,-14)
      CALL LIGHT (RMDUMY,1,1,ISECT,1,NECON,5,-16)
  390 IF (ISTOP.GT.0) GO TO 720
C
C
C
      DO 640 K=1,NTOUCH
      KSN = INODE(K)
  640 INODE(K) = IABS( KSN )
C
      DO 660 J=1,NECON
      I10 = 0
      I30 = 0
      DO 680 L=1,4
      KNODE=NCA(L,J)
      KSN  = INODE(KNODE)
      IF (KSN.EQ.10) I10 = 1
  680 IF (KSN.EQ.30) I30 = 1
      IF (I10.EQ.1) GO TO 700
C
      IF (ISECT(J).EQ.-10) GO TO 700
      IF (ISECT(J).GT. 10) GO TO 660
      ISECT(J)=20
      IF (I30.EQ.1) ISECT(J)=30
      GO TO 660
C
  700 ISECT(J)=10
  660 CONTINUE
C
      IF (IDBA.EQ.0) GO TO 720
      CALL LIGHT (RMDUMY,1,1,INODE,1,NTOUCH,5,-15)
      CALL LIGHT (RMDUMY,1,1,ISECT,1,NECON,5,-17)
      CALL LIGHT (RMDUMY,1,1,ITS,1,NTOUCH,5,-18)
      CALL LIGHT (RMDUMY,1,1,ITSP,1,NTOUCH,5,-19)
      CALL LIGHT (DELTA,3,NTOUCH,IMDUMY,1,1,5,8)
C
  720 IF (ISTOP.EQ.0) RETURN
      write(66,1150)
      DO 730 K=1,NSNOD
      write(66,1160) NODSF(K),(XYZS(L,K),L=1,3),(XYZ(L,K),L=1,3)
  730 CONTINUE
      write(66,1200)
      STOP
C
 1000 FORMAT(/,10X,34HCONDITIONS OF CONTACT INADMISSIBLE,/,
     1         11X,33HDUE TO DETECTION OF GROSS OVERLAP,/,
     1         14X,26HAT THE FOLLOWING LOCATIONS,/)
 1100 FORMAT(2X,I5,22H). CONTACTOR NODE NO.= ,I5,5X,
     1                16HTARGET NODE NO.= ,I5)
 1150 FORMAT(/,30X,27HCURRENT SURFACE COORDINATES,28X,
     1             27HINITIAL SURFACE COORDINATES,/,
     1          6X,6HGLOBAL,8X,12HX-COORDINATE,5X,12HY-COORDINATE,
     1          5X,12HZ-COORDINATE,9X,12HX-COORDINATE,
     1          5X,12HY-COORDINATE,5X,12HZ-COORDINATE,/,
     1          5X,8HNODE NO.,/)
 1160 FORMAT(I10,5X,3(5X,E12.4),4X,3(5X,E12.4))
 1200 FORMAT(//,10X,34HSTOPPED IN THE 3-D CONTACT OVERLAY   )
C
C
 2000 FORMAT(/,1X,7HPHASE 1,2X,5HKCON=,I3,2X,6HKNEAR=,I3 )
 2010 FORMAT(  1X,7HPHASE 1,2X,5HKCON=,I3,2X,6HKNEAR=,I3,2X,3HJN=,I3,
     1         2X,6HLOCAL=,I3,2X,7HNCYCLE=,4I3,2X,5HVOLA=,E12.4,
     1         2X,5HVOLC=,E12.4 )
 2020 FORMAT(/,1X,21HTERMINATED AT PHASE 1,2X,5HKCON=,I3,2X,
     1             6HKNEAR=,I3 )
 2030 FORMAT(/,1X,7HPHASE 2,2X,5HKCON=,I3,2X,6HKNEAR=,I3,2X,
     1            5HLOPP=,I3,2X,5HKOPP=,I3 )
 2040 FORMAT(  1X,7HPHASE 2,2X,5HKCON=,I3,3X,5HKOPP=,I3,2X,3HJN=,I3,
     1         3X,5HLOPP=,I3,2X,7HNCYCLE=,4I3,2X,5HVOLC=,E12.4,
     1         2X,5HVOLA=,E12.4 )
 2050 FORMAT(/,1X,21HTERMINATED AT PHASE 2,2X,5HKCON=,I3,3X,
     1            5HKOPP=,I3 )
 2060 FORMAT(/,1X,7HPHASE 3,2X,5HKCON=,I3,2X,6HJSAVE=,I3,2X,
     1             7HNORMAL=,3E12.4,3X,7HCENTER=,3E12.4 )
 2070 FORMAT(  1X,7HPHASE 3,2X,5HKCON=,I3,2X,6HJSAVE=,I3,2X,