トリリニア補間
公開しよう思って書いていたのですが，あまりにもソースが長くなってしまったので，少し綺麗に修正していました．それでも長いですが．
これってかなり一般的な形状にも対応しているため複雑ですが，立方体の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
