Questions on simplicial interpolation

Help
1

Dear all:

Hopefully, this is not too far away from Fortran.

Recently I have read several papers to understand the efficiency and accuracy of simplicial interpolation compared with trilinear interpolation. It turns out that in three dimensions, there are multiple ways to decompose a cube into triangles.

I start with Weiser and Zarantonello (1988), which describes the benefit of simplicial interpolation (only requires D+1 points in D dimensional problem) and recursive implementation of simplicial interpolation. Its implementation looks like this:

    function serp_3d(m, loc, wgt, vmat) result(fval)
        integer, intent(in) :: m, loc(m)
        real(wp), intent(in) :: vmat(:, :, :), wgt(m)
        integer :: xloc(m), xloc1(m), hvec(m), i
        real(wp) :: xwgt(m)
        real(wp) :: fval, f0, a0, f1, a1

        xwgt = 1.0_wp - wgt

        call insertionSort(xwgt, hvec, m)

        xloc = loc + 1 ! starting from point (1, 1, 1)

        f0 = vmat(xloc(1), xloc(2), xloc(3))
        a0 = f0
        do i = 1, m, 1
            xloc1 = xloc
            xloc1(hvec(i)) = xloc1(hvec(i)) - 1
            f1 = vmat(xloc1(1), xloc1(2), xloc1(3))
            a1 = a0 + (1.0_wp - xwgt(i))*(f1 - f0)
            xloc = xloc1
            a0 = a1
            f0 = f1
        enddo

        fval = a1

    end function serp_3d

This method traces the unit cube from (1, 1, 1) gradually walks back to (0, 0, 0), and uses the sorted index on the right-hand side weight to figure out which tetrahedron the desired point is located. In Kasson, Plouffe, and Nin (1993) and Kasson (1994), this is called “diagonal extraction”, as all the tetrahedron are separated following the x = y = z diagonal line.

A similar algorithm appears in Rovatti, Borgatti, and Guerrieri (1998), the implementation now starts from point (0, 0, 0) and gradually walks to (1, 1, 1):

    function serp_3d(d, loc, wgt, fmat) result(fval)
        integer, intent(in) :: d, loc(d)
        real(wp), intent(in) :: fmat(:, :, :), wgt(d)
        integer :: xloc(d), xloc1(d), hvec(d), i
        real(wp) :: xwgt(d)
        real(wp) :: fval, f0, a0, f1, a1

        xwgt = 1.0_wp - wgt
        xloc = loc

        call insertionSort_reverse(xwgt, hvec, 3)

        fval = fmat(xloc(1), xloc(2), xloc(3))*(1.0_wp - xwgt(1))
        xloc(hvec(1)) = xloc(hvec(1)) + 1
        do i = 1, d-1, 1
            fval = fval + fmat(xloc(1), xloc(2), xloc(3))*(xwgt(i) - xwgt(i+1))
            xloc(hvec(i+1)) = xloc(hvec(i+1)) + 1
        enddo
        fval = fval + fmat(xloc(1), xloc(2), xloc(3))*xwgt(d)

    end function serp_3d

I tested both with a function

pure function frost(x,y,z)
real(wp):: x, y, z, frost
intent(in):: x, y, z
frost = -1.0_wp*x + dlog(y) + 2.0_wp*(y+z)**2.0_wp
end function frost

and they give the same error compared with the true value, so I think both implementations are the same.

What’s problematic for me is the accuracy of simplicial interpolation. Here is a comparison with linear interpolation:

 absolute error with cube:    8.6531889356791680E-002
 absolute error with simplex:   0.14751192565441329

And here is the elapsed time between cube and simplex over 1000 randomly generated interpolating points:

  error from   1000 approximations is     0.092995968
      recursive elapsed time:         0.0000358630 seconds.
  error from   1000 approximations is     0.132277268
simplex interpolation elapsed      0.0000864090 seconds.
             Extra error caused by serp compared with lerp:        42.2397886346%

Somehow it is slower and less accurate in 3D.

Will there be any other way to bisect a cube into tetrahedrons that give higher accuracy and less elapsed time? If so, how to implement them in Fortran?

Thank you so much!

Hui-Jun Chen