Moving bodies

What it does

A rigidly moving body needs its geometry at every new pose. Rebuilding the tree every time step is the cost to avoid. FOSSIL offers two mechanisms, and the kind of query decides which one to use:

Mechanism	Use it for	Cost	Agreement with a mesh moved and rebuilt
Body-frame queries: pass pose=	distance, IB queries, inside tests (CPU and device)	one 3 × 3 product per point	distances ≤ 1e-16 (L + |t|), signs equal away from the surface
Moved copy: move(pose) with a tree refit	grid-aligned queries: box tagging, line crossings, cut cells	O(N) (a refit, not a rebuild)	bit for bit (feet at edge and vertex ties to 4 ulps)

state_changes flags the points (cell centres) that change state between two poses: fresh (solid → fluid) and dead (fluid → solid).

swept_cells_2d / swept_cells_3d give the swept geometry a conservative moving cut-cell update needs over a step (see Swept geometry).

Usage

program moving
use fossil, only : surface_stl_object, rigid_pose_t, pose_from_axis_angle, pose_identity
use penf,   only : R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: body, moved
type(rigid_pose_t)       :: p0, p1
type(vector_R8P)         :: centres(2)
real(R8P)                :: d(2)
logical                  :: fresh(2), dead(2)

call body%load_from_file(file_name='body.stl', guess_format=.true.)
p0 = pose_identity()
p1 = pose_from_axis_angle(axis=[0._R8P, 0._R8P, 1._R8P], angle=0.1_R8P, translation=[0.05_R8P, 0._R8P, 0._R8P])
centres = [vector_R8P(0.5_R8P, 0.5_R8P, 0.5_R8P), vector_R8P(1.2_R8P, 0.3_R8P, 0.4_R8P)]

! body-frame query: the points are in the world frame, the tree stays as built
call body%distance_many(points=centres, distances=d, is_signed=.true., pose=p1)
! cells that change state between the two poses
call body%state_changes(pose0=p0, pose1=p1, points=centres, fresh=fresh, dead=dead)
! a moved copy for grid-aligned queries (box tagging, line crossings, cut cells)
moved = body
call moved%move(p1)
endprogram moving

Poses

• rigid_pose_t holds r(3, 3) and t(3), with x_world = R x_body + t.
  • Constructors: pose_identity, pose_from_axis_angle(axis, angle, translation), pose_from_quaternion(q, translation).
  • Algebra: pose_compose(a, b) (apply b, then a) and pose_inverse.
  • Maps: to_body, to_world, rotate, which are device-callable.
• rigid_pose_2d_t (pose_2d(angle, translation)) is the in-plane pose of a curve.
• pose_interpolate(p0, p1, s) (2-D and 3-D) is the screw motion between two poses: the translation and the rotation angle, about the fixed axis of R₁R₀ᵀ, are linear in s ∈ [0, 1]. The end poses are returned exactly.

The pose argument is on:

• distance_many, query_many, is_point_inside and their device twins distance_many_device, query_many_device;
• for surfaces and curves alike.

Closest points, second foot points and normals come back in the world frame.

How it works

• Body-frame queries. The points are mapped by x̃ = Rᵀ (x − t); the results are mapped back. Distances are invariant. The device twins map on the host with the CPU code, so the posed device results equal the posed CPU results bit for bit.
• Move and refit.
  • The vertices move, and everything derived from coordinates is recomputed exactly as analyze does.
  • The SAH BVH keeps its topology: node boxes are recomputed bottom-up as unions of facet boxes, and the payloads are rebuilt. Curves refit their segment tree the same way.
  • Topology, shells and bodies are kept. An octree is rebuilt.
  • The facets equal those of the mesh moved and analyzed from scratch, so the queries agree with a rebuild.
  • Exception: a foot on an edge or vertex shared by tied facets is parametrized differently by each; which facet is reported depends on the tree, so its bits can differ by a few ulps.
• Fresh and dead points.
  • A point at the two body-frame positions b₀ and b₁ can change state only if the segment b₀b₁ reaches the surface.
  • So when |φ(b₀)| > |b₀ − b₁| the second query is skipped. That filter is exact.
  • States come from the signed distance with the chosen sign algorithm.

Swept geometry

A conservative cut-cell update over a step with a moving wall needs, per cell, the old and new fluid volumes. Per grid face it needs the fluid aperture averaged over the step (Falcovitz et al. 1997). Three decisions follow from the fixed grid (issue #41):

• Nothing sweeps through a fixed face.
  • The volume the wall sweeps through a cell is its fluid volume change, by Reynolds transport.
  • So swept = (α₁ − α₀)·|cell| comes exactly from the cut cells at the two poses, in 2-D and 3-D.
  • The geometric conservation law holds by construction.
  • Natarajan et al. (2022) omit this term and lose order.
• The only time integral is the face aperture.
  • face_u / face_v (2-D) and face_x / face_y / face_z (3-D) are the fluid length or area of every grid face, averaged over the step.
  • wall = −Σ_f Ā_f n_f is the time-averaged projected wall area of the cell, which carries the wall pressure and work. It uses the same quantities as the face fluxes.
• Motion. By default the motion is pose_interpolate. A caller that knows the true motion, such as a 6-DOF integrator, can instead pass the poses at the Gauss nodes (gauss_legendre_01).

call swept_cells_2d(curve, pose0, pose1, x, y, sw, nq, poses, exact, status, check_input)
call swept_cells_3d(surface, pose0, pose1, x, y, z, sw, nq, poses, exact, status, check_input)

Integration.

• 2-D pure translation: exact. Both poses share their rotation.
  • Each face is an axis-aligned segment moving against the curve, so its fluid length is piecewise linear in time.
  • The breakpoints are the instants when a wall vertex crosses the face's line, or when an end node of the face crosses a wall segment.
  • Two-point Gauss on every piece is exact.
• 3-D pure translation: exact. The same idea, one grid plane at a time.
  • In the frame of the body, a grid plane moves by −s d with its grid lines. Every vertex of a face's fluid polygon moves linearly in time while the combinatorics hold, so the area is quadratic in time.
  • The breakpoints of a plane are the instants when a mesh vertex crosses it, a mesh edge crosses one of its grid lines, or one of its grid nodes crosses a triangle.
  • Two-point Gauss on every piece is exact. One evaluation (a section of the shifted plane and the 2-D cut cells on the shifted grid) gives every face of the plane.
• exact=.false. forces the Gauss path for a translation (both dimensions); sw%exact reports which path ran and sw%evaluations how many aperture evaluations it took (per face in 2-D, per plane in 3-D).
• 2-D rotation: exact (issue #69, see Exact rotation below). The screw motion with a rotation takes certified breakpoints per face and Gauss–Legendre doubled to round-off on every piece; the residual falls to round-off (1.7e-16 to 6.3e-16 on the cases below, where 4-node Gauss left 1.7e-4 to 1.4e-2). sw%evaluations is the largest number of evaluations per face (92 for the square, 1838 for the 200-vertex NACA 0012).
• 3-D rotation: exact (issue #69). The same scheme one grid plane at a time: a plane's breakpoints are a mesh vertex on it, a mesh edge meeting one of its grid lines, one of its grid nodes on a triangle; one evaluation slices a moved copy and gives every face of the plane. Residual 1.4e-16 to 1.8e-15 on the cases below, where 4-node Gauss left 1.0e-3 to 9.3e-3.
• Otherwise (explicit node poses, or exact=.false.): Gauss–Legendre at nq nodes (default 4), on copies moved to the node poses. With the screw motion, error_estimate is the largest face change between nq and 2 nq nodes.

Exact rotation (design, issue #69). The Gauss path converges algebraically under rotation (rate 2–3, see the validation table), because a face's fluid measure has kinks inside the step. The exact path removes them, as the translation path does, and is built as follows.

• Motion. The screw motion of pose_interpolate has the angle linear in s ∈ [0, 1] about a fixed axis (the plane normal in 2-D) and the translation linear in s: a body point q is at X(s) = R_k(sΔθ) R₀ q + t₀ + s Δt. It is not a rotation about a fixed pivot, so the events below are transcendental in general.

• Events. A face's fluid measure is smooth between the instants when

  • a wall vertex crosses the face's line (2-D) or plane (3-D);
  • an end node of the face crosses a wall segment (2-D), or a grid node of the plane crosses a triangle (3-D);
  • a mesh edge crosses one of the plane's grid lines (3-D).

  Each is a zero of an orientation of moving and fixed points, and with X(s) as above every one has the form

  g(s) = (p₀ + p₁ s) cos φ + (q₀ + q₁ s) sin φ + k₀ + k₁ s, φ = φ₀ + s Δθ,

  with coefficients from the input points (the rotation acts on differences of body points, the translation enters linearly), plus a check that the event lies on the face (the vertex inside the span, the node's foot inside the segment or triangle).

• Certified isolation. |g''(s)| ≤ M = Δθ² max|(p, q)| + 2|Δθ| |(p₁, q₁)|, the maximum over [0, 1]. On an interval of length h whose end values have the same sign and exceed M h²/8 in magnitude, g has no zero (g lies above its chord minus M h²/8). Otherwise the interval is halved; a sign change is refined by bisection to one ulp of s. No pair of zeros is missed. An interval narrower than 1e-10 without a certificate (a tangency, a zero at an end of the step) gives one breakpoint; an interval on which |g| stays at round-off throughout (a zero zone: a mesh face starting on a grid plane) gets breakpoints at its two ends only; breakpoints within 1e-10 of each other merge. A kink that close to a breakpoint changes its piece's Gauss value by ~slope × 1e-20. If the 4096-interval budget of one event runs out, sw%exact is false.

• Integration. On each piece between breakpoints the fluid measure is analytic, so Gauss–Legendre converges exponentially. The node count on a piece doubles from 2 until two successive values agree to round-off of the face measure (pieces shorter than 1e-6 to 64 ε of the measure times 1e-6: slices within 1e-10 of a mesh face carry ~1e-7 relative noise); sw%evaluations reports the largest count (per face in 2-D, per plane in 3-D). A piece that has not converged at 64 nodes marks the result sw%exact = .false., never silently. Each evaluation moves a copy of the body to the node's pose and measures the face exactly as the cut cells do, so apertures and fractions share one convention.

• What changes. The swept volumes are exact already (cut-cell volumes at the two poses), and wall closes the uniform flow by construction, so a conservative update is consistent either way (Falcovitz's requirement). The exact path makes the time-averaged apertures exact under rotation too, and with them the rigid-field residual, which falls to round-off. Explicit node poses (poses=) keep the Gauss path; exact=.false. forces it.

Check. residual is the volume balance of the body's rigid velocity field: swept volume plus the time integral of the field's outflow through the fluid part of every face. It vanishes when a uniform flow moving with the body is reproduced exactly, Falcovitz's design requirement. It is not computed (0) when the caller passes the node poses. The residual also checks the input: the exact paths assume a wall that does not cross itself (see the limitation below).

cleared and covered flag the cells with no fluid at one pose and fluid at the other.

Results

Case	Result
body-frame vs moved and rebuilt: cube (also translated by 108), cone, L-prism, bunny, NACA 0012, annulus	distances ≤ 1e-16 (L + |t|), closest points ≤ 6.5e-15 (L + |t|), signs and inside flags equal
posed device vs posed CPU	every output bit for bit (device path run on the host)
move (refit) vs moved and rebuilt: cube, cone, L-prism, bunny, octree cube, NACA 0012, annulus	geometry, distances, inside flags, box codes, line crossings, cut cells bit for bit
refit time / rebuild time	0.03 (bunny) to 0.4 (tiny meshes)
fresh / dead: cube translated, rotated about an off-centre axis, general poses; square, pitching NACA 0012	0 errors against the analytic change of membership; the band skips 9–60 % of the second queries
swept, 2-D translation: square (every face vs analytic), non-convex heptagon, heptagon at a rotated pose	faces to 8.3e-17; uniform-flow residual ≤ 1.7e-16 in every cell; walls sum to ≤ 1.4e-16
swept, 2-D negative control: Gauss with 4 nodes, no breakpoints	residual 9.0e-3
swept, 2-D rotation, Gauss (exact=.false.): off-centre square, pitching NACA 0012; nq = 2, 4, 8, 16	residual 5.6e-3 → 1.6e-4 (rate ≈ 2.3) and 1.5e-3 → 1.6e-6 (rate ≈ 3.4)
swept, 3-D rotation, exact (issue #69): cube 0.3 rad about (1, 1, 1) + shift (faces starting on grid planes), L-prism 0.5 rad about a tilted axis + shift, 32-face cone pitching 0.4 rad	residual 1.4e-16, 1.8e-15, 5.0e-16 (6³ validation: 4.3e-17, 8.3e-16, 5.0e-16); walls sum ≤ 1.8e-15; Gauss nq = 8 ≥ 4× closer than nq = 2
swept, 2-D rotation, exact (issue #69): square 0.4 rad off-centre + shift, square 1.2 rad + shift, NACA 0012 pitching 0.15 rad	residual 1.7e-16, 1.8e-16, 6.3e-16; walls sum ≤ 1.9e-16; Gauss nq = 16 ≥ 4× closer to the exact faces than nq = 4
swept, 3-D translation: unit cube (every face vs analytic)	faces to 5.6e-17; uniform-flow residual 2.5e-17; walls sum to ≤ 8.3e-17
swept, 3-D translation: cube and L-prism at a pose rotated about (1, 1, 1), 32-face cone	uniform-flow residual 7.6e-17, 3.0e-16, 1.9e-16
swept, 3-D translation without one breakpoint class (vertex / edge-line / node-triangle)	residual up to 5.7e-4 / 1.4e-3 / 1.3e-5: each class is needed
swept, 3-D negative control: unit cube translating, Gauss, nq = 2, 4, 8, 16	faces vs analytic 2.5e-2 → 2.3e-3; residual 1.3e-2 → 7.1e-4
swept, 3-D: cube rotating about (1, 1, 1) while translating; nq = 2, 4, 8	residual 1.1e-2 → 1.6e-3
swept: explicit node poses = screw poses	face apertures bit for bit

Limitations

• Swept geometry: a self-intersecting wall breaks the balance. The swept volumes come from the cut cells (Green's theorem: a region counts by its winding number) and the face apertures from line crossings (parity). Where the wall crosses itself the two count the knot differently, so the cells around the knot get a residual that is not round-off, even on the exact translation path. Measured on a section of the (non-manifold) dragon that crosses itself once at a fold: residual 5.3e-6 in the cell the knot enters; with the knot removed, 7.8e-18 over the same 25 shifts. Pass check_input=.true. to swept_cells_2d / swept_cells_3d: the input is checked first (one self_intersections, or find_self_intersections in 3-D, which costs more), sw%self_intersecting is set and a warning goes to stderr; the result is computed anyway.
• Swept geometry: rotations converge slowly across jumps.
  • A face's fluid length (area) jumps when the face's line (plane) leaves the body. Gauss across a jump is O(1/nq).
  • Rotation in 2-D gives rates of about 2–3; in 3-D about first order (cube: residual 1.1e-2 → 1.6e-3 from nq = 2 to 8).
  • Only translations are exact: under a rotation the polygon vertices no longer move linearly, so there is no finite set of pieces on which a fixed rule is exact.
• Swept geometry: the cost. The 2-D exact path costs O((faces + nodes) × segments). The 3-D exact path costs, per grid plane, a scan of the triangles in the slab the plane sweeps, plus 2 (breakpoints + 1) sections and 2-D cut-cell passes; the breakpoints grow with the mesh vertices and edges in that slab, so a fine mesh and a long step make it expensive (sw%evaluations reports it). The Gauss paths cost one moved copy and one face pass per node, plus as many again for the error estimate.
• Rigid bodies only. Deforming bodies, and refitting under deformation, are out of scope.
• Points near the surface. state_changes decides states with the point-query sign, so a point within round-off of the wall can differ from the cut-cell parity.
• GPU not validated yet. The posed device paths run on the host in tests-gnu; the GPU benches print posed device vs serial mismatches for the manual gate.