11. Moving bodies

The dragon now moves: over one time step it pitches by 10 degrees about y and moves 0.02 along x. Rebuilding the tree at every step is the cost to avoid. FOSSIL either queries the body in its own frame or moves a copy and refits its tree, and it gives what a moving cut-cell update needs: the cells that change state and the swept geometry.

use fossil, only : surface_stl_object, curve_object, rigid_pose_t, rigid_pose_2d_t, pose_identity, pose_from_axis_angle, &
                   pose_2d, swept_cells_2d_t, swept_cells_2d, inertia_to_world, SIGN_FAST_WINDING_NUMBER, STATUS_OK
use penf,   only : I4P, R8P
use vecfor, only : vector_R8P, ex_R8P, ey_R8P, ez_R8P

implicit none
integer(I4P), parameter       :: n = 24                        ! cell centres per direction
real(R8P),    parameter       :: pi = 4._R8P * atan(1._R8P)
type(surface_stl_object)      :: surface, moved, wrap, box
type(curve_object)            :: section
type(rigid_pose_t)            :: p0, p1
type(rigid_pose_2d_t)         :: q0, q1
type(swept_cells_2d_t)        :: sw
type(vector_R8P), allocatable :: centres(:)
real(R8P),        allocatable :: d_posed(:), d_moved(:)
logical,          allocatable :: fresh(:), dead(:)
real(R8P)                     :: h, x(49), mass, centre(3), moments(3), inertia(3,3), inertia1(3,3)
integer(I4P)                  :: i, j, k, c, checked, status, fa, fb
real(R8P)                     :: gap, gap1
logical                       :: hit

call surface%load_from_file(file_name='dragon-placed.stl', guess_format=.true.)
h = 1.2_R8P / n
allocate(centres(n**3), d_posed(n**3), d_moved(n**3), fresh(n**3), dead(n**3))
c = 0
do k = 1, n ; do j = 1, n ; do i = 1, n
   c = c + 1
   centres(c) = (-0.6_R8P + (i - 0.5_R8P) * h) * ex_R8P + (-0.6_R8P + (j - 0.5_R8P) * h) * ey_R8P + &
                (-0.6_R8P + (k - 0.5_R8P) * h) * ez_R8P
enddo ; enddo ; enddo
p0 = pose_identity()
p1 = pose_from_axis_angle(axis=[0._R8P, 1._R8P, 0._R8P], angle=pi / 18._R8P, translation=[0.02_R8P, 0._R8P, 0._R8P])

• A rigid_pose_t holds a rotation r(3, 3) and a translation t(3), with x_world = R x_body + t. pose_identity() is the rest pose, pose_from_axis_angle(axis=, angle=, translation=) builds one from an axis, an angle in radians and a translation (default 0); pose_from_quaternion, pose_compose and pose_inverse are the rest of the algebra.

call surface%distance_many(points=centres, distances=d_posed, is_signed=.true., is_square_root=.true., &
                           sign_algorithm=SIGN_FAST_WINDING_NUMBER, pose=p1)    ! the tree stays as built
moved = surface
call moved%move(p1)                                                       ! a moved copy, tree refitted
call moved%distance_many(points=centres, distances=d_moved, is_signed=.true., is_square_root=.true., &
                         sign_algorithm=SIGN_FAST_WINDING_NUMBER)
print '(A,L1)', 'posed = moved to 1e-14: ', maxval(abs(d_posed - d_moved)) < 1.e-14_R8P
call surface%state_changes(pose0=p0, pose1=p1, points=centres, fresh=fresh, dead=dead, &
                           sign_algorithm=SIGN_FAST_WINDING_NUMBER, checked=checked)
print '(A,I0,A,I0,A,I0,A,I0)', 'cells: fresh ', count(fresh), ', dead ', count(dead), ', second query ', checked, &
                               ' of ', size(centres)

• pose= on distance_many (also on query_many, is_point_inside and the device batches) maps the world points into the body frame, x̃ = Rᵀ (x − t), and queries the tree as built: one 3 × 3 product per point.
• move(pose) moves the facets of a copy and refits its BVH: node boxes recomputed bottom-up, topology kept, O(N) instead of a rebuild. Grid-aligned queries (box tagging, line crossings, cut cells) need the moved copy.
• state_changes(pose0=, pose1=, points=, fresh=, dead=) flags the points that go from solid to fluid (fresh) and from fluid to solid (dead) between two poses. A point farther from the wall than it moves cannot change state, so its second query is skipped; checked counts those that needed it. The sign algorithm defaults to the pseudo-normal, as for the point queries.

call section%load_from_file(file_name='dragon-section.dat')
x = [(-0.6_R8P + 0.025_R8P * i, i = 0, 48)]
q0 = pose_2d(angle=0._R8P)
q1 = pose_2d(angle=0._R8P, translation=[0.02_R8P, 0._R8P])            ! a pure translation, 0.8 of a cell
call swept_cells_2d(section, q0, q1, x, x, sw)
print '(A,L1,A,I0,A,I0)', 'section swept: exact ', sw%exact, ', cells with swept area ', count(sw%swept /= 0._R8P), &
                         ', cleared ', count(sw%cleared)
print '(A,L1,A,I0)', '  sum of swept = 0 to 1e-14: ', abs(sum(sw%swept)) < 1.e-14_R8P, &
                    ', cells with residual > 1e-14: ', count(abs(sw%residual) > 1.e-14_R8P)
print '(A,ES9.2,A,2I3,A)', '  largest residual ', maxval(abs(sw%residual)), ' in cell', maxloc(abs(sw%residual)), &
                           ': the knot of the section (chapter 7)'

• swept_cells_2d(curve, pose0, pose1, x, y, sw) (and swept_cells_3d(surface, pose0, pose1, x, y, z, sw)) are free routines with the swept geometry of a step on a fixed grid. A rigid_pose_2d_t comes from pose_2d(angle=, translation=).
  • sw%swept is the fluid area change of each cell, (α₁ − α₀)|cell|, from the cut cells at the two poses.
  • sw%face_u, sw%face_v are the fluid lengths of the faces averaged over the step; sw%wall is the averaged projected wall.
  • sw%cleared and sw%covered flag cells with no fluid at one pose and fluid at the other.
  • sw%residual is the volume balance of the rigid velocity field in each cell: round-off when the face averages are exact and the wall does not cross itself.
• A pure translation takes the exact path (sw%exact): the face lengths are integrated piecewise between the instants when a vertex crosses a grid line or a node crosses a segment. A rotation is integrated by Gauss-Legendre at nq nodes (default 4), with sw%error_estimate.

Mass and inertia

A body moved by the fluid, not along prescribed poses, needs its mass properties to integrate Newton–Euler.

call surface%mass_properties(mass=mass, status=status)
print '(A,L1)', 'the placed dragon is a closed body: ', status == STATUS_OK               ! chapter 2: it is not
call surface%alpha_wrap(alpha=0.03_R8P, offset=0.003_R8P, wrapped=wrap, status=status)     ! a closed surrogate
call wrap%mass_properties(mass=mass, center=centre, inertia=inertia, principal_moments=moments, &
                          density=1000._R8P, status=status)                      ! a body of density 1000
if (status /= STATUS_OK) error stop 'the wrap is not closed'
print '(A,ES11.4,A,3F8.4)', 'mass ', mass, ', centre of mass', centre
print '(A,3ES11.4)', 'principal moments', moments
call wrap%move(p1)
call wrap%mass_properties(inertia=inertia1, density=1000._R8P)
print '(A,ES9.2)', 'moved wrap against inertia_to_world: ', maxval(abs(inertia1 - inertia_to_world(inertia, p1))) / &
                                                           maxval(abs(inertia))

• mass_properties(mass=, center=, inertia=, principal_moments=, principal_axes=, density=, about=, status=) gives the mass, the centre of mass, the inertia tensor about the centre (or about about) and its principal moments and axes. The integrals are exact; see mass properties.
• It needs a closed body, and the repaired dragon is not one (chapter 2). Its alpha wrap is closed by construction, a surrogate 0.003 outside the surface, so the inertia is the wrap's.
• inertia_to_world(inertia, pose) is the tensor at a pose, R I Rᵀ: the moved wrap gives the same tensor to round-off.

The animation integrates Euler's equations with the inertia tensor mass_properties returns for a T-handle: a spin about the intermediate principal axis is unstable, and the handle flips over periodically (the tennis-racket effect).

Distance between bodies

Two bodies in the same flow need the gap between them: for contact, for a lubrication correction when it falls to a few cells, for a collision.

call box%load_from_file(file_name='cube-01.stl', guess_format=.true.)  ! a unit cube beside the dragon
call box%translate(x=0.7_R8P, y=-0.5_R8P, z=-0.5_R8P)
call box%analyze
call surface%distance_to(box, distance=gap, facet_self=fa, facet_other=fb, intersecting=hit, pose_self=p1)
print '(A,F8.5,A,I0,A,I0,A,L1)', 'gap to the cube ', gap, ' (facets ', fa, ', ', fb, '), touching ', hit
call moved%distance_to(box, distance=gap1)                              ! the moved copy of the queries
print '(A,L1)', 'posed = moved, bit for bit: ', gap == gap1

• distance_to(other, distance=, point_self=, point_other=, facet_self=, facet_other=, intersecting=) is the least distance over all pairs of facets of the two bodies, with the pair that realizes it; 0 and intersecting when they touch or cross, decided exactly. The two trees are descended together, best first, and the result is bit-identical to the loop over all pairs.
• pose_self= and pose_other= take the bodies where move would put them, with the trees as built: the result equals the one of the moved copy, bit for bit. proximity_pairs(cutoff=) lists the pairs of bodies closer than a cutoff. See body–body distance.

Running it

$ tut11_moving_bodies
posed = moved to 1e-14: T
cells: fresh 96, dead 98, second query 992 of 13824
section swept: exact T, cells with swept area 131, cleared 16
  sum of swept = 0 to 1e-14: T, cells with residual > 1e-14: 2
  largest residual  5.30E-06 in cell 23 22: the knot of the section (chapter 7)
the placed dragon is a closed body: F
mass  6.0274E+01, centre of mass -0.0203 -0.0306 -0.1039
principal moments 2.0049E+00 3.4901E+00 4.9239E+00
moved wrap against inertia_to_world:  3.63E-16
gap to the cube  0.50539 (facets 1077, 48), touching F
posed = moved, bit for bit: T

• The body-frame distances equal those of the moved and refitted copy to 1e-14.
• The gap to the cube with pose_self= equals the gap of the moved copy bit for bit.
• Of the 13824 cell centres, 96 become fluid and 98 become solid; only 992 needed the second query.
• The section translated by 0.02 sweeps 131 cells, 16 of which go from solid to fluid. The swept areas sum to zero, as a rigid translation keeps the area.
• The residual is round-off in every cell but two, and 5.3e-6 in cell (23, 22). That is the knot of chapter 7: the section crosses itself once, at a fold. The swept areas come from the cut cells (Green's theorem, a region counted by its winding number) and the face averages from line crossings (parity), which count the knot differently; the residual is where they disagree. With the knot cut out of the section the largest residual over the same kind of shifts is 7.8e-18. Pass check_input=.true. to have swept_cells_2d check the curve first (sw%self_intersecting, a warning on stderr), or check it with check_validity(check_self_intersections=.true.).

What you learned

pose= queries a moving body without touching its tree, move refits a copy for grid-aligned queries, state_changes finds fresh and dead cells, swept_cells_2d/swept_cells_3d give the swept geometry of a step, and distance_to the gap between two bodies. Reference: Moving bodies, Body–body distance.

Next: 12. On the GPU.