Fortran Coder

标题: 老代码翻译 [打印本页]

作者: little_kar 时间: 2017-3-8 16:31
标题: 老代码翻译
各位前辈，在其他地方找来一段代码，是有限元解二维热传导方程的，代码比较老，看着很费劲，想劳烦谁能适当翻译一下，不用全部，比如C -----------------  HEAT2D  -----------------------
   DIMENSION X(200,2),NOC(300,3),MAT(300),PM(10,3),
   $ EHS(300),NU(100),U(100),F(400),BT(2,3),S(400,95)
   CHARACTER*16 FILE1,FILE2,FILE3
   CHARACTER*81 DUMMY,TITLE
   PRINT *, '***************************************'
   PRINT *, '*       PROGRAM  HEAT2D          *'
   PRINT *, '* HEAT 2-D  WITH 3-NODED TRIANGLES  *'
   PRINT *, '* T.R.Chandrupatla and A.D.Belegundu  *'
   PRINT *, '***************************************'
C    IMAX = FIRST DIMENSION OF THE S-MATRIX, IX=1ST DIM. OF X-MATRIX, Etc.
   IMAX = 400
   IX=200
   INOC=300
   IPM=10
   PRINT *,'Input Data File Name <DOS file name>'
   READ '(A)',FILE1
   LINP=10
   OPEN(UNIT=10,FILE=FILE1,STATUS='UNKNOWN')
   PRINT *,'Output Data File Name <DOS file name>'
   READ '(A)',FILE2
   LOUT=11
   OPEN(UNIT=11,FILE=FILE2,STATUS='UNKNOWN')
   READ(LINP,'(A)')DUMMY
   READ(LINP,'(A)')TITLE
   READ(LINP,'(A)')DUMMY
   READ(LINP,*) NN, NE, NM, NDIM, NEN, NDN
   READ(LINP,'(A)')DUMMY
   READ(LINP,*) ND, NL, NMPC
C --- Material property Thermal Conductivity
   NPR = 1
   NMPC = 0
   NDN = 1
   NDIM=2
   NEN=3
C    ELEMENT HEAT SOURCE, EHS(I),I=1,...,NE
C    --- ND = NO. OF SPECIFIED TEMPERATURES
C    --- NL = NO. OF NODAL HEAT SOURCES
C    ---  NPR =1 (THERMAL CONDUCTIVITY) AND NMPC = 0
   PRINT *,'PLOT CHOICE'
   PRINT *,'  1) No Plot Data'
   PRINT *,'  2) Create Data File for Temperatures'
   PRINT *,'Choose 1 or 2 '
   READ(5,*) IPL
   IF (IPL .LT. 1 .OR. IPL .GT. 2)IPL = 1
C    --- default is no data
   IF (IPL .GT. 1) THEN
      PRINT *,'Give File Name for Plot Data'
READ '(A)',FILE3
      LOUT2=12
OPEN(UNIT=12,FILE=FILE3,STATUS='UNKNOWN')
   ENDIF
C    ----- Coordinates
   READ(LINP,'(A)')DUMMY
   DO 39 I = 1, NN
      READ(LINP,*) N ,(X(N,J),J=1,NDIM)
39 CONTINUE
C    ----- Connectivity, Material, Element Heat Source
   READ(LINP,'(A)')DUMMY
   DO 41 I = 1,NE
      READ(LINP,*) N, (NOC(N,J),J=1,NEN),MAT(N),EHS(N)
41 CONTINUE
C    ----- Temperature BC
   READ(LINP,'(A)')DUMMY
   IF (ND.GT.0)THEN
      DO 50 I=1,ND
      READ(LINP,*)NU(I), U(I)
50    CONTINUE
   ENDIF
C    ----- Nodal Heat Sources
   DO 42 I=1,NN
      F(I)=0.
42 CONTINUE
   READ(LINP,'(A)')DUMMY
   IF (NL.GT.0)THEN
      DO 51 I=1,NL
      READ(LINP,*) N, F(N)
51    CONTINUE
   ENDIF
C    ----- Thermal Conductivity
   READ(LINP,'(A)')DUMMY
   READ (LINP,*)(N,(PM(N,J),J=1,NPR),I=1,NM)
C    --- BAND WIDTH CALCULATION ---
   NBW = 0
   DO 5 I = 1,NE
      NMIN = NOC(I, 1)
      NMAX = NOC(I, 1)
      DO 3 J = 2,3
         IF (NMIN .GT. NOC(I, J)) NMIN = NOC(I, J)
         IF (NMAX .LT. NOC(I, J)) NMAX = NOC(I, J)
3       CONTINUE
      NTMP = NDN * (NMAX - NMIN + 1)
      IF (NBW .LT. NTMP) NBW = NTMP
5    CONTINUE
   WRITE(6,*)'THE BAND WIDTH IS',NBW
C    --- INITIALIZATION OF CONDUCTIVITY MATRIX
   DO 8 I = 1, NN
   DO 8 J = 1,NBW
      S(I, J) = 0.
8    CONTINUE
   READ(LINP,'(A)')DUMMY
   READ (LINP,*) NHF
   IF (NHF .GT. 0) THEN
      DO 11 I=1,NHF
         READ(LINP,*) N1, N2, V
         ELEN = SQRT((X(N1, 1) - X(N2, 1)) ** 2 + (X(N1, 2) -
   $          X(N2, 2)) ** 2)
         F(N1) = F(N1) - ELEN * V / 2
         F(N2) = F(N2) - ELEN * V / 2
11 CONTINUE
   END IF
   READ(LINP,'(A)')DUMMY
   READ(LINP,*)NCONV
   IF (NCONV .GT. 0) THEN
      DO 13 I = 1,NCONV
         READ(LINP,*)N1, N2, H, TINF
         ELEN = SQRT((X(N1, 1) - X(N2, 1))**2 +
   $          (X(N1, 2) - X(N2, 2)) ** 2)
         F(N1) = F(N1) + ELEN * H * TINF / 2
         F(N2) = F(N2) + ELEN * H * TINF / 2
         S(N1, 1) = S(N1, 1) + H * ELEN / 3
         S(N2, 1) = S(N2, 1) + H * ELEN / 3
         IF (N1 .GE. N2) THEN
            N3 = N1
            N1 = N2
            N2 = N3
         END IF
         S(N1, N2 - N1 + 1) = S(N1, N2 - N1 + 1) + H * ELEN / 6
13 CONTINUE
   END IF
   CLOSE(LINP)
C    --- CONDUCTIVITY MATRIX
   DO 15 I = 1, NE
      I1 = NOC(I, 1)
      I2 = NOC(I, 2)
      I3 = NOC(I, 3)
      X32 = X(I3, 1) - X(I2, 1)
      X13 = X(I1, 1) - X(I3, 1)
      X21 = X(I2, 1) - X(I1, 1)
      Y23 = X(I2, 2) - X(I3, 2)
      Y31 = X(I3, 2) - X(I1, 2)
      Y12 = X(I1, 2) - X(I2, 2)
      DETJ = X13 * Y23 - X32 * Y31
      AREA = .5 * ABS(DETJ)
C       --- ELEMENT HEAT SOURCES
      IF (EHS(I) .NE.0.) THEN
         C = EHS(I) * AREA / 3
         F(I1) = F(I1) + C
         F(I2) = F(I2) + C
         F(I3) = F(I3) + C
      END IF
      BT(1, 1) = Y23
      BT(1, 2) = Y31
      BT(1, 3) = Y12
      BT(2, 1) = X32
      BT(2, 2) = X13
      BT(2, 3) = X21
      DO 19 II = 1, 3
         DO 17 JJ = 1,2
            BT(JJ, II) = BT(JJ, II) / DETJ
17       CONTINUE
19    CONTINUE
      DO 21 II = 1,3
         DO 23 JJ = 1,3
            II1 = NOC(I, II)
            II2 = NOC(I, JJ)
            IF (II1 .LE.II2) THEN
               SUM = 0.
               DO 25 J = 1,2
                  SUM = SUM + BT(J, II) * BT(J, JJ)
25             CONTINUE
               IC = II2 - II1 + 1
               S(II1, IC) = S(II1, IC) + SUM * AREA * PM(MAT(I), 1)
            END IF
23       CONTINUE
21    CONTINUE
15 CONTINUE
   IF (ND .GT. 0) THEN
C'---- - Decide Penalty Parameter CNST -----
      CNST = 0.
      DO 27 I = 1, NN
      IF (CNST .LT. S(I, 1)) CNST = S(I, 1)
27    CONTINUE
      CNST = CNST * 10000.
C'---- - Modify for Boundary Conditions -----
C       --- Temperature BC ---
      DO 29 I = 1, ND
      N = NU(I)
      S(N, 1) = S(N, 1) + CNST
      F(N) = F(N) + CNST * U(I)
29    CONTINUE
   END IF
C    --- EQUATION SOLVING
   CALL BANSOL(NN,NBW,IMAX,S,F)
   WRITE (LOUT,'(1X,2A)') 'Output for Input Data File: ',FILE1
   WRITE(LOUT,'(1X,A)') TITLE
   PRINT *, TITLE
   WRITE(LOUT,*)'NODE# Temperature'
   PRINT *,'NODE# Temperature'
   DO 33 I = 1,NN
   WRITE(*,'(I5,E15.5)') I, F(I)
   WRITE(LOUT,'(I5,E15.5)') I, F(I)
33 CONTINUE
   IF (IPL .GT. 1) THEN
WRITE(LOUT2,*)'(Node Temp.) for Data in File ',FILE1
      DO 60 I=1,NN
   WRITE(LOUT2,*) I,F(I)
60    CONTINUE
   END IF
   CLOSE (LOUT2)
   WRITE(LOUT,*)' -- CONDUCTION HEAT FLOW PER UNIT AREA IN EACH ELEME
   $NT -- '
   WRITE(LOUT,*)'ELEMENT# QX= -K*DT/DX QY= -K*DT/DY '
   DO 45 I = 1,NE
      I1 = NOC(I, 1)
      I2 = NOC(I, 2)
      I3 = NOC(I, 3)
      X32 = X(I3, 1) - X(I2, 1)
      X13 = X(I1, 1) - X(I3, 1)
      X21 = X(I2, 1) - X(I1, 1)
      Y23 = X(I2, 2) - X(I3, 2)
      Y31 = X(I3, 2) - X(I1, 2)
      Y12 = X(I1, 2) - X(I2, 2)
      DETJ = X13 * Y23 - X32 * Y31
      BT(1, 1) = Y23
      BT(1, 2) = Y31
      BT(1, 3) = Y12
      BT(2, 1) = X32
      BT(2, 2) = X13
      BT(2, 3) = X21
      DO 35 II = 1,3
         DO 37 JJ = 1,2
            BT(JJ, II) = BT(JJ, II) / DETJ
37       CONTINUE
35    CONTINUE
      QX = BT(1, 1) * F(I1) + BT(1, 2) * F(I2) + BT(1, 3) * F(I3)
      QX = -QX * PM(MAT(I), 1)
      QY = BT(2, 1) * F(I1) + BT(2, 2) * F(I2) + BT(2, 3) * F(I3)
      QY = -QY * PM(MAT(I), 1)
   WRITE(LOUT, '(I5,2E15.5)') I, QX, QY
45 CONTINUE
   CLOSE(LOUT)
   PRINT *,'Complete results are in file ',FILE2
   IF (IPL .GT. 1) THEN
   CLOSE (LOUT2)
   PRINT *,'Element Stress Data in file ',FILE3
   PRINT *,'Run CONTOURA or CONTOURB to plot ISOTHERMS'
   END IF
   END

   SUBROUTINE BANSOL(NQ,NBW,IMAX,S,F)
   DIMENSION S(IMAX,1),F(1)
   N = NQ
C    ----- Forward Elimination -----
   DO 39 K=1,N-1
      NBK = N - K + 1
      IF ((N - K + 1) .GT. NBW) NBK = NBW
      DO 43 I=K+1, NBK+K-1
         I1 = I - K + 1
         C = S(K, I1) / S(K, 1)
         DO 41 J=I, NBK+K-1
            J1 = J - I + 1
            J2 = J - K + 1
            S(I, J1) = S(I, J1) - C * S(K, J2)
41       CONTINUE
         F(I) = F(I) - C * F(K)
43       CONTINUE
39 CONTINUE
C    ----- Back Substitution -----
   F(N) = F(N) / S(N, 1)
   DO 47 II=1,N-1
      I = N - II
      NBI = N - I + 1
      IF ((N - I + 1) .GT. NBW) NBI = NBW
      SUM = 0.
      DO 45 J=2,NBI
         SUM = SUM + S(I, J) * F(I + J - 1)
45    CONTINUE
      F(I) = (F(I) - SUM) / S(I, 1)
47 CONTINUE
   RETURN
   END

作者: little_kar 时间: 2017-3-8 16:34
没写完，但又不知道怎么删，好尴尬啊，就比如PRINT *,'Input Data File Name <DOS file name>'
   READ '(A)',FILE1
   LINP=10
   OPEN(UNIT=10,FILE=FILE1,STATUS='UNKNOWN')

作者: little_kar 时间: 2017-3-8 16:36
没写完，但又不知道怎么删，好尴尬啊，就比如这一段
PRINT *,'Input Data File Name <DOS file name>'
   READ '(A)',FILE1
   LINP=10
   OPEN(UNIT=10,FILE=FILE1,STATUS='UNKNOWN')
中，A是个啥，linp是干嘛的
还有后面怎么老是出现这一行“READ(LINP,'(A)')DUMMY”
我应该就贴几行的，实在抱歉，没及时找到怎么删

作者: vvt 时间: 2017-3-8 17:11
你可以直接上传文件。

作者: little_kar 时间: 2017-3-8 18:14

vvt 发表于 2017-3-8 17:11
你可以直接上传文件。

多谢提醒，我重新发

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)