趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

トリリニア補間

トリリニア補間

公開しよう思って書いていたのですが,あまりにもソースが長くなってしまったので,少し綺麗に修正していました.それでも長いですが.

これってかなり一般的な形状にも対応しているため複雑ですが,立方体の3辺がxyz軸に沿っている場合,数行で書けることに途中で気付きました.
例えば,Trilinear interpolation - Wikipedia, the free encyclopedia
の場合です.IB法の計算ではこれで十分だったのでせっかく書いたソースは不要でした.

Fortran2003

今回構造体を使ったのですが,f90だとメソッドを構造体に入れられないことが分かりました.可読性の観点からも残念です.
Object-oriented programming in Fortran Wiki
を見ていると2003を使いたくなりますね.ただ,gfortran,g95ともにまだ対応してないので今回は断念します.


! pseudo-code from NASAを基に作成
!http://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/trilinear.html

MODULE Class_Trilinear
    IMPLICIT NONE
    PRIVATE
    PUBLIC :: Trilinear_ShapeCoff

    TYPE, PUBLIC :: shape
        REAL(8),PUBLIC :: x,y,z
    END TYPE shape
CONTAINS
    SUBROUTINE Trilinear_ShapeCoff(c_, xans_, yans_, zans_, bata_)
        REAL(8),INTENT(IN) :: xans_,yans_,zans_
        TYPE(shape),INTENT(IN) :: c_(1:8)
        REAL(8),INTENT(OUT) :: bata_(1:8)
        REAL(8) :: fmat(1:3),dfmat(1:3,1:3),rdfmat(1:3,1:3),amat(1:3),damat(1:3), x(1:8), fx(0:7), fy(0:7), fz(0:7)
        INTEGER(4) :: n
!Set_fdata
        x(1:8) = c_(1:8)%x
        CALL Set_felement(xans_,fx)
        x(1:8) = c_(1:8)%y
        CALL Set_felement(yans_,fy)
        x(1:8) = c_(1:8)%z
        CALL Set_felement(zans_,fz)
!High order Newton Method
        amat(1:3)=0.0
        DO n = 1, 100
!Set_matrix
            CALL Set_matrix(fx,1)
            CALL Set_matrix(fy,2)
            CALL Set_matrix(fz,3)
!Cal revers matrix
            CALL Rev_mat3
            damat = MATMUL(rdfmat , fmat)
            PRINT* , n, amat(1:3), MAXVAL(damat(1:3))
        IF(MAXVAL(damat(1:3)) < 0.00001 ) EXIT
            amat(1:3) = amat(1:3) - damat(1:3)
        END DO

        bata_(1) = (1-amat(3))*(1-amat(2))*(1-amat(1))
        bata_(2) = (1-amat(3))*(1-amat(2))*(1+amat(1))
        bata_(3) = (1-amat(3))*(1+amat(2))*(1-amat(1))
        bata_(4) = (1-amat(3))*(1+amat(2))*(1+amat(1))
        bata_(5) = (1+amat(3))*(1-amat(2))*(1-amat(1))
        bata_(6) = (1+amat(3))*(1-amat(2))*(1+amat(1))
        bata_(7) = (1+amat(3))*(1+amat(2))*(1-amat(1))
        bata_(8) = (1+amat(3))*(1+amat(2))*(1+amat(1))

    CONTAINS
!----------------------------------------------------------------------
        SUBROUTINE Set_felement(x0,f)
            REAL(8),INTENT(IN) :: x0
            REAL(8),INTENT(OUT) :: f(0:7)
            f(0) = (x(8)+x(7)+x(6)+x(5)+x(4)+x(3)+x(2)+x(1))/8-x0
            f(1) = (x(8)-x(7)+x(6)-x(5)+x(4)-x(3)+x(2)-x(1))/8
            f(2) = (x(8)+x(7)-x(6)-x(5)+x(4)+x(3)-x(2)-x(1))/8
            f(3) = (x(8)+x(7)+x(6)+x(5)-x(4)-x(3)-x(2)-x(1))/8
            f(4) = (x(8)-x(7)-x(6)+x(5)+x(4)-x(3)-x(2)+x(1))/8
            f(5) = (x(8)-x(7)+x(6)-x(5)-x(4)+x(3)-x(2)+x(1))/8
            f(6) = (x(8)+x(7)-x(6)-x(5)-x(4)-x(3)+x(2)+x(1))/8
            f(7) = (x(8)-x(7)-x(6)+x(5)-x(4)+x(3)+x(2)-x(1))/8
        END SUBROUTINE Set_felement
!----------------------------------------------------------------------
        SUBROUTINE Set_matrix(f,itmp)
            REAL(8),INTENT(IN) :: f(0:7)
            INTEGER(4),INTENT(IN) :: itmp
            fmat(itmp) = f(0) + f(1)*amat(1) + f(2)*amat(2) + f(3)*amat(3)&
                       + f(4)*amat(1)*amat(2) + f(5)*amat(1)*amat(3) + f(6)*amat(2)*amat(3)&
                       + f(7)*amat(1)*amat(2)*amat(3)

            dfmat(itmp, 1) = f(1) + f(4)*amat(2) + f(5)*amat(3) + f(7)*amat(2)*amat(3)
            dfmat(itmp, 2) = f(2) + f(4)*amat(1) + f(6)*amat(3) + f(7)*amat(1)*amat(3)
            dfmat(itmp, 3) = f(3) + f(5)*amat(1) + f(6)*amat(2) + f(7)*amat(1)*amat(2)
        END SUBROUTINE Set_matrix
!----------------------------------------------------------------------
        SUBROUTINE Rev_mat3
            REAL(8) :: det
            det = dfmat(1,1)*dfmat(2,2)*dfmat(3,3)+dfmat(2,1)*dfmat(3,2)*dfmat(1,3)+dfmat(3,1)*dfmat(1,2)*dfmat(2,3)&
                 -dfmat(1,1)*dfmat(3,2)*dfmat(2,3)-dfmat(3,1)*dfmat(2,2)*dfmat(1,3)-dfmat(2,1)*dfmat(1,2)*dfmat(3,3)
            IF(det == 0 )STOP 'error reverse matrix'
            det = 1.0/det
            rdfmat(1,1) = det*(dfmat(2,2)*dfmat(3,3)-dfmat(2,3)*dfmat(3,2))
            rdfmat(1,2) = det*(dfmat(1,3)*dfmat(3,2)-dfmat(1,2)*dfmat(3,3))
            rdfmat(1,3) = det*(dfmat(1,2)*dfmat(2,3)-dfmat(1,3)*dfmat(2,2))
            rdfmat(2,1) = det*(dfmat(2,3)*dfmat(3,1)-dfmat(2,1)*dfmat(3,3))
            rdfmat(2,2) = det*(dfmat(1,1)*dfmat(3,3)-dfmat(1,3)*dfmat(3,1))
            rdfmat(2,3) = det*(dfmat(1,3)*dfmat(2,1)-dfmat(1,1)*dfmat(2,3))
            rdfmat(3,1) = det*(dfmat(2,1)*dfmat(3,2)-dfmat(2,2)*dfmat(3,1))
            rdfmat(3,2) = det*(dfmat(1,2)*dfmat(3,1)-dfmat(1,1)*dfmat(3,2))
            rdfmat(3,3) = det*(dfmat(1,1)*dfmat(2,2)-dfmat(1,2)*dfmat(2,1))
        END SUBROUTINE Rev_mat3
    END SUBROUTINE Trilinear_ShapeCoff
END MODULE Class_Trilinear



PROGRAM Main
    USE Class_Trilinear
    IMPLICIT NONE
    INTEGER(4) :: i,j,k,n
    REAL(8), DIMENSION(0:1,0:1,0:1) :: x,y,z,v
    REAL(8) :: xans, yans, zans, vans, beta(1:8)
    TYPE(shape) :: c(1:8)

    i = 0; j = 0; k = 0

    x(i,j,k)=0 ; x(i+1,j,k)=1 ; x(i,j+1,k)=0 ; x(i+1,j+1,k)=1
    x(i,j,k+1)=0 ; x(i+1,j,k+1)=1 ; x(i,j+1,k+1)=0 ; x(i+1,j+1,k+1)=1

    y(i,j,k)=0 ; y(i+1,j,k)=0 ; y(i,j+1,k)=1 ; y(i+1,j+1,k)=1
    y(i,j,k+1)=0 ; y(i+1,j,k+1)=0 ; y(i,j+1,k+1)=1 ; y(i+1,j+1,k+1)=1

    z(i,j,k)=0 ; z(i+1,j,k)=0 ; z(i,j+1,k)=0 ; z(i+1,j+1,k)=0
    z(i,j,k+1)=1 ; z(i+1,j,k+1)=1 ; z(i,j+1,k+1)=1 ; z(i+1,j+1,k+1)=1

    v(i,j,k)=10 ; v(i+1,j,k)=10 ; v(i,j+1,k)=10 ; v(i+1,j+1,k)=10
    v(i,j,k+1)=20 ; v(i+1,j,k+1)=20 ; v(i,j+1,k+1)=20 ; v(i+1,j+1,k+1)=20

!Set data
    n = 1
    DO k = 0, 1 ; DO j = 0, 1 ; DO i = 0, 1
        c(n)%x = x(i,j,k) ; c(n)%y = y(i,j,k) ; c(n)%z = z(i,j,k)
        n=n+1
    END DO ; END DO ; END DO

!Cal_Trilinear Interpalation
    xans=0.1 ; yans=0.6 ; zans=0.2
    CALL Trilinear_ShapeCoff(c, xans, yans, zans, beta)

    vans = 0.0
    n = 1
    DO k = 0, 1 ; DO j = 0, 1 ; DO i = 0, 1
        vans = vans + v(i,j,k)*beta(n)
        n=n+1
    END DO ; END DO ; END DO
    vans = vans/8.0
    PRINT* , vans

END PROGRAM Main