a37.for

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

     1          NODSF(*),NCA(4,NEALL),
     1          INODE(*),ISECT(*),
     1          ITS(*),ITSP(*),JOIN(NJOIN,NSNOD),IDBUG(16)
      DIMENSION FCOFF(3,NSURFP),XYZS(3,NSNOD),
     1          CPR(3,NTOUCH),DELTA(3,NTOUCH),GUSTAV(4,NECON),
     1          VN(3,NEALL),CFOR(3,NSNOD)
      DIMENSION PS(4,3),PN(4,3),PT(4,3),NCYCLE(8)
      DIMENSION IMDUMY(1,1),RMDUMY(1,1)
      EQUIVALENCE ( NPAR(2),NSURF ),( NPAR(13),ICPRNT )
C
      DATA LONELY / 100 / , SMALL / 1.0D-12 / , PRCENT / 0.01D0 /
C
C
C
C
C
      IDBA = IDBUG(7)
      IDBB = IDBUG(8)
      KCPRI=1
      IF (IDBA.LE.1) GO TO 10
      CALL LIGHT (RMDUMY,1,1,INODE,1,NTOUCH,4,-14)
      CALL LIGHT (RMDUMY,1,1,ISECT,1,NECON ,4,-16)
      CALL LIGHT (CPR,3,NTOUCH,IMDUMY,1,1,4,6)
   10 IF (KPRI.EQ.0 .OR. IDBB.GE.1) KCPRI=0
      IF (KPRI.EQ.0) CALL PRINT3 (IPS,ISV,0)
C
      DO 12 K=1,NTOUCH
      DO 12 L=1,3
   12 CFOR(L,K)=0.0D0
C
      DO 100 I=1,NSURF
      IS=I
      CALL TRGET3 (ISURFP,IFSN,IFSE,NSURF,NSURFP,1)
      IF (ISR.EQ.0) GO TO 100
      IF (KCPRI.EQ.0) CALL PRINT3 (IPS,ISV,1)
C
      DO 120 J=JFS,JLS
      JSE  = J
      JSS  = ISECT(J)
      JSSA = IABS(JSS)
C
      KA = NCA(1,J)
      KB = NCA(2,J)
      KC = NCA(3,J)
      KD = NCA(4,J)
C
      NSE(1) = NODSF(KA)
      NSE(2) = NODSF(KB)
      NSE(3) = NODSF(KC)
      NSE(4) = NODSF(KD)
C
      PNAV = 0.0D0
      PTAV = 0.0D0
      RNALL = 0.0D0
      RTALL = 0.0D0
      DO 130 L=1,3
      RN(L)=0.0D0
      RT(L)=0.0D0
      RS(L)=0.0D0
      VNJ(L) = VN(L,J)
      DO 132 N=1,4
      PN(N,L)=0.0D0
      PT(N,L)=0.0D0
  132 PS(N,L)=0.0D0
  130 CONTINUE
      IF (JSSA.EQ.10) GO TO 145
C
C
C
      XJ1=GUSTAV(1,J)
      XJ2=GUSTAV(2,J)
      XJ3=GUSTAV(3,J)
      XJ4=GUSTAV(4,J)
      AREA = XJ1 + XJ2 + XJ3 + XJ4
C
      S1  = C1*XJ1 + C3*( XJ2 + XJ4 ) + C5*XJ3
      S5  = C1*XJ2 + C3*( XJ3 + XJ1 ) + C5*XJ4
      S8  = C1*XJ3 + C3*( XJ4 + XJ2 ) + C5*XJ1
      S10 = C1*XJ4 + C3*( XJ1 + XJ3 ) + C5*XJ2
C
      S2  = C2*( XJ1 + XJ2 ) + C4*( XJ3 + XJ4 )
      S6  = C2*( XJ2 + XJ3 ) + C4*( XJ4 + XJ1 )
      S9  = C2*( XJ3 + XJ4 ) + C4*( XJ1 + XJ2 )
      S4  = C2*( XJ4 + XJ1 ) + C4*( XJ2 + XJ3 )
C
      S3  = C3*AREA
      S7  = S3
C
      SA = S1 + S2 + S3 + S4
      SB = S2 + S5 + S6 + S7
      SC = S3 + S6 + S8 + S9
      SD = S4 + S7 + S9 + S10
C
C
C
      DO 150 L=1,3
      PS(1,L) = CPR(L,KA)
      PS(2,L) = CPR(L,KB)
      PS(3,L) = CPR(L,KC)
      PS(4,L) = CPR(L,KD)
      RS(L)   = SA*PS(1,L) + SB*PS(2,L) + SC*PS(3,L) + SD*PS(4,L)
  150 CONTINUE
C
      RNALL = RS(1)*VNJ(1) + RS(2)*VNJ(2) + RS(3)*VNJ(3)
C
      DO 152 L=1,3
      VTEM = VNJ(L)
      RN(L) = RNALL*VTEM
      RT(L) = RS(L) - RN(L)
  152 CONTINUE
C
      DO 154 N=1,4
      PNN = PS(N,1)*VNJ(1) + PS(N,2)*VNJ(2) + PS(N,3)*VNJ(3)
      DO 155 L=1,3
      VTEM = VNJ(L)
      PN(N,L) = PNN*VTEM
  155 PT(N,L) = PS(N,L) - PN(N,L)
  154 CONTINUE
C
      RTSQ = RT(1)*RT(1) + RT(2)*RT(2) + RT(3)*RT(3)
      IF (RTSQ.GT.0.0D0) RTALL = SQRT(RTSQ)
C
      PNAV = RNALL/AREA
      PTAV = RTALL/AREA
C
C
C
      CHECK = PNAV
      IF (CHECK.GE.0.0D0) GO TO 180
      PNAV=0.0D0
      PTAV=0.0D0
      RNALL=0.0D0
      RTALL=0.0D0
      DO 165 L=1,3
      RN(L)=0.0D0
      RT(L)=0.0D0
      RS(L)=0.0D0
      DO 168 N=1,4
      PN(N,L)=0.0D0
      PT(N,L)=0.0D0
  168 PS(N,L)=0.0D0
  165 CONTINUE
      ISECT(J)=-10
      GO TO 145
C
C
C
  180 JTSECT=ITS(KA)
      CALL TRGET3 (ISURFP,IFSN,IFSE,NSURF,NSURFP,2)
      PCF = FCOFF(1,IPAIR)
C
      PFAV = PCF*PNAV
      REDUCE = 1.0D0
      IF (JSSA.EQ.30) REDUCE = 1.0 - PRCENT
      CHECK = REDUCE*PFAV - PTAV
      IF (CHECK.GE.0.0D0) GO TO 200
C
C
      ISECT(J)=30
      PFAV = PCF*PNAV
      PTMIN = MIN(PTAV,PFAV)
      PTMAX = MAX(PTAV,SMALL)
      RATIO = PTMIN/PTMAX
      PTAV  = RATIO*PTAV
      DO 172 L=1,3
      RTEM = RATIO*RT(L)
      PTEM = RTEM/AREA
      RT(L) = RTEM
      RS(L) = RTEM + RN(L)
      DO 174 N=1,4
      PT(N,L) = PTEM
  174 PS(N,L) = PTEM + PN(N,L)
  172 CONTINUE
      GO TO 300
C
  200 ISECT(J)=20
C
C
C
  300 DO 302 L=1,3
      PA = PS(1,L)
      PB = PS(2,L)
      PC = PS(3,L)
      PD = PS(4,L)
      CFOR(L,KA) = CFOR(L,KA) + PA*S1 + PB*S2 + PC*S3 + PD*S4
      CFOR(L,KB) = CFOR(L,KB) + PA*S2 + PB*S5 + PC*S6 + PD*S7
      CFOR(L,KC) = CFOR(L,KC) + PA*S3 + PB*S6 + PC*S8 + PD*S9
      CFOR(L,KD) = CFOR(L,KD) + PA*S4 + PB*S7 + PC*S9 + PD*S10
  302 CONTINUE
  145 JSS = ISECT(J)
      IF (KCPRI.EQ.0) CALL PRINT3 (IPS,ISV,2)
C
  120 CONTINUE
C
C
C
C
      NOSOL=0
      ITITLE=1
      IF (KCPRI.EQ.0) ITITLE=0
      NTOTAL = KLS - KFS + 1
C
      DO 330 K=KFS,KLS
      KA=K
      RNALL = 0.0D0
      RTALL = 0.0D0
      KSN   = INODE(K)
      KSNA  = IABS( KSN )
      NSE(1)= NODSF(KA)
C
      DO 335 L=1,3
      RN(L) =0.0D0
      RT(L) =0.0D0
      RS(L) =0.0D0
  335 VNJ(L)=0.0D0
C
      IF (KSNA.EQ.10) GO TO 332
      DO 350 N=1,NJOIN
      JN=JOIN(N,K)
      CALL JOINTS (IFSN,IFSE,KA,JN)
      IF (JN.EQ.0) GO TO 350
      JSS = ISECT(JN)
      IF (JSS.NE.10) GO TO 332
      DO 360 L=1,3
  360 VNJ(L) = VNJ(L) + VN(L,JN)
  350 CONTINUE
C
      IF (ITITLE.EQ.0) CALL PRINT3 (IPS,ISV,3)
      ITITLE=1
C
      VNJM = SQRT( VNJ(1)*VNJ(1) + VNJ(2)*VNJ(2) + VNJ(3)*VNJ(3) )
C
      DO 380 L=1,3
      RTEM = CPR(L,K)
      RS(L) = RTEM
      VNJ(L)=VNJ(L)/VNJM
  380 RNALL = RNALL + RTEM*VNJ(L)
C
      DO 385 L=1,3
      RN(L) = RNALL*VNJ(L)
  385 RT(L) = RS(L) - RN(L)
      RTALL = SQRT( RT(1)*RT(1) + RT(2)*RT(2) + RT(3)*RT(3) )
C
      CHECK = RNALL
      IF (CHECK.GE.0.0D0) GO TO 390
      INODE(K) = -10*LONELY
      DO 400 L=1,3
      RN(L)=0.0D0
      RT(L)=0.0D0
      RS(L)=0.0D0
  400 CFOR(L,K)=0.0D0
      GO TO 332
C
  390 JTSECT = ITS(K)
      CALL TRGET3 (ISURFP,ISFN,IFSE,NSURF,NSURFP,2)
      PCF = FCOFF(1,IPAIR)
      RESIST = PCF*RNALL
C
      REDUCE = 1.0D0
      IF (KSNA.EQ.30) REDUCE = 1.0 - PRCENT
      CHECK = REDUCE*RESIST - RTALL
      IF (CHECK.GE.0.0D0) GO TO 430
C
      INODE(K) = 30*LONELY
      RESIST= PCF*RNALL
      RTMIN = MIN(RTALL,RESIST)
      RTMAX = MAX(RTALL,SMALL)
      RATIO = RTMIN/RTMAX
      DO 420 L=1,3
      RS(L) = RN(L) + RATIO*RT(L)
  420 CFOR(L,K) = RS(L)
      GO TO 340
C
  430 INODE(K) = 20*LONELY
      DO 440 L=1,3
  440 CFOR(L,K) = RS(L)
C
  340 IF (KCPRI.EQ.0) CALL PRINT3 (IPS,ISV,5)
      GO TO 330
C
  332 NOSOL = NOSOL + 1
C
  330 CONTINUE
C
      IF (NOSOL.LT.NTOTAL) GO TO 345
      IF (KCPRI.EQ.0) CALL PRINT3 (IPS,ISV,4)
C
  345 IF (KPRI.EQ.0) GO TO 100
C
C
C
      DO 435 K=KFS,KLS
      KA=K
      KSN=INODE(K)
      KSNA=IABS(KSN)
      IF (KSNA.LT.LONELY) GO TO 445
      INODE(K)=KSN/LONELY
      GO TO 435
C
  445 INODE(K)=10
      DO 450 NCASE=1,3
      DO 460 N=1,NJOIN
      JN = JOIN(N,K)
      CALL JOINTS (IFSN,IFSE,KA,JN)
      IF (JN.EQ.0) GO TO 460
      JSS=ISECT(JN)
      JSSA=IABS(JSS)
C
      IF (NCASE-2) 470,480,490
C
  470 IF (JSSA.EQ.20) INODE(K)=20
      GO TO 500
  480 IF (JSSA.EQ.30) INODE(K)=30
      GO TO 500
  490 IF (JSS.EQ.-10) INODE(K)=-10
C
  500 IF (INODE(K).NE.10) GO TO 435
  460 CONTINUE
  450 CONTINUE
C
  435 CONTINUE
  100 CONTINUE
C
C
C
      DO 510 K=1,NTOUCH
      DO 510 L=1,3
  510 CPR(L,K)=CFOR(L,K)
C
      IF (IDBA.EQ.0) GO TO 520
      CALL LIGHT (CPR,3,NTOUCH,IMDUMY,1,1,4,7)
      CALL LIGHT (RMDUMY,1,1,INODE,1,NTOUCH,4,-15)
      CALL LIGHT (RMDUMY,1,1,ISECT,1,NECON,4,-17)
C
C
C
  520 IF (KPRI.GT.0 .OR. ICPRNT.EQ.2) RETURN
C
      DO 525 K=1,NSNOD
      DO 525 L=1,3
  525 CFOR(L,K)=0.0D0
C
      DO 540 I=1,NSURF
      IS=I
      CALL TRGET3(ISURFP,IFSN,IFSE,NSURF,NSURFP,1)
      IF (ISR.EQ.0) GO TO 540
      DO 550 K=KFS,KLS
      KCON = K
      KSNA = IABS( INODE(K) )
      IF (KSNA.EQ.10) GO TO 550
C
      LA=ITSP(K)
      JTSECT=ITS(K)
      DO 555 L=1,4
      LNODE=NCA(L,JTSECT)
      NCYCLE(L)   = LNODE
  555 NCYCLE(L+4) = LNODE
      KLA = NCYCLE(LA)
      KLB = NCYCLE(LA+1)
      KLC = NCYCLE(LA+2)
      KLD = NCYCLE(LA+3)
C
      CALL CPOINT (XYZS,DELTA,ITSP,NCYCLE,NTOUCH,NSNOD)
C
      GA = (1.0 - ALFA - BETA)/4.0
      HA = ALFA + GA
      HB = BETA + GA
      HC = GA
      HD = GA
C
      DO 560 L=1,3
      CPRCON = CPR(L,K)
      CFOR(L,K)   = CFOR(L,K)   +    CPRCON
      CFOR(L,KLA) = CFOR(L,KLA) - HA*CPRCON
      CFOR(L,KLB) = CFOR(L,KLB) - HB*CPRCON
      CFOR(L,KLC) = CFOR(L,KLC) - HC*CPRCON
      CFOR(L,KLD) = CFOR(L,KLD) - HD*CPRCON
  560 CONTINUE
  550 CONTINUE
  540 CONTINUE
C
      DO 570 I=1,NSURF
      IS=I
      CALL TRGET3 (ISURFP,IFSN,IFSE,NSURF,NSURFP,1)
      CALL PRINT3 (IPS,ISV,6)
      DO 580 K=KFS,KLS
      KA=K
      NSE(1)= NODSF(KA)
      RS(1) = CFOR(1,KA)
      RS(2) = CFOR(2,KA)
      RS(3) = CFOR(3,KA)
      CALL PRINT3 (IPS,ISV,7)
  580 CONTINUE
  570 CONTINUE
C
      RETURN
      END
      SUBROUTINE TRGET3 (ISURFP,IFSN,IFSE,NSURF,NSURFP,IFLAG)
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /SURF3/ IC,KFC,KLC,JFC,JLC
      COMMON /TRGT3/ IT,KFT,KLT,JFT,JLT
      COMMON /MATCH3/ ISR,IPAIR,JTSECT
C
      DIMENSION ISURFP(2,NSURFP),IFSN(*),IFSE(*)
C
      GO TO ( 50,100,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
  100 DO 120 I=1,NSURF
      IF (I.EQ.IC) GO TO 120
      ITRGET=I
      JFSECT=IFSE(I)
      JLSECT=IFSE(I+1)-1
      IF (JFSECT.LE.JTSECT .AND. JTSECT.LE.JLSECT) GO TO 110
  120 CONTINUE
C
C
C
  110 DO 150 L=1,NSURFP
      IPAIR=L
      IF (ISURFP(1,L).EQ.ITRGET .AND. ISURFP(1,L).EQ.IC) GO TO 160
  150 CONTINUE
  160 RETURN
C
  200 KFT=IFSN(IT)
      KLT=IFSN(IT+1)-1
      JFT=IFSE(IT)
      JLT=IFSE(IT+1)-1
  250 KFC=IFSN(IC)
      KLC=IFSN(IC+1)-1
      JFC=IFSE(IC)
      JLC=IFSE(IC+1)-1
      RETURN
C
      END
      SUBROUTINE JOINTS (IFSN,IFSE,NODE,JOINT)
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
      EQUIVALENCE (NPAR(2),NSURF)
      DIMENSION IFSN(*),IFSE(*)
C
      IF (JOINT.EQ.0) RETURN
      DO 100 I=1,NSURF
      ISURF=I
      KFS=IFSN(I)
      KLS=IFSN(I+1)-1
      IF (KFS.LE.NODE .AND. NODE.LE.KLS) GO TO 120
  100 CONTINUE
C
  120 JFS = IFSE(ISURF)
      JLS = IFSE(ISURF+1)-1
      KJOIN = JOINT