diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs
new file mode 100644
index 0000000..2d48bb7
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/body.rs
@@ -0,0 +1,159 @@
+//! An embedded body: its signed distance `φ(x, y, z, t)` (φ ≥ 0 = fluid),
+//! the velocity of its surface, its unit normal by central differences,
+//! and a surface sampler (quadrature points with normals and areas) for the
+//! surface-stress load route.
+
+type SdfFn = Box<dyn Fn(f64, f64, f64, f64) -> f64 + Send + Sync>;
+type VelFn = Box<dyn Fn(f64, f64, f64, f64) -> (f64, f64, f64) + Send + Sync>;
+type SamplerFn = Box<dyn Fn(f64) -> Vec<SurfaceSample> + Send + Sync>;
+
+/// One surface quadrature point: position, unit normal out of the solid,
+/// the area it represents.
+#[derive(Debug, Clone, Copy)]
+pub struct SurfaceSample {
+    pub x: f64,
+    pub y: f64,
+    pub z: f64,
+    pub nx: f64,
+    pub ny: f64,
+    pub nz: f64,
+    pub area: f64,
+}
+
+pub struct Body {
+    phi: SdfFn,
+    velocity: Option<VelFn>,
+    sampler: Option<SamplerFn>,
+}
+
+impl Body {
+    pub fn from_sdf<F>(phi: F) -> Self
+    where
+        F: Fn(f64, f64, f64, f64) -> f64 + Send + Sync + 'static,
+    {
+        Self {
+            phi: Box::new(phi),
+            velocity: None,
+            sampler: None,
+        }
+    }
+
+    /// A sphere of radius `r` whose centre is `centre(t)`; its samples (at
+    /// `t = 0`) are a Fibonacci lattice of the requested spacing.
+    pub fn sphere<C>(centre: C, r: f64) -> Self
+    where
+        C: Fn(f64) -> (f64, f64, f64) + Send + Sync + 'static,
+    {
+        let centre = std::sync::Arc::new(centre);
+        let c1 = centre.clone();
+        let mut body = Self::from_sdf(move |x, y, z, t| {
+            let (cx, cy, cz) = c1(t);
+            ((x - cx).powi(2) + (y - cy).powi(2) + (z - cz).powi(2)).sqrt() - r
+        });
+        body.sampler = Some(Box::new(move |ds| {
+            let (cx, cy, cz) = centre(0.0);
+            let n = ((4.0 * std::f64::consts::PI * r * r / (ds * ds)).ceil() as usize).max(32);
+            let golden = std::f64::consts::PI * (3.0 - 5.0_f64.sqrt());
+            (0..n)
+                .map(|k| {
+                    let zz = 1.0 - 2.0 * (k as f64 + 0.5) / n as f64;
+                    let rr = (1.0 - zz * zz).sqrt();
+                    let theta = golden * k as f64;
+                    let (nx, ny, nz) = (rr * theta.cos(), rr * theta.sin(), zz);
+                    SurfaceSample {
+                        x: cx + r * nx,
+                        y: cy + r * ny,
+                        z: cz + r * nz,
+                        nx,
+                        ny,
+                        nz,
+                        area: 4.0 * std::f64::consts::PI * r * r / n as f64,
+                    }
+                })
+                .collect()
+        }));
+        body
+    }
+
+    /// A cylinder along z of radius `r` at `(cx, cy)` (no samples: use
+    /// [`Self::extruded`] for loads).
+    pub fn cylinder_z(cx: f64, cy: f64, r: f64) -> Self {
+        Self::from_sdf(move |x, y, _z, _t| ((x - cx).powi(2) + (y - cy).powi(2)).sqrt() - r)
+    }
+
+    /// A 2D body extruded along z over `[0, lz]`; the samples are the 2D
+    /// samples at z levels of the requested spacing (the end caps carry no
+    /// samples — they lie on the z walls or the periodic seam).
+    pub fn extruded(body2: crate::solvers::incompressible::EmbeddedBody, lz: f64) -> Self {
+        let body2 = std::sync::Arc::new(body2);
+        let (b1, b2, b3) = (body2.clone(), body2.clone(), body2);
+        let mut body = Self::from_sdf(move |x, y, _z, t| b1.phi(x, y, t));
+        body.velocity = Some(Box::new(move |x, y, _z, t| {
+            let (u, v) = b2.surface_velocity(x, y, t);
+            (u, v, 0.0)
+        }));
+        body.sampler = Some(Box::new(move |ds| {
+            let nz = ((lz / ds).ceil() as usize).max(1);
+            let dzs = lz / nz as f64;
+            let mut out = Vec::new();
+            for s in b3.surface_samples(ds) {
+                for k in 0..nz {
+                    out.push(SurfaceSample {
+                        x: s.x,
+                        y: s.y,
+                        z: (k as f64 + 0.5) * dzs,
+                        nx: s.nx,
+                        ny: s.ny,
+                        nz: 0.0,
+                        area: s.ds * dzs,
+                    });
+                }
+            }
+            out
+        }));
+        body
+    }
+
+    #[must_use]
+    pub fn with_surface_velocity<F>(mut self, f: F) -> Self
+    where
+        F: Fn(f64, f64, f64, f64) -> (f64, f64, f64) + Send + Sync + 'static,
+    {
+        self.velocity = Some(Box::new(f));
+        self
+    }
+
+    #[inline]
+    #[must_use]
+    pub fn phi(&self, x: f64, y: f64, z: f64, t: f64) -> f64 {
+        (self.phi)(x, y, z, t)
+    }
+
+    #[inline]
+    #[must_use]
+    pub fn surface_velocity(&self, x: f64, y: f64, z: f64, t: f64) -> (f64, f64, f64) {
+        self.velocity
+            .as_ref()
+            .map_or((0.0, 0.0, 0.0), |f| f(x, y, z, t))
+    }
+
+    /// Unit normal out of the solid (the SDF gradient by central differences).
+    #[must_use]
+    pub fn normal(&self, x: f64, y: f64, z: f64, t: f64, eps: f64) -> (f64, f64, f64) {
+        let gx = (self.phi(x + eps, y, z, t) - self.phi(x - eps, y, z, t)) / (2.0 * eps);
+        let gy = (self.phi(x, y + eps, z, t) - self.phi(x, y - eps, z, t)) / (2.0 * eps);
+        let gz = (self.phi(x, y, z + eps, t) - self.phi(x, y, z - eps, t)) / (2.0 * eps);
+        let norm = (gx * gx + gy * gy + gz * gz).sqrt();
+        if norm > 0.0 {
+            (gx / norm, gy / norm, gz / norm)
+        } else {
+            (1.0, 0.0, 0.0)
+        }
+    }
+
+    /// Surface samples at roughly spacing `ds`; empty for a bare SDF body.
+    #[must_use]
+    pub fn surface_samples(&self, ds: f64) -> Vec<SurfaceSample> {
+        self.sampler.as_ref().map_or_else(Vec::new, |s| s(ds))
+    }
+}
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs
new file mode 100644
index 0000000..0a65ae0
--- /dev/null
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/cut.rs
@@ -0,0 +1,261 @@
+//! The cut geometry of an embedded body on the grid: φ at the cell
+//! corners; inside every cell the interface is the LINEAR interpolant on a
+//! fixed Kuhn split into six tetrahedra (each face into two triangles along
+//! the same diagonal from both sides), so face apertures and cell volumes
+//! are exact for the interpolant, continuous in the corner values, and
+//! consistent across shared faces — the wall polygon's vector area by
+//! closure (`A_w n_w = −Σ_f A_f n_f`) telescopes exactly over a closed body.
+
+use super::Grid;
+use super::body::Body;
+
+/// The cut data of one instant.
+#[derive(Debug, Clone)]
+pub struct CutGeometry {
+    pub grid: Grid,
+    /// φ at the corners, `(nx + 1) × (ny + 1) × (nz + 1)`.
+    pub phi: Vec<f64>,
+    /// Fluid area fraction of every u / v / w face (the staggered layouts).
+    pub a_u: Vec<f64>,
+    pub a_v: Vec<f64>,
+    pub a_w: Vec<f64>,
+    /// Fluid volume fraction of every cell.
+    pub vol: Vec<f64>,
+    /// The wall polygon's vector area per cell (outward from the fluid),
+    /// `−Σ_f A_f n_f · face area`, in area units.
+    pub wall: Vec<[f64; 3]>,
+    /// φ at the face centres (the wall distance of a near-wall face).
+    pub d_u: Vec<f64>,
+    pub d_v: Vec<f64>,
+    pub d_w: Vec<f64>,
+}
+
+impl CutGeometry {
+    #[inline]
+    fn node(g: Grid, k: usize, j: usize, i: usize) -> usize {
+        (k * (g.ny + 1) + j) * (g.nx + 1) + i
+    }
+
+    /// Build the cut data of `body` at time `t`.
+    pub fn build(body: &Body, grid: Grid, t: f64) -> Self {
+        let g = grid;
+        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
+        let mut phi = vec![0.0; (nx + 1) * (ny + 1) * (nz + 1)];
+        for k in 0..=nz {
+            for j in 0..=ny {
+                for i in 0..=nx {
+                    phi[Self::node(g, k, j, i)] =
+                        body.phi(i as f64 * dx, j as f64 * dy, k as f64 * dz, t);
+                }
+            }
+        }
+        let corner = |k: usize, j: usize, i: usize| phi[Self::node(g, k, j, i)];
+        // Face apertures: a face's two triangles along the diagonal from its
+        // (0, 0) to its (1, 1) corner in the face's own (a, b) order — the
+        // Kuhn split's diagonals: for an x-face (y, z), a y-face (x, z), a
+        // z-face (x, y); the same triangles seen from either cell.
+        let tri_area_fraction = |p0: f64, p1: f64, p2: f64| -> f64 {
+            let v = [p0, p1, p2];
+            let pos = v.iter().filter(|&&q| q >= 0.0).count();
+            match pos {
+                0 => 0.0,
+                3 => 1.0,
+                1 => {
+                    let a = v.iter().position(|&q| q >= 0.0).unwrap();
+                    let (b, c) = ((a + 1) % 3, (a + 2) % 3);
+                    (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]))
+                }
+                _ => {
+                    let a = v.iter().position(|&q| q < 0.0).unwrap();
+                    let (b, c) = ((a + 1) % 3, (a + 2) % 3);
+                    1.0 - (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]))
+                }
+            }
+        };
+        // Quad corners in (a, b) order: q00, q10, q01, q11; triangles
+        // (q00, q10, q11) and (q00, q11, q01).
+        let quad_fraction = |q00: f64, q10: f64, q01: f64, q11: f64| -> f64 {
+            0.5 * (tri_area_fraction(q00, q10, q11) + tri_area_fraction(q00, q11, q01))
+        };
+        let mut a_u = vec![0.0; (nx + 1) * ny * nz];
+        let mut a_v = vec![0.0; nx * (ny + 1) * nz];
+        let mut a_w = vec![0.0; nx * ny * (nz + 1)];
+        let mut d_u = vec![0.0; (nx + 1) * ny * nz];
+        let mut d_v = vec![0.0; nx * (ny + 1) * nz];
+        let mut d_w = vec![0.0; nx * ny * (nz + 1)];
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 0..=nx {
+                    // x-face at i: corners (j, k), (j+1, k), (j, k+1), (j+1, k+1)
+                    let (q00, q10, q01, q11) = (
+                        corner(k, j, i),
+                        corner(k, j + 1, i),
+                        corner(k + 1, j, i),
+                        corner(k + 1, j + 1, i),
+                    );
+                    a_u[g.uface(k, j, i)] = quad_fraction(q00, q10, q01, q11);
+                    d_u[g.uface(k, j, i)] = 0.25 * (q00 + q10 + q01 + q11);
+                }
+            }
+        }
+        for k in 0..nz {
+            for j in 0..=ny {
+                for i in 0..nx {
+                    // y-face at j: corners (i, k), (i+1, k), (i, k+1), (i+1, k+1)
+                    let (q00, q10, q01, q11) = (
+                        corner(k, j, i),
+                        corner(k, j, i + 1),
+                        corner(k + 1, j, i),
+                        corner(k + 1, j, i + 1),
+                    );
+                    a_v[g.vface(k, j, i)] = quad_fraction(q00, q10, q01, q11);
+                    d_v[g.vface(k, j, i)] = 0.25 * (q00 + q10 + q01 + q11);
+                }
+            }
+        }
+        for k in 0..=nz {
+            for j in 0..ny {
+                for i in 0..nx {
+                    // z-face at k: corners (i, j), (i+1, j), (i, j+1), (i+1, j+1)
+                    let (q00, q10, q01, q11) = (
+                        corner(k, j, i),
+                        corner(k, j, i + 1),
+                        corner(k, j + 1, i),
+                        corner(k, j + 1, i + 1),
+                    );
+                    a_w[g.wface(k, j, i)] = quad_fraction(q00, q10, q01, q11);
+                    d_w[g.wface(k, j, i)] = 0.25 * (q00 + q10 + q01 + q11);
+                }
+            }
+        }
+        // Cell volumes by the Kuhn split: the six tetrahedra around the
+        // diagonal (0,0,0)–(1,1,1) in unit-cube coordinates.
+        let mut vol = vec![0.0; g.cells()];
+        let mut wall = vec![[0.0; 3]; g.cells()];
+        const KUHN: [[[usize; 3]; 4]; 6] = [
+            [[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]],
+            [[0, 0, 0], [1, 0, 0], [1, 0, 1], [1, 1, 1]],
+            [[0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 1, 1]],
+            [[0, 0, 0], [0, 1, 0], [0, 1, 1], [1, 1, 1]],
+            [[0, 0, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1]],
+            [[0, 0, 0], [0, 0, 1], [0, 1, 1], [1, 1, 1]],
+        ];
+        for k in 0..nz {
+            for j in 0..ny {
+                for i in 0..nx {
+                    let mut fluid = 0.0;
+                    for tet in &KUHN {
+                        let pts: Vec<[f64; 3]> = tet
+                            .iter()
+                            .map(|c| [c[0] as f64, c[1] as f64, c[2] as f64])
+                            .collect();
+                        let vals: Vec<f64> = tet
+                            .iter()
+                            .map(|c| corner(k + c[2], j + c[1], i + c[0]))
+                            .collect();
+                        fluid += tet_fluid_volume(&pts, &vals);
+                    }
+                    // The six tets fill the unit cube (volume 1).
+                    let idx = g.cell(k, j, i);
+                    vol[idx] = fluid;
+                    let ax = dy * dz;
+                    let ay = dx * dz;
+                    let az = dx * dy;
+                    // Outward normals of the cell's faces times their fluid area,
+                    // summed; the wall closes the fluid part of the cell.
+                    let sx = (a_u[g.uface(k, j, i + 1)] - a_u[g.uface(k, j, i)]) * ax;
+                    let sy = (a_v[g.vface(k, j + 1, i)] - a_v[g.vface(k, j, i)]) * ay;
+                    let sz = (a_w[g.wface(k + 1, j, i)] - a_w[g.wface(k, j, i)]) * az;
+                    wall[idx] = [-sx, -sy, -sz];
+                }
+            }
+        }
+        Self {
+            grid,
+            phi,
+            a_u,
+            a_v,
+            a_w,
+            vol,
+            wall,
+            d_u,
+            d_v,
+            d_w,
+        }
+    }
+
+    /// Total fluid volume.
+    #[must_use]
+    pub fn fluid_volume(&self) -> f64 {
+        let g = self.grid;
+        self.vol.iter().sum::<f64>() * g.dx * g.dy * g.dz
+    }
+
+    /// Σ over cells of |wall vector area| (the wall's area up to the
+    /// non-planarity of the per-cell polygon) and the closure vector Σ wall.
+    #[must_use]
+    pub fn wall_area_and_closure(&self) -> (f64, [f64; 3]) {
+        let mut area = 0.0;
+        let mut sum = [0.0; 3];
+        for w in &self.wall {
+            area += (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
+            sum[0] += w[0];
+            sum[1] += w[1];
+            sum[2] += w[2];
+        }
+        (area, sum)
+    }
+}
+
+fn det3(a: [f64; 3], b: [f64; 3], c: [f64; 3]) -> f64 {
+    a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0])
+        + a[2] * (b[0] * c[1] - b[1] * c[0])
+}
+
+fn tet_volume(p: [[f64; 3]; 4]) -> f64 {
+    let e = |a: [f64; 3], b: [f64; 3]| [b[0] - a[0], b[1] - a[1], b[2] - a[2]];
+    det3(e(p[0], p[1]), e(p[0], p[2]), e(p[0], p[3])).abs() / 6.0
+}
+
+fn lerp(a: [f64; 3], b: [f64; 3], t: f64) -> [f64; 3] {
+    [
+        a[0] + t * (b[0] - a[0]),
+        a[1] + t * (b[1] - a[1]),
+        a[2] + t * (b[2] - a[2]),
+    ]
+}
+
+/// The volume of `{φ ≥ 0}` in a tetrahedron with the linear interpolant of
+/// the corner values `v` (φ = 0 at a corner counts as fluid). Continuous in
+/// `v`: every case's cut points move continuously and the cases agree on
+/// their boundaries.
+fn tet_fluid_volume(pts: &[[f64; 3]], v: &[f64]) -> f64 {
+    let total = tet_volume([pts[0], pts[1], pts[2], pts[3]]);
+    let pos: Vec<usize> = (0..4).filter(|&q| v[q] >= 0.0).collect();
+    let neg: Vec<usize> = (0..4).filter(|&q| v[q] < 0.0).collect();
+    let cut = |a: usize, b: usize| lerp(pts[a], pts[b], v[a] / (v[a] - v[b]));
+    match pos.len() {
+        0 => 0.0,
+        4 => total,
+        1 => {
+            let a = pos[0];
+            let (b, c, d) = (neg[0], neg[1], neg[2]);
+            tet_volume([pts[a], cut(a, b), cut(a, c), cut(a, d)])
+        }
+        3 => {
+            let a = neg[0];
+            let (b, c, d) = (pos[0], pos[1], pos[2]);
+            total - tet_volume([pts[a], cut(a, b), cut(a, c), cut(a, d)])
+        }
+        _ => {
+            // Two fluid corners A, B; the fluid wedge {A, B, P_AC, P_AD, P_BC, P_BD}
+            // as three tetrahedra.
+            let (a, b) = (pos[0], pos[1]);
+            let (c, d) = (neg[0], neg[1]);
+            let (pac, pad, pbc, pbd) = (cut(a, c), cut(a, d), cut(b, c), cut(b, d));
+            tet_volume([pts[a], pts[b], pbc, pbd])
+                + tet_volume([pts[a], pac, pbc, pbd])
+                + tet_volume([pts[a], pac, pad, pbd])
+        }
+    }
+}
diff --git a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs
index 2a6036e..a0f9f97 100644
--- a/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs
+++ b/crates/specialized/rtx-cfd/src/solvers/incompressible/embedded3/mod.rs
@@ -6,11 +6,15 @@
 //! Layout: cells `(k, j, i)` row-major, `cell = (k·ny + j)·nx + i`;
 //! u faces on `(nx + 1)·ny·nz`, v on `nx·(ny + 1)·nz`, w on `nx·ny·(nz + 1)`.
 
+pub mod body;
+pub mod cut;
 pub mod field;
 pub mod grid;
 pub mod poisson;
 pub mod step;
 
+pub use body::{Body, SurfaceSample};
+pub use cut::CutGeometry;
 pub use field::Field;
 pub use grid::Grid;
 pub use step::{Boundaries, Fluid, Parameters, Side, Solver, StepResult};
diff --git a/crates/specialized/rtx-cfd/tests/embedded3_cut_geometry.rs b/crates/specialized/rtx-cfd/tests/embedded3_cut_geometry.rs
new file mode 100644
index 0000000..ebf63d7
--- /dev/null
+++ b/crates/specialized/rtx-cfd/tests/embedded3_cut_geometry.rs
@@ -0,0 +1,172 @@
+//! embedded3 gate 8: the cut geometry. A sphere and a z-cylinder on
+//! n = 16 / 32 / 64: the fluid volume and the wall area converge at second
+//! order; the wall closure Σ_c A_w n_w vanishes to rounding for a body
+//! inside the box; a continuity sweep of the sphere's centre across one
+//! cell (200 positions): bounded difference quotient of every aperture and
+//! volume, no jump > 1e-3 between neighbouring positions; and the
+//! translating sphere's discrete volume change per step, MEASURED (the
+//! plan's "to rounding" clause is checked, not assumed).
+
+use rtx_cfd::solvers::incompressible::embedded3::{Body, CutGeometry, Grid};
+use std::f64::consts::PI;
+
+fn cube(n: usize) -> Grid {
+    let h = 1.0 / n as f64;
+    Grid {
+        nx: n,
+        ny: n,
+        nz: n,
+        dx: h,
+        dy: h,
+        dz: h,
+    }
+}
+
+#[test]
+fn volume_and_area_converge_at_second_order_and_the_wall_closes() {
+    let r: f64 = 0.3;
+    let sphere_v = 4.0 / 3.0 * PI * r.powi(3);
+    let sphere_a = 4.0 * PI * r * r;
+    let cyl_v = PI * r * r * 1.0;
+    let cyl_a = 2.0 * PI * r * 1.0;
+    for (name, exact_v, exact_a, body) in [
+        (
+            "sphere",
+            1.0 - sphere_v,
+            sphere_a,
+            Body::sphere(|_t| (0.5, 0.5, 0.5), r),
+        ),
+        (
+            "z-cylinder",
+            1.0 - cyl_v,
+            cyl_a,
+            Body::cylinder_z(0.5, 0.5, r),
+        ),
+    ] {
+        let mut ev = Vec::new();
+        let mut ea = Vec::new();
+        for n in [16usize, 32, 64] {
+            let g = CutGeometry::build(&body, cube(n), 0.0);
+            let v = g.fluid_volume();
+            let (a, closure) = g.wall_area_and_closure();
+            let closure_norm =
+                (closure[0].powi(2) + closure[1].powi(2) + closure[2].powi(2)).sqrt();
+            ev.push((v - exact_v).abs());
+            ea.push((a - exact_a).abs());
+            println!(
+                "  {name} n {n}: fluid volume {v:.8} (exact {exact_v:.8}, err {:.2e}); wall area {a:.6} (exact {exact_a:.6}, err {:.2e}); closure |Σ A_w n_w| {closure_norm:.2e}",
+                ev.last().unwrap(),
+                ea.last().unwrap()
+            );
+            // The z-cylinder touches the z walls: its closure includes the
+            // end caps' missing area only through the cell walls, so the
+            // closure holds for the sphere; for the cylinder the z component
+            // is the two caps (equal and opposite) and x, y close.
+            if name == "sphere" {
+                assert!(
+                    closure_norm < 1e-12,
+                    "{name} n {n}: closure {closure_norm:.3e}"
+                );
+            } else {
+                assert!(
+                    closure[0].abs() < 1e-12 && closure[1].abs() < 1e-12,
+                    "{name} n {n}: closure x/y {closure:?}"
+                );
+            }
+        }
+        for (label, e) in [("volume", &ev), ("area", &ea)] {
+            for w in e.windows(2) {
+                let order = (w[0] / w[1]).log2();
+                println!("  {name} {label} order {order:.2}");
+                assert!(order > 1.5, "{name} {label}: order {order:.2} below second");
+            }
+        }
+    }
+}
+
+/// Continuity in the body position. A face aperture where the interface is
+/// tangent to the face changes at a rate of order `r / h` per cell width of
+/// shift (the cap's area grows linearly in the shift, on a face of area
+/// h²), so an O(1) Lipschitz bound is the wrong premise; the discriminating
+/// test is that the largest change between neighbouring positions falls in
+/// proportion when the sweep is refined tenfold — a discontinuity would not.
+fn sweep(positions: usize) -> f64 {
+    let n = 16;
+    let g = cube(n);
+    let h = g.dx;
+    let r: f64 = 0.3;
+    let mut prev: Option<CutGeometry> = None;
+    let mut worst_jump = 0.0_f64;
+    for s in 0..=positions {
+        let shift = h * s as f64 / positions as f64;
+        let body = Body::sphere(
+            move |_t| (0.5 + shift, 0.5 + 0.37 * shift, 0.5 + 0.11 * shift),
+            r,
+        );
+        let cut = CutGeometry::build(&body, g, 0.0);
+        if let Some(p) = &prev {
+            let jump = |a: &[f64], b: &[f64]| {
+                a.iter()
+                    .zip(b)
+                    .fold(0.0_f64, |m, (x, y)| m.max((x - y).abs()))
+            };
+            let j = jump(&cut.a_u, &p.a_u)
+                .max(jump(&cut.a_v, &p.a_v))
+                .max(jump(&cut.a_w, &p.a_w))
+                .max(jump(&cut.vol, &p.vol));
+            worst_jump = worst_jump.max(j);
+        }
+        prev = Some(cut);
+    }
+    worst_jump
+}
+
+#[test]
+fn apertures_and_volumes_are_continuous_in_the_body_position() {
+    let j200 = sweep(200);
+    let j2000 = sweep(2000);
+    let ratio = j2000 / j200;
+    println!(
+        "  sphere over one cell: largest neighbour change {j200:.3e} at 200 positions, {j2000:.3e} at 2000 (ratio {ratio:.3}; 0.1 = Lipschitz, 1 = a jump)"
+    );
+    assert!(
+        ratio < 0.2,
+        "the largest change does not fall with the sweep resolution (ratio {ratio:.3}): a discontinuity"
+    );
+    assert!(
+        ratio > 0.05,
+        "ratio {ratio:.3} below the Lipschitz expectation — check the sweep"
+    );
+}
+
+#[test]
+fn the_translating_sphere_volume_change_is_measured() {
+    let n = 32;
+    let g = cube(n);
+    let r: f64 = 0.3;
+    let dt = 1e-3;
+    let speed = 1.0; // one cell width in n·dt... 0.03125 m per 31 steps
+    let mut worst_rel = 0.0_f64;
+    let mut prev = CutGeometry::build(
+        &Body::sphere(move |t| (0.5 + speed * t, 0.5, 0.5), r),
+        g,
+        0.0,
+    );
+    let swept_per_step = PI * r * r * speed * dt; // the sphere's cross-section swept
+    for s in 1..=31 {
+        let t = s as f64 * dt;
+        let cut = CutGeometry::build(&Body::sphere(move |t| (0.5 + speed * t, 0.5, 0.5), r), g, t);
+        let dv = cut.fluid_volume() - prev.fluid_volume();
+        worst_rel = worst_rel.max(dv.abs() / swept_per_step);
+        prev = cut;
+    }
+    let body_v = 4.0 / 3.0 * PI * r.powi(3);
+    println!(
+        "  translating sphere n {n}: largest |Σ ΔV_c| per step = {worst_rel:.3e} of the swept cross-section volume per step ({:.2e} of the body volume)",
+        worst_rel * swept_per_step / body_v
+    );
+    // Recorded; the compatibility of the moving-body projection is decided
+    // on this number (the plan's clause "to rounding" is not true for the
+    // piecewise-linear interpolant — the cut cells' volume error moves with
+    // the body).
+}