6. A distance field

Chapter 5 computed the exact distance at every point. A solver usually wants the signed distance at every node of its grid, as an array, and does not need the last digit far from the body. distance_field_object builds that field on a rectilinear grid, stretched toward the far field as solver grids are.

use fossil, only : surface_stl_object, distance_field_object, SIGN_FAST_WINDING_NUMBER, STATUS_OK
use penf,   only : I4P, R8P
use vecfor, only : vector_R8P, ex_R8P, ey_R8P, ez_R8P

implicit none
type(surface_stl_object)    :: surface
type(distance_field_object) :: field
real(R8P), allocatable      :: x(:)
real(R8P)                   :: h, d, err_band, err_far
integer(I4P)                :: i, j, k, m, ni, nj, nk, sweeps, fallback, status

call surface%load_from_file(file_name='dragon-placed.stl', guess_format=.true.)
x = [(-0.6_R8P + 0.04_R8P * i, i = 0, 30)]                    ! uniform core: 31 points, spacing 0.04
h = 0.04_R8P
do m = 1, 6                                                     ! 6 points per side, spacing growing by 1.2
   h = 1.2_R8P * h
   x = [x(1) - h, x, x(size(x)) + h]
enddo
call field%initialize(x=x, y=x, z=x, status=status)

• initialize(x=, y=, z=) takes the three node coordinate arrays, strictly increasing, so each axis may be stretched. Here each axis has a uniform core of 31 points with spacing 0.04 over [-0.6, 0.6] and 6 more points per side, the spacing growing by 1.2. ng= declares halo layers for a block decomposition (default 0). A coordinate array that is not strictly increasing gives status = STATUS_INVALID_INPUT.

call field%compute(surface=surface, sign_algorithm=SIGN_FAST_WINDING_NUMBER, sweeps=sweeps, fallback_cells=fallback, &
                   status=status)
if (status /= STATUS_OK) error stop 'distance field failed'
call field%get_dims(ni=ni, nj=nj, nk=nk)
print '(A,I0,A,I0,A,I0)', 'grid         ', ni, ' x ', nj, ' x ', nk
print '(A,F8.4,A,I0,A)',  'band         ', field%get_band_width(), ' half-width, ', field%get_band_cells(), ' cells'
print '(A,I0,A,I0)',      'sweeps       ', sweeps, ', point-query signs ', fallback
print '(A,I0)',           'inside cells ', count(field%phi < 0._R8P)

• compute(surface=) runs three stages, also callable one by one:
  • compute_band: exact signed distances, by one distance_many batch, on the nodes within the band of the surface. band= is the half-width; the default is twice the largest grid spacing around the surface.
  • compute_signs: the signs outside the band, by a flood fill over the components of non-band nodes. A component whose sign cannot be read from the band gets one point query per node; fallback_cells counts them.
  • sweep: the magnitudes outside the band by the Fast Sweeping Method, which solves |∇φ| = 1 outward from the band. Band values stay fixed. tol (default 1e-12, relative to the domain diagonal) and max_sweeps (default 400) stop it; sweeps returns the directional sweeps done.
• sign_algorithm= is passed to the band and point queries (default SIGN_PSEUDO_NORMAL); the open dragon needs SIGN_FAST_WINDING_NUMBER (chapter 5).
• The result is the public array field%phi(i, j, k), the signed distance at (x(i), y(j), z(k)), negative inside. field%in_band marks the band nodes. get_dims, get_band_width and get_band_cells describe the field.

err_band = 0._R8P ; err_far = 0._R8P
do k = 1, nk ; do j = 1, nj ; do i = 1, ni
   d = surface%distance(point=x(i) * ex_R8P + x(j) * ey_R8P + x(k) * ez_R8P, is_signed=.true., is_square_root=.true., &
                        sign_algorithm=SIGN_FAST_WINDING_NUMBER)
   if (field%in_band(i, j, k)) then
      err_band = max(err_band, abs(field%phi(i, j, k) - d))
   else
      err_far = max(err_far, abs(field%phi(i, j, k) - d))
   endif
enddo ; enddo ; enddo
print '(A,ES10.3)', 'max error in the band ', err_band
print '(A,ES10.3)', 'max error elsewhere   ', err_far

The check compares every node with the exact distance of chapter 4, separately in and out of the band.

Running it

$ tut06_distance_field
grid         43 x 43 x 43
band           0.0800 half-width, 3667 cells
sweeps       32, point-query signs 0
inside cells 849
max error in the band  0.000E+00
max error elsewhere    7.609E-02

• The band has half-width 0.08, twice the core spacing, and holds 3667 of the 79507 nodes. Its values are the exact distances, bit for bit.
• No node needed a point query for its sign, and 32 directional sweeps (4 cycles of 8 directions) converged.
• Outside the band the sweeping is a first-order approximation of the distance: the largest error is 0.076, on the far nodes of the stretched grid. A solver that needs exact values farther out widens band.
• 849 nodes are inside. With 0.04³ per node that is a volume of 0.054, the volume of the dragon from chapter 3.

Each node is drawn as its dual cell. Left, the colour is sign(φ)·√(|φ| / 0.6), so the thin inside shows; right, it is φ / 0.08, the band half-width.

What you learned

distance_field_object fills a rectilinear grid with the signed distance: exact in a band, signed by flood fill, swept beyond. Reference: Signed-distance field.

Next: 7. Curves in 2-D.