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-сборке) здесь |
