Поиск:

Ответ в темуСоздание новой темы Создание опроса
> [General] Конвертер из Fortran в C, Помогите разобраться 
:(
    Опции темы
Aleskert
Дата 21.4.2009, 15:58 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 14
Регистрация: 18.10.2008

Репутация: нет
Всего: нет



Доброго времени суток.
Я прошу у Вас помощи. Дело в том, что не могу разобраться как конвертировать из Fortran в C с помощью программы f2c. Было бы неплохо посмотреть step-by-step.
Если вы имеете опыт в работе с этой программой или более-менее разбираетесь, прошу помочь.

Если Вам не составит труда: гляньте код программы
Суть проги: рассчитывает поток воды в криволинейных областях
Код

1
2    C
3    C    LAGEN APPLIES THE FINITE DIFFERENCE METHOD TO LAPLACES
4    C    EQUATION IN GENERALISED COORDINATES ON A HODIFIED POLAR GRID
5    C    THE DISCRETISED EQUATION IS SOLVED BY SOR
6    C
7    DIMENSION XG(23,23),YG(23,23),GWW(21,21),GWT(21,21),GTT(21,2l),
8    lDELZI(21,21),DELET(21,21),PHI(21,21),PHIX(21,21)
9    COMMON /GRIDP/XG,YG,PHIX,PHI
10    COMMON /TRAPP/GVV,GVT,GTT,DELZI,DEL*ET
11    C
12    OPEN(1,FILE='LAGEN.DAT')
13    OPEN(6,FILE='LAGEN.OUT')
14    READ(1,1)JMAX,KMAX,NMAX,IEX,IPR
15    READ(l,2)RV,RX,RY,RZ,THEB,THEN,EPS,OM
16    READ(l,2)BXI,BXXI,CXI,CXXI
17    1 FORMAT(8I5)
18    2 FORMAT(8E10.3)
19    C
20    WRITE(6,3)
21    WRITE(6,4)JMAX,KMAX,NMAX,IEX,EPS,OK 
22    WRITE(6,5)RW,RX,RY,RZ,THEB,THEN
23    3 FORMAT(' LAPLACE EQUATION BY GEN. COORD. FDM',//)
24    4 FORMAT(' JMAX',M2,' KMAX=',I2,' NMAX=',I5,' IEX-',I2,
25    1'     EPS=',E10.3,'     OM=',F5.3)
26    5 FORMAT('     RW=',F5.3,'     RX=',F5.3,'     RY=',F5.3,'     RZ=',F5.3.
27    15X,'   THEB=',F5.1,'   THEN=',F5.1,//)
28    AKM = KMAX - 2
29    AJM = JMAX - 2
30    C
31    C    SET UP GRID
32    C
33    CALL GRID(JMAX,KMAX,THEB,THEN,RW,RX,RY,RZ)
34    C
35    C    SET BOUNDARY VALUES OF PHI
36    C
37    DO 8 J =1,JMAX
38    PHI(J,1) = 0.
39    PHI(J,KMAX) = PHIX(J,KMAX)
40    8 CONTINUE
41    DO 9 K = 1,KMAX
42    PHI(1,K) = PHIX(1,K)
43    PHI(JMAX,K) s PHIX(JMAX1K)
44    9 CONTINUE
45    C
46    C    SET GRID RELATED PARAMETERS
47    C
48    CALL TRAPA(JMAX,KMAX)
49    C
50    C    ITERATE USING SOR
51    C
52    CALL ITER(JMAX,KMAX,NMAX,N,OM,EPS)
53    C
54    C    COMPARE SOLUTION VITH EXACT
55    C
56    10 SUM = 0.
57    DO 15 K = 1,KMAX
58    WRITE(6,11)K
59    11 FORMAT(/,' K=',12)
60    DO 12 J = 1,JMAX
61    DIF = PHI(J,K) - PHIX(J,K)
62    SUM = SUM + DIF*DIF
63    12 CONTINUE
64    WRITE(6,13)(PHI(J,K),J=1,JMAX)
65    WRITE(6,14)(PHIX(J,K),J=1,JMAX)
66    13 FORMAT(' PHI=',10F7.4)
67    14 FORMAT(' PHX=',10F7.4)
68    15 CONTINUE
69    RMS = SQRT(SUM/AJM/AKM)
70    WRITE(6,16)N,RMS
71    16 FORMAT(/,' CONVERGED AFTER ',13,' STEPS,     RMS=',E12.5)
72    17 CONTINUE
73    STOP
74    END

Код

1
2    SUBROUTINE GRID(JMAX,KMAX,THEB,THEN,RW,RX,RY,RZ)
3    C
4    C    SET THE AUGMENTED GRID,  INITIAL AND EXACT PHI
5    C
6    DIMENSION XG(23.23),YG(23.23),PHIX(21,21),PHI(21,21)
7    COMMON /GRIDP/XG,YG,PHIX,PHI
8    JMAP = JMAX - 1
9    KMAP = KMAX - 1
10    AJM = JMAP
11    AKM = KMAP
12    DRWX = (RX - RW)/AJM
13    DRZY = (RY - RZ)/AKM
14    DTH = (THEN-THEB)/AKM
15    PI = 3.1415927
16    KPP = KMAX + 2
17    PP = JMAX + 2
18    C
19    C    SE; XG, YG. EXACT AKD INITIAL PHI
20    C
21    DO 7 K = 1,KPP
22    AK = K - 2
23    THK = (THEB + AK*DTH)*PI/180.
24    CK = COS(THK)
25    SK = SIN(THK)
26    DR = DRWX + (DRZY - DRWX)*AK/AKM
27    RWZ = RW + (RZ - RW)*AK/AKM
28    DO 6 J = 1,JPP
29    AJ = J - 2
30    R = RWZ + AJ*DR
31    XG(J,K) = R*CK
32    YG(J,K) = R*SK
33    IF(K .EQ. 1 .OR. K .EQ. KPP)GOTO 6
34    IF(J .EQ. 1 .OR. J .EQ. JPP)GOTO 6
35    JM = J-1
36    KM = K - 1
37    PHIX(JM,KM) = SK/R
38    PHI(JM,KM) = PHIX(JM,KM)
39    6 CONTINUE
40    7 COTTTINUE
41    RETURN
42    END

Код

1
2    SUBROUTINE TRAPA(JMAX,KMAX)
3    C
4    C    FROM GRID COORDINATES CALCULATES TRANSFORM PARAMETERS
5    C
€    DIMENSION XG(23,23),YG(23,23),PHIX(21,21),PHI(21,21)
7    DIMENSION GWW(21,21),GWT(21,21),GTT(21,21),DELZI(21,21),
8    1DELET(21,21)
9    COMMON /GRIDP/XG,YG,PHIX,PHI
10    COMMON /TRAPP/GWW,GWT,GTT,DELZI,DELET
11    DO 2 K = 1,KMAX
12    KP = K+l
13    KPP = K + 2
14    DO 1 J = 1,JMAX
15    JP = J + 1
16    JPP = J + 2
17    C
18    C    BASIC TRANSFORM PARAMETERS
19    C
20    XZI = 0.5*(XG(JPP,KP) - XG(J,KP))
21    YZI = 0.5*(YG(JPP,KP) - YG(J,KP))
22    XET = 0.5*(XG(JP,KPP) - XG(JP,K))
23    YET = 0.5*(YG(JP,KPP) - YG(JP,K))
24    XZZ = XG(J,KP) - 2.*XG(JP,KP) + XG(JPP,KP)
25    YZZ = YG(J,KP) - 2.*YG(JP,KP) + YG(JPP,KP)
26    XEE = XG(JP,K) - 2.*XG(JP,KP) + XG(JP,KPP)
27    YEE = YG(JP,K) - 2.*YG(JP,KP) + YG(JP,KPP)
28    XZE = 0.25*(XG(JPP,KPP)-XG(J,KPP)+XG(J,K)-XG(JPP,K))
29    YZE = 0.25*(YG(JPP,KPP)-YG(J,KPP)+YG(J,K)-YG(JPP,K))
30    AJ = XZI*YET - XET*YZI
31    C
32    C    MODIFIED METRIC TENSOR COEFFICIENTS, G11,G12,G22
33    C
34    GWW(J,K) = (XZI*XZI + YZI*YZI)/AJ
35    GWT(J,K) = -2.*(XZI*XET + YZI*YET>/AJ
36    GTT(J,K) " (XET*XET + YET*YET)/AJ
37    C
38    C    MODIFIED DEL**2ZI AND DEL**2ETA
39    C
40    DUM = GTT(J,K)*(XET*YZZ-YET*XZZ)/AJ
41    DUM = DUM + GWT(J,K)*(XET*YZE-YET*XZE)/AJ
42    DELZI(J,K) = DUM + GWW(J,K)*(XET*YEE-YET*XEE)/AJ
43    DUM * GTT(J,K)*(YZI*XZZ - XZI*YZZ)/AJ
44    DUM = DUM + GVT(J,K)*(YZI*XZE-XZI*YZE)/AJ
45    DELET(J,K) = DUM + GVV(J,K)*(YZI*XEE-XZI*YEE)/AJ
46    1 CONTINUE
47    2 CONTINUE
48    RETURN
49    END

Код

1
2    SUBROUTINE ITER(JMAX,KMAX,NMAX,N,OM,EPS)
3    C
4    C    ITERATE USING SOR APPLIED TO DISCRETISED EQUATIOHS
5    C
6    DIMENSION GWW(21,21),GWT(21,21),GTT(21,21),DELZI(21,21),
7    1DELET(21.21),PHI(21,21),PHIX(21,21),XG(23,23),YG(23,23)
8    COHNON /GRIDP/XG,YG,PHIX,PHI
9    COMHON /TRAPP/GWW,GWT,GTT,DELZI,DELET
10    KMAP = KMAX-1
11    JMAP = JMAX-1
12    AKM = KMAP-1
13    AJM = JMAP-1
14    C
15    DO 3 N = l,NMAX
16    SUM = 0.
17    DO 2 K = 2,KMAP
18    KM = K - 1
19    KP = K + 1
20    DO 1 J = 2,JMAP
21    JM = J - 1
22    JP = J + 1
23    PHD = -0.5*(DELZI(JP,K)*PHI(JP,K)-DELZI(JM,K)*PHI(JM,K))
24    PHD = PHD-0.5*(DELET(J,KP)*PHI(J,KP)-DELET(J,KM)*PHI(J,KM)*
25    PHD = PHD + GTT(JM.K)*PHI(JH,K) + GTT(JP.K)*PHI(JP,K)
26    PHD = PHD + GWW(J,KM)*PHI(J,KM) + GWW(J,KP)*PHI(J,KP)
27    PHD = PHD + 0.25*(GWT(JP,KP)*PHI(JP,KP)-GWT(JM.KP)*PHI(JM,KPP
28    1 + GVT(JM,KM)*PHI(JM,KM) - GVT(JP,KH)*PHI(JP,KM))
29    PHD = 0.5*PHD/(GTT(J,K)+GWW(J,K))
30    DIF = PHD - PHI(J,K)
31    SUM = SUM + DIF*DIF
32    PHI(J,K) " PHI(J,K) + OM*DIF
33    1 CONTINUE
34    2 CONTINUE
35    RMS = SQRT(SUM/AJM/AKM)
36    IF(RMS .LT. EPS)RETURN
37    3 CONTINUE
38    WRITE(6,4)NMAX,RMS
39    4 FORMAT(' CONVERGENCE NOT ACHIEVED IN',I5,' STEPS,    RMS',
40    1E12.5)
41    RETURN
42    END

PM MAIL   Вверх
FCM
Дата 25.4.2009, 09:53 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 461
Регистрация: 30.3.2009

Репутация: 8
Всего: 9



Насчет f2c ничего не скажу.
Но в данной программе слишком много опечаток (много больше 10) + ошибки форматирования исходника (относительно стандартного фиксированного формата F77) .
Кстати, а файл 'LAGEN.DAT' у тебя есть?
PM MAIL   Вверх
Aleskert
Дата 26.4.2009, 13:41 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 14
Регистрация: 18.10.2008

Репутация: нет
Всего: нет



LAGEN.DAT - файл с исходными данными
LAGEN.OUT - файл с обработанными данными
Если Вас не затруднит, залейте сюда код на С и скомпилированный выполняемый файл.
Ибо в Фортране и f2c не смыслю ровным счетом ничего.
Заранее спасибо!
PM MAIL   Вверх
FCM
Дата 26.4.2009, 18:22 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 461
Регистрация: 30.3.2009

Репутация: 8
Всего: 9



Насчет f2c, как уже выше писал, ничего не скажу. 
Но, судя по характеру ошибок в листинге, его перепечатывали в далекие времена на машинке, и помимо явных опечаток там могут быть и скрытые недоразумения, например неправильно указанные метки (например 10 и 17). 
По моему, если знаешь алгоритм, надежнее сразу написать на С/C++ или современном фортране.





Это сообщение отредактировал(а) FCM - 26.4.2009, 19:19
PM MAIL   Вверх
Aleskert
Дата 26.4.2009, 18:42 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 14
Регистрация: 18.10.2008

Репутация: нет
Всего: нет



В том то и дело, что я в Фортране ничего не смыслю!
Я не могу понять даже что и куда там обращается - только общий смысл улавливаю.
А общий смысл мне не надо smile
PM MAIL   Вверх
FCM
Дата 26.4.2009, 19:20 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 461
Регистрация: 30.3.2009

Репутация: 8
Всего: 9



Есть другой путь: собрать из данного фортран-исходника библиотеку и затем вызывать ее из С/С++,
но опять же все упирается в опечатки исходника, надежное устранение которых требует понимания алгоритма.

А тестовые значения для LAGEN.DAT и LAGEN.OUT имеются?

Это сообщение отредактировал(а) FCM - 26.4.2009, 19:22
PM MAIL   Вверх
Aleskert
Дата 28.4.2009, 23:56 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 14
Регистрация: 18.10.2008

Репутация: нет
Всего: нет



LAGEN.OUT
Код

LAPLACE EQUATION BY GEN. COORD. FDM

JMAX= 6 KMAX= 6 NMAX=100 IEX= 0 EPS= .100E-04 OM=1.500
RW= .100 RX= 1.000 RY= 1.000 RZ= .100 THEB= .0 THEN= 90.0

K= 1
PHI= .0000 .0000 .0000 .0000 .0000 .0000  
PHX= .0000 .0000 .0000 .0000 .0000 .0000

K= 2
PHI= 3.0902 1.2454 .7453 .5209 .3931 .3090
PHX= 3.0902 1.1036 .6718 .4828 .3768 .3090

K= 3
PHI= 5.8779 2.3498 1.4048 .9834 .7444 .5878
PHX= 5.8779 2.0992 1.2778 .9184 .7168 .5878

K= 4
PHI= 8.0902 3.1813 1.9002 1.3351 1.0166 .8090
PHX= 8.0902 2.8893 1.7587 1.2641 .9866 .8090

K= 5
PHI= 9.5106 3.6228 2.1688 1.5352 1.1804 .9511
PHX= 9.5106 3.3966 2.0675 1.4860 1.1598 .9511

K= 6
PHI= 10.0000 3.5714 2.1739 1.5625 1.2195 1.0000
PHX= 10.0000 3.5714 2.1739 1.5625 1.2195 1.0000

CONVERGED AFTER 15 STEPS, RMS= .13384E+00


LAGEN.DAT
Код

6 
6 
100 
0 
.100E-04 
1.500
.100 
1.000 
1.000 
.100 
.0 
90.0


PM MAIL   Вверх
FCM
Дата 29.4.2009, 19:06 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 461
Регистрация: 30.3.2009

Репутация: 8
Всего: 9



Файл LAGEN.DAT действительно такой короткий?
ПРоверь на всякий случай.

PM MAIL   Вверх
Aleskert
Дата 30.4.2009, 13:29 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 14
Регистрация: 18.10.2008

Репутация: нет
Всего: нет



Да, в этом файле находятся исходные значения переменных, который используются в программе.
Они потом дублируются в файле результатов в 3 и 4 строчке, только уже с объяснениями что какая переменная делает.
PM MAIL   Вверх
FCM
Дата 30.4.2009, 15:03 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 461
Регистрация: 30.3.2009

Репутация: 8
Всего: 9



Фортран-программа заработала.
В принципе далее можно попробовать сделать из нее библиотеку и вызывать ее из C++ программы, если конечно не нужно вмешиваться в алгоритм.



Это сообщение отредактировал(а) FCM - 30.4.2009, 18:14
PM MAIL   Вверх
Aleskert
Дата 2.5.2009, 14:31 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 14
Регистрация: 18.10.2008

Репутация: нет
Всего: нет



Можно код в С++ ?
Или экзешку выложить?
PM MAIL   Вверх
FCM
Дата 2.5.2009, 18:52 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 461
Регистрация: 30.3.2009

Репутация: 8
Всего: 9



1)Кода в С++ у меня нет. Но выкладываю подправленный вариант фортран-программы (преобразованной также к стандарту Fortran-90), по которому сам можешь сделать С++ аналог (не так это и трудно в данном случае - нужно просто распечать нижеприведенные исходники и постепенно картина станет ясной).

Код

!    LAGEN APPLIES THE FINITE DIFFERENCE METHOD TO LAPLACES
!    EQUATION IN GENERALISED COORDINATES ON A HODIFIED POLAR GRID
!    THE DISCRETISED EQUATION IS SOLVED BY SOR
PROGRAM LAGEN
    
    USE DATA_GRIDP  ! contains defs of XG, YG, PHI, PHIX
    USE DATA_TRAPP  ! contains defs of GWW, GWT, GTT, DELZI, DELET
    IMPLICIT NONE
    INTEGER :: JMAX, KMAX, NMAX, IEX, &
               N, K, J
    REAL*4  :: EPS, OM, RW, RX, RY, RZ, THEB, THEM,  &
               AKM, AJM, SUM, DIF, RMS
    OPEN(11, FILE = 'LAGEN.DAT') 
    OPEN(12, FILE = 'LAGEN.OUT')
    READ(11,*) JMAX, KMAX, NMAX, IEX, EPS, OM   
    READ(11,*) RW, RX, RY, RZ, THEB, THEM
!
    WRITE(12,101)
    WRITE(12,102) JMAX, KMAX, NMAX, IEX, EPS, OM 
    WRITE(12,103) RW, RX, RY, RZ, THEB, THEM
101 FORMAT(' LAPLACE EQUATION BY GEN. COORD. FDM',//)
102 FORMAT(' JMAX',I2,' KMAX=',I2,' NMAX=',I5,' IEX=',I2,  &
           '     EPS=', E10.3,'     OM=', F5.3 )
103 FORMAT('     RW=', F5.3,'     RX=',F5.3,'     RY=',F5.3, &
           '    RZ=',F5.3, 5X,'   THEB=',F5.1,'   THEM=',F5.1,//)
    ! Allocation of arrays defined in DATA_GRIDP module
    ALLOCATE( XG(JMAX+2,KMAX+2), YG(JMAX+2,KMAX+2),  &
              PHI(JMAX,KMAX), PHIX(JMAX,KMAX) )
    ! Allocation of arrays defined in DATA_TRAPP module
    ALLOCATE( GWW(JMAX,KMAX), GWT(JMAX,KMAX), GTT(JMAX,KMAX), &
              DELZI(JMAX,KMAX), DELET(JMAX,KMAX) )
           
    AKM = KMAX - 2
    AJM = JMAX - 2
        
!    SET UP GRID

    CALL GRID( JMAX, KMAX, THEB, THEM, RW, RX, RY, RZ )

!    SET BOUNDARY VALUES OF PHI

    PHI(1:JMAX,1) = 0.
    PHI(1:JMAX,KMAX) = PHIX(1:JMAX,KMAX)
    PHI(1,1:KMAX) = PHIX(1,1:KMAX)
    PHI(JMAX,1:KMAX) = PHIX(JMAX,1:KMAX)
    
!    SET GRID RELATED PARAMETERS

    CALL TRAPA( JMAX, KMAX )

!    ITERATE USING SOR

    CALL ITER( JMAX, KMAX, NMAX, N, OM, EPS )

!   COMPARE SOLUTION VITH EXACT

    SUM = 0.
    DO  K = 1, KMAX
        WRITE(12,104) K
        DO  J = 1, JMAX
              DIF = PHI(J,K) - PHIX(J,K)
              SUM = SUM + DIF*DIF
        END DO
        WRITE(12,105) PHI(1:JMAX,K)  !(PHI(J,K), J=1,JMAX)
        WRITE(12,106) PHIX(1:JMAX,K) !(PHIX(J,K), J=1,JMAX)
    END DO
    RMS = SQRT( SUM/AJM/AKM ) 
    WRITE(12,107) N, RMS
    
104 FORMAT(/,' K=',I2)
105 FORMAT(' PHI =', 10F10.5)
106 FORMAT(' PHIX=', 10F10.5)
107 FORMAT(/,' CONVERGED AFTER ',I3,' STEPS,     RMS=',E12.5)

    ! Deallocation of arrays 
    DEALLOCATE( XG, YG, PHI, PHIX )
    DEALLOCATE( GWW, GWT, GTT, DELZI, DELET )

END PROGRAM LAGEN


где

Код

MODULE DATA_GRIDP
    REAL*4, ALLOCATABLE ::  XG(:,:), YG(:,:), PHI(:,:), PHIX(:,:)
END MODULE DATA_GRIDP


Код

MODULE DATA_TRAPP
    REAL*4, ALLOCATABLE ::  GWW(:,:), GWT(:,:),  GTT(:,:), DELZI(:,:), DELET(:,:)
END MODULE DATA_TRAPP


Код

SUBROUTINE GRID( JMAX, KMAX, THEB, THEM, RW, RX, RY, RZ )
!    SET THE AUGMENTED GRID,  INITIAL AND EXACT PHI

    USE DATA_GRIDP  ! contains defs of XG, YG, PHI, PHIX
    IMPLICIT NONE
    ! Parameters of procedure
    INTEGER, INTENT(IN) ::  JMAX, KMAX
    REAL*4,  INTENT(IN) ::  THEB, THEM, RW, RX, RY, RZ
    ! Local vars
    INTEGER :: K, KM, KPP, KMAP, J, JM, JPP, JMAP
    REAL*4  :: AK, AKM, AJ, AJM, DRWX, DRZY, DR, RWZ, R, DTH, PI, THK, CK, SK 
    
    JMAP = JMAX - 1
    KMAP = KMAX - 1
    AJM = JMAP
    AKM = KMAP
    DRWX = (RX - RW)/AJM
    DRZY = (RY - RZ)/AKM
    DTH = (THEM-THEB)/AKM
    PI = 3.1415927
    KPP = KMAX + 2
    JPP = JMAX + 2
!
!    SE; XG, YG. EXACT AKD INITIAL PHI
!
    DO  K = 1,KPP
      
        AK = K - 2
        THK = (THEB + AK*DTH)*PI/180.
        CK = COS(THK)
        SK = SIN(THK)
        DR = DRWX + (DRZY - DRWX)*AK/AKM
        RWZ = RW + (RZ - RW)*AK/AKM
        
        DO  J = 1,JPP
          
            AJ = J - 2
            R = RWZ + AJ*DR
            XG(J,K) = R*CK
            YG(J,K) = R*SK
            IF(K .EQ. 1 .OR. K .EQ. KPP) CYCLE   !GO TO NEXT J-iteration
            IF(J .EQ. 1 .OR. J .EQ. JPP) CYCLE   !GO TO NEXT J-iteration
            JM = J - 1
            KM = K - 1
            PHIX(JM,KM) = SK/R
            PHI(JM,KM) = PHIX(JM,KM)
            
         END DO
    
     END DO
    
END SUBROUTINE GRID


Код

SUBROUTINE ITER( JMAX, KMAX, NMAX, N, OM, EPS )

!      ITERATE USING SOR APPLIED TO DISCRETISED EQUATIOHS

    USE DATA_GRIDP  ! contains defs of XG, YG, PHI, PHIX
    USE DATA_TRAPP  ! contains defs of GWW, GWT, GTT, DELZI, DELET
    IMPLICIT NONE
    ! Parameters of procedure
    INTEGER, INTENT(IN)  :: JMAX, KMAX, NMAX
    INTEGER, INTENT(OUT) :: N
    REAL*4,  INTENT(IN)  :: OM, EPS
    ! Local vars
    INTEGER :: K, KM, KP, KMAP, J, JM, JP, JMAP
    REAL*4  :: AKM, AJM, SUM, DIF, PHD, RMS
    
    KMAP = KMAX-1
    JMAP = JMAX-1
    AKM = KMAP-1
    AJM = JMAP-1

    DO  N = 1, NMAX
        SUM = 0.
        DO  K = 2, KMAP
            KM = K - 1
            KP = K + 1
            DO  J = 2, JMAP
                JM = J - 1
                JP = J + 1
                PHD = -0.5*( DELZI(JP,K)*PHI(JP,K) -   &
                             DELZI(JM,K)*PHI(JM,K) )
                PHD = PHD - 0.5*( DELET(J,KP)*PHI(J,KP) -   &
                                  DELET(J,KM)*PHI(J,KM) )
                PHD = PHD + GTT(JM,K)*PHI(JM,K) + GTT(JP,K)*PHI(JP,K)
                PHD = PHD + GWW(J,KM)*PHI(J,KM) + GWW(J,KP)*PHI(J,KP)
                PHD = PHD + 0.25*   &
                         ( GWT(JP,KP)*PHI(JP,KP) - GWT(JM,KP)*PHI(JM,KP) + &
                           GWT(JM,KM)*PHI(JM,KM) - GWT(JP,KM)*PHI(JP,KM) )
                PHD = 0.5*PHD/(GTT(J,K)+GWW(J,K))
                DIF = PHD - PHI(J,K)
                SUM = SUM + DIF*DIF
                PHI(J,K) = PHI(J,K) + OM*DIF
            END DO
        END DO
        RMS = SQRT(SUM/AJM/AKM)
        IF(RMS < EPS) RETURN
    END DO
    
    WRITE(6,108) NMAX,RMS
108 FORMAT(' CONVERGENCE NOT ACHIEVED IN',I5,' STEPS,    RMS', E12.5)
    RETURN
END SUBROUTINE ITER



Код

SUBROUTINE TRAPA(JMAX,KMAX)

!    FROM GRID COORDINATES CALCULATES TRANSFORM PARAMETERS

    USE DATA_GRIDP  ! contains defs of XG, YG, PHI, PHIX
    USE DATA_TRAPP  ! contains defs of GWW, GWT, GTT, DELZI, DELET
    IMPLICIT NONE
    ! Parameters of procedure
    INTEGER, INTENT(IN) ::  JMAX, KMAX
    ! Local vars
    INTEGER ::  J, JP, JPP, K, KP, KPP
    REAL*4  ::  XZI, YZI, XET, YET, XZZ, YZZ, XEE, YEE, XZE, YZE, AJ, DUM 
    
    DO K = 1,KMAX
        KP = K + 1
        KPP = K + 2
        DO J = 1,JMAX
            JP = J + 1
            JPP = J + 2

!   BASIC TRANSFORM PARAMETERS

            XZI = 0.5*(XG(JPP,KP) - XG(J,KP))
            YZI = 0.5*(YG(JPP,KP) - YG(J,KP))
            XET = 0.5*(XG(JP,KPP) - XG(JP,K))
            YET = 0.5*(YG(JP,KPP) - YG(JP,K))
            XZZ = XG(J,KP) - 2.*XG(JP,KP) + XG(JPP,KP)
            YZZ = YG(J,KP) - 2.*YG(JP,KP) + YG(JPP,KP)
            XEE = XG(JP,K) - 2.*XG(JP,KP) + XG(JP,KPP)
            YEE = YG(JP,K) - 2.*YG(JP,KP) + YG(JP,KPP)
            XZE = 0.25*(XG(JPP,KPP)-XG(J,KPP)+XG(J,K)-XG(JPP,K))
            YZE = 0.25*(YG(JPP,KPP)-YG(J,KPP)+YG(J,K)-YG(JPP,K))
            AJ = XZI*YET - XET*YZI

!    MODIFIED METRIC TENSOR COEFFICIENTS, G11,G12,G22

            GWW(J,K) = (XZI*XZI + YZI*YZI)/AJ
            GWT(J,K) = -2.*(XZI*XET + YZI*YET)/AJ
            GTT(J,K) = (XET*XET + YET*YET)/AJ

!    MODIFIED DEL**2ZI AND DEL**2ETA

            DUM = GTT(J,K)*(XET*YZZ-YET*XZZ)/AJ
            DUM = DUM + GWT(J,K)*(XET*YZE-YET*XZE)/AJ
            DELZI(J,K) = DUM + GWW(J,K)*(XET*YEE-YET*XEE)/AJ
            DUM = GTT(J,K)*(YZI*XZZ - XZI*YZZ)/AJ
            DUM = DUM + GWT(J,K)*(YZI*XZE-XZI*YZE)/AJ
            DELET(J,K) = DUM + GWW(J,K)*(YZI*XEE-XZI*YEE)/AJ
        END DO
    END DO
    
END SUBROUTINE TRAPA



Добавлено @ 19:05
2) Другой вариант. Можно вызвать библиотеку, скомпонованную под фортраном, из С++ программ.
В данном случае это в первом приближении может выглядеть так.

Код


extern "C" void __cdecl
LAGEN_F90( const int& JMAX, const int& KMAX, const int& NMAX, const int& IEX,
           const float& EPS, const float& OM, const float& RW, 
           const float& RX, const float& RY, const float& RZ, 
           const float& THEB, const float& THEM );

int main()
{
    int    JMAX = 6, KMAX = 6, NMAX = 100, IEX = 0;
    float  EPS = 1.e-5f, OM = 1.5f, RW = 0.1f, 
           RX = 1.f, RY = 1.f, RZ = 0.1f, THEB = 0.f, THEM = 90.f;
    LAGEN_F90( JMAX, KMAX, NMAX, IEX,
               EPS, OM, RW, RX, RY, RZ, THEB, THEM );
    return 0;

}


где процедура LAGEN_F90 и используемые ею оформляется предварительно в проекте типа стат.библиотека под Intel Fortran. Отличие от вышеприведенных исходников будет только в первом листинге

Код

SUBROUTINE LAGEN_F90( JMAX, KMAX, NMAX, IEX, EPS, OM, &
                           RW, RX, RY, RZ, THEB, THEM )
    
    USE DATA_GRIDP  ! contains defs of XG, YG, PHI, PHIX
    USE DATA_TRAPP  ! contains defs of GWW, GWT, GTT, DELZI, DELET
    
    IMPLICIT NONE
    INTEGER, INTENT(IN) ::  JMAX, KMAX, NMAX, IEX
    REAL(4), INTENT(IN) ::  EPS, OM, RW, RX, RY, RZ, THEB, THEM
    INTEGER ::  N, K, J
    REAL(4) ::  AKM, AJM, SUM, DIF, RMS
    OPEN(12, FILE = 'LAGEN.OUT')
  
    WRITE(12,101)
    WRITE(12,102) JMAX, KMAX, NMAX, IEX, EPS, OM 
    WRITE(12,103) RW, RX, RY, RZ, THEB, THEM
101 FORMAT(' LAPLACE EQUATION BY GEN. COORD. FDM',//)
102 FORMAT(' JMAX',I2,' KMAX=',I2,' NMAX=',I5,' IEX=',I2,  &
           '     EPS=', E10.3,'     OM=', F5.3 )
103 FORMAT('     RW=', F5.3,'     RX=',F5.3,'     RY=',F5.3, &
           '    RZ=',F5.3, 5X,'   THEB=',F5.1,'   THEM=',F5.1,//)
    ! Allocation of arrays defined in DATA_GRIDP module
    ALLOCATE( XG(JMAX+2,KMAX+2), YG(JMAX+2,KMAX+2),  &
              PHI(JMAX,KMAX), PHIX(JMAX,KMAX) )
    ! Allocation of arrays defined in DATA_TRAPP module
    ALLOCATE( GWW(JMAX,KMAX), GWT(JMAX,KMAX), GTT(JMAX,KMAX), &
              DELZI(JMAX,KMAX), DELET(JMAX,KMAX) )
           
    AKM = KMAX - 2
    AJM = JMAX - 2
        
!    SET UP GRID

    CALL GRID( JMAX, KMAX, THEB, THEM, RW, RX, RY, RZ )

!    SET BOUNDARY VALUES OF PHI

    PHI(1:JMAX,1) = 0.
    PHI(1:JMAX,KMAX) = PHIX(1:JMAX,KMAX)
    PHI(1,1:KMAX) = PHIX(1,1:KMAX)
    PHI(JMAX,1:KMAX) = PHIX(JMAX,1:KMAX)
    
!    SET GRID RELATED PARAMETERS

    CALL TRAPA( JMAX, KMAX )

!    ITERATE USING SOR

    CALL ITER( JMAX, KMAX, NMAX, N, OM, EPS )

!   COMPARE SOLUTION VITH EXACT

    SUM = 0.
    DO  K = 1, KMAX
        WRITE(12,104) K
        DO  J = 1, JMAX
              DIF = PHI(J,K) - PHIX(J,K)
              SUM = SUM + DIF*DIF
        END DO
        WRITE(12,105) PHI(1:JMAX,K)  !(PHI(J,K), J=1,JMAX)
        WRITE(12,106) PHIX(1:JMAX,K) !(PHIX(J,K), J=1,JMAX)
    END DO
    RMS = SQRT( SUM/AJM/AKM ) 
    WRITE(12,107) N, RMS
    
104 FORMAT(/,' K=',I2)
105 FORMAT(' PHI =', 10F10.5)
106 FORMAT(' PHIX=', 10F10.5)
107 FORMAT(/,' CONVERGED AFTER ',I8,' STEPS,     RMS=',E12.5)

    ! Deallocation of arrays 
    DEALLOCATE( XG, YG, PHI, PHIX )
    DEALLOCATE( GWW, GWT, GTT, DELZI, DELET )

END SUBROUTINE LAGEN_F90

 
Вывод будет вестись в LAGEN.OUT в текущей рабочей директории для С++ программы.
Если речь идет  МVC++ или Intel C++(под Windows ), то никаких проблем со статикой не должно быть. 
В др. случаях возможно потребуется строить динамическую библиотеку.

Сама стат. библиотека (в Debug-сборке) здесь 

Это сообщение отредактировал(а) FCM - 4.5.2009, 09:37

Присоединённый файл ( Кол-во скачиваний: 2 )
Присоединённый файл  LAGEN_F90.lib 103,85 Kb
PM MAIL   Вверх
FCM
Дата 2.5.2009, 19:38 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Опытный
**


Профиль
Группа: Участник
Сообщений: 461
Регистрация: 30.3.2009

Репутация: 8
Всего: 9



ps
 по приведенному ранее тобой LAGEN.OUT не совсем понятно почему есть некоторое несовпадение по промежуточным точкам - так и должно быть ?



Это сообщение отредактировал(а) FCM - 2.5.2009, 19:44
PM MAIL   Вверх
Aleskert
Дата 4.5.2009, 22:54 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 14
Регистрация: 18.10.2008

Репутация: нет
Всего: нет



Огромное спасибо!
Несовпадение может получаться. Задача построена для потока воды в криволинейных областях. То есть могут быть "завихрения" в течении. Но как они представлены в программе - буду выяснять (еще не до конца разобрался в представлении переменных и их обработке).
Главное что программа считает все точки более-менее правильно smile
Еще раз спасибо!
PM MAIL   Вверх
Aleskert
Дата 4.5.2009, 23:45 (ссылка) | (нет голосов) Загрузка ... Загрузка ... Быстрая цитата Цитата


Новичок



Профиль
Группа: Участник
Сообщений: 14
Регистрация: 18.10.2008

Репутация: нет
Всего: нет



Огромное спасибо!
Несовпадение может получаться. Задача построена для потока воды в криволинейных областях. То есть могут быть "завихрения" в течении. Но как они представлены в программе - буду выяснять (еще не до конца разобрался в представлении переменных и их обработке).
Главное что программа считает все точки более-менее правильно smile
Еще раз спасибо!
PM MAIL   Вверх
  
Ответ в темуСоздание новой темы Создание опроса
0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема | Fortran | Следующая тема »


 




[ Время генерации скрипта: 0.0575 ]   [ Использовано запросов: 22 ]   [ GZIP включён ]


Реклама на сайте     Информационное спонсорство

 
По вопросам размещения рекламы пишите на vladimir(sobaka)vingrad.ru
Отказ от ответственности     Powered by Invision Power Board(R) 1.3 © 2003  IPS, Inc.