//! The apertured cut-cell wall (`WallScheme::CutCell`, item 10 — the
//! "AM-wall"): the cut geometry classifies the grid (a cell is fluid where
//! its fluid volume is positive; an interior face is an unknown where its
//! aperture is positive, prescribed the surface velocity otherwise — no
//! ghost faces), and the wall enters the operators through the apertures:
//! continuity `Σ_f A_f u_f·n_f + U_b·W_c = 0` per cell, the projection
//! coefficient `dt·A_f/δ`, the momentum control volume of an unknown face
//! `V_u = A_f h` with its own faces' apertures averaged from the two
//! adjacent cells (mass fluxes averaged, so the momentum volume conserves
//! mass exactly when the cells do), the wall closing it (`W = −Σ A n`),
//! an implicit wall shear `μ A_w (u_f − U_b)/d_f`, and an inertia floor of
//! 0.1 in the time derivative only. Every prescribed value is the limit of
//! the computed one as the aperture closes (the shear coefficient grows as
//! `1/A_f`), which is what makes the wall smooth in the interface position.

use super::body::Body;
use super::cut::CutGeometry;
use super::exchange::in_load_window;
use super::field::Field;
use super::step::{Boundaries, Side};
use super::wall::{FaceKind, Mask};
use super::Grid;

/// The inertia floor: the momentum volume's fraction in the time
/// derivative is at least this.
pub(super) const INERTIA_FLOOR: f64 = 0.1;
/// The wall-distance floor of a face, in units of the smallest spacing.
pub(super) const DISTANCE_FLOOR: f64 = 0.05;
/// The fine floor (S2-6, `Parameters::distance_floor_fine`).
pub(super) const DISTANCE_FLOOR_FINE: f64 = 0.01;
/// Virtual merging: a cell whose fluid fraction (at either end of the
/// step) stays below this shares its pressure unknown with a neighbour.
pub(super) const MERGE_FRACTION: f64 = 0.1;

/// `RTX_E3_MOVING_PROFILE` (read once).
pub(super) fn profile_on() -> bool {
    static ON: std::sync::OnceLock<bool> = std::sync::OnceLock::new();
    *ON.get_or_init(|| std::env::var("RTX_E3_MOVING_PROFILE").is_ok())
}

/// Lattice addressing of faces and cells with the periodic wrap in z as
/// data: a face of component `c` at `p = [i, j, k]` (its own coordinate is
/// the face index, the others the cell's), a cell at `[i, j, k]`.
#[derive(Debug, Clone, Copy)]
pub(super) struct Lattice {
    pub(super) g: Grid,
    pub(super) periodic_z: bool,
}

impl Lattice {
    fn wrap_z(&self, k: i64, planes: i64) -> Option<usize> {
        if self.periodic_z {
            Some(k.rem_euclid(self.g.nz as i64) as usize)
        } else if (0..planes).contains(&k) {
            Some(k as usize)
        } else {
            None
        }
    }

    /// The index of the face of component `c` at `p`, `None` outside.
    pub(super) fn face(&self, c: usize, p: [i64; 3]) -> Option<usize> {
        let (nx, ny, nz) = (self.g.nx as i64, self.g.ny as i64, self.g.nz as i64);
        let (i, j) = (p[0], p[1]);
        let (ni, nj, nk) = match c {
            0 => (nx + 1, ny, nz),
            1 => (nx, ny + 1, nz),
            _ => (nx, ny, nz + 1),
        };
        if !(0..ni).contains(&i) || !(0..nj).contains(&j) {
            return None;
        }
        let k = self.wrap_z(p[2], nk)?;
        let (i, j) = (i as usize, j as usize);
        Some(match c {
            0 => self.g.uface(k, j, i),
            1 => self.g.vface(k, j, i),
            _ => self.g.wface(k, j, i),
        })
    }

    /// The index of the cell at `p`, `None` outside.
    pub(super) fn cell(&self, p: [i64; 3]) -> Option<usize> {
        let (nx, ny, nz) = (self.g.nx as i64, self.g.ny as i64, self.g.nz as i64);
        if !(0..nx).contains(&p[0]) || !(0..ny).contains(&p[1]) {
            return None;
        }
        let k = self.wrap_z(p[2], nz)?;
        Some(self.g.cell(k, p[1] as usize, p[0] as usize))
    }

    /// The centre of the face of component `c` at `p`.
    pub(super) fn face_position(&self, c: usize, p: [i64; 3]) -> [f64; 3] {
        let h = [self.g.dx, self.g.dy, self.g.dz];
        let mut x = [0.0; 3];
        for d in 0..3 {
            let off = if d == c { 0.0 } else { 0.5 };
            x[d] = (p[d] as f64 + off) * h[d];
        }
        x
    }
}

/// The geometry of an unknown face's momentum control volume.
#[derive(Debug, Clone, Copy)]
pub(super) struct CvGeometry {
    /// The face's own aperture.
    pub(super) alpha: f64,
    /// The control volume's face apertures `[direction][minus, plus]`.
    pub(super) ap: [[f64; 2]; 3],
    /// The wall's vector area closing the control volume (into the body).
    pub(super) wall: [f64; 3],
    /// The wall distance of the face (floored).
    pub(super) distance: f64,
}

impl Mask {
    /// Classify the grid against `body` at `t` by its cut geometry.
    pub fn build_cut(body: &Body, g: Grid, t: f64, b: Boundaries) -> Result<Self, String> {
        Self::build_cut_from(body, g, t, b, None)
    }

    /// As [`Self::build_cut`], re-evaluating φ only within `band` of the
    /// previous geometry moved by at most `motion` (see
    /// [`CutGeometry::build_from`]).
    pub fn build_cut_from(
        body: &Body,
        g: Grid,
        t: f64,
        b: Boundaries,
        prev: Option<(&CutGeometry, f64, f64)>,
    ) -> Result<Self, String> {
        let lap = std::time::Instant::now();
        let cut = CutGeometry::build_from(body, g, t, prev);
        if std::env::var("RTX_E3_MOVING_PROFILE").is_ok() {
            eprintln!(
                "    mask laps: cut geometry (host, within build_mask) {:.0} ms",
                lap.elapsed().as_secs_f64() * 1e3
            );
        }
        Self::from_cut(cut, g, b)
    }

    /// R6-1: the mask of a cut geometry built elsewhere (the device's
    /// `DeviceGeom` mirror); `build_cut_from` is this after `CutGeometry::build_from`.
    pub fn from_cut(cut: CutGeometry, g: Grid, b: Boundaries) -> Result<Self, String> {
        let lap = std::time::Instant::now();
        let (nx, ny, nz) = (g.nx, g.ny, g.nz);
        let periodic = b.z0 == Side::Periodic;
        let allowed = |side: Side| matches!(side, Side::Velocity | Side::Periodic | Side::SlipWall);
        // P1-3 (f): the whole-grid classifications as parallel per-entry
        // maps (the same values; the anchor is the smallest fluid index, the
        // serial loop's first).
        use rayon::prelude::*;
        let nxy = nx * ny;
        let cell_fluid: Vec<bool> = cut.vol.par_iter().map(|&v| v > 0.0).collect();
        let fluid_cells = cell_fluid.par_iter().filter(|&&f| f).count();
        let anchor = (0..g.cells())
            .into_par_iter()
            .find_first(|&idx| cell_fluid[idx]);
        let touching = (0..g.cells()).into_par_iter().find_first(|&idx| {
            if cell_fluid[idx] {
                return false;
            }
            let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx);
            (i == 0 && !allowed(b.x0))
                || (i + 1 == nx && !allowed(b.x1))
                || (j == 0 && !allowed(b.y0))
                || (j + 1 == ny && !allowed(b.y1))
                || (k == 0 && !allowed(b.z0))
                || (k + 1 == nz && !allowed(b.z1))
        });
        if let Some(idx) = touching {
            let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx);
            return Err(format!(
                "embedded body reaches a domain side that is not a Velocity/Periodic side at cell ({k}, {j}, {i})"
            ));
        }
        let Some(anchor) = anchor else {
            return Err("embedded body covers the whole domain".into());
        };
        let kind = |a: f64| {
            if a > 0.0 {
                FaceKind::Fluid
            } else {
                FaceKind::Solid
            }
        };
        // Domain-side faces keep `Fluid` (the serial loops skipped them).
        let u_kind: Vec<FaceKind> = (0..g.n_ufaces())
            .into_par_iter()
            .map(|f| {
                let i = f % (nx + 1);
                if i == 0 || i == nx {
                    FaceKind::Fluid
                } else {
                    kind(cut.a_u[f])
                }
            })
            .collect();
        let v_kind: Vec<FaceKind> = (0..g.n_vfaces())
            .into_par_iter()
            .map(|f| {
                let j = (f / nx) % (ny + 1);
                if j == 0 || j == ny {
                    FaceKind::Fluid
                } else {
                    kind(cut.a_v[f])
                }
            })
            .collect();
        let w_kind: Vec<FaceKind> = (0..g.n_wfaces())
            .into_par_iter()
            .map(|f| {
                let k = f / nxy;
                if !periodic && (k == 0 || k == nz) {
                    FaceKind::Fluid
                } else {
                    kind(cut.a_w[f])
                }
            })
            .collect();
        let mut mask = Self::from_parts(
            g,
            periodic,
            cell_fluid,
            [u_kind, v_kind, w_kind],
            anchor,
            fluid_cells,
            cut,
        );
        let l_class = lap.elapsed();
        mask.compute_merging(None);
        if profile_on() {
            eprintln!(
                "    from_cut laps: classification {:.0} ms, merging {:.0} ms",
                l_class.as_secs_f64() * 1e3,
                (lap.elapsed() - l_class).as_secs_f64() * 1e3
            );
        }
        Ok(mask)
    }

    /// The mask struct of a cut classification (every closure flag at its
    /// default; the caller sets them).
    pub(super) fn from_parts(
        g: Grid,
        periodic: bool,
        cell_fluid: Vec<bool>,
        kinds: [Vec<FaceKind>; 3],
        anchor: usize,
        fluid_cells: usize,
        cut: CutGeometry,
    ) -> Self {
        let [u_kind, v_kind, w_kind] = kinds;
        Self {
            grid: g,
            periodic_z: periodic,
            cell_fluid,
            u_kind,
            v_kind,
            w_kind,
            u_ghosts: Vec::new(),
            v_ghosts: Vec::new(),
            w_ghosts: Vec::new(),
            anchor,
            fluid_cells,
            cut: Some(cut),
            step_apertures: None,
            step_open: None,
            merge_master: Vec::new(),
            scheme: crate::solvers::incompressible::ConvectionScheme::Upwind,
            density: 1.0,
            wall_order: 1,
            wall_distance_oblique: false,
            diffusion_transverse: false,
            distance_floor_fine: false,
            wall_advancing: false,
            exchange_convection_off: false,
            wall_exchange_axis: false,
            changed_cells: None,
            cv_sides_exact: false,
            wall_order2_centroid: false,
            wall_exchange_foot: false,
            conv_sides_exact: false,
            wall_flux_true_normal: false,
            wall_foot_centroid: false,
            grad_weights: None,
            diffusion_centroid: false,
            face_shifts: None,
        }
    }

    /// The virtual merging map: a small cell (fraction < `MERGE_FRACTION`
    /// at both ends of the step) takes as master its active face
    /// neighbour of largest fraction that is not small itself; a small
    /// cell without such a neighbour keeps its own row.
    pub(super) fn compute_merging(&mut self, old: Option<&Mask>) {
        let Some(cut) = self.cut.as_ref() else {
            return;
        };
        let g = self.grid;
        let n = g.cells();
        let lat = self.lattice();
        let frac = |idx: usize| {
            let v = cut.vol[idx];
            old.and_then(|o| o.cut.as_ref())
                .map_or(v, |oc| v.max(oc.vol[idx]))
        };
        // `RTX_E3_MERGE_FRACTION` overrides the threshold (a study knob).
        let threshold = std::env::var("RTX_E3_MERGE_FRACTION")
            .ok()
            .and_then(|v| v.parse().ok())
            .unwrap_or(MERGE_FRACTION);
        let small: Vec<bool> = {
            use rayon::prelude::*;
            (0..n)
                .into_par_iter()
                .map(|idx| self.cell_active(idx) && frac(idx) < threshold)
                .collect()
        };
        let mut master = vec![usize::MAX; n];
        for idx in (0..n).filter(|&i| small[i]) {
            let (k, j, i) = g.kji(idx);
            let p = [i as i64, j as i64, k as i64];
            let mut best: Option<(f64, usize)> = None;
            for d in 0..3 {
                for side in [-1i64, 1] {
                    let mut q = p;
                    q[d] += side;
                    if let Some(nb) = lat.cell(q) {
                        let v = cut.vol[nb];
                        if self.cell_active(nb) && !small[nb] && best.is_none_or(|b| v > b.0) {
                            best = Some((v, nb));
                        }
                    }
                }
            }
            if let Some((_, m)) = best {
                master[idx] = m;
            }
        }
        self.merge_master = master;
    }

    /// Set the step-averaged apertures and the space-time classification
    /// from the previous mask's geometry: the trapezoid `½(αⁿ + αⁿ⁺¹)`,
    /// or with `inner` intermediate geometries the composite trapezoid
    /// over the step (the exact time integral of the space-time cut cell,
    /// arXiv 2512.23358, approached as the sub-sampling refines).
    pub fn set_step_apertures(&mut self, old: &mut Mask) {
        self.set_step_apertures_with(old, &[]);
    }

    /// As [`Self::set_step_apertures`] with the apertures of the
    /// intermediate geometries `inner` (in time order) inside the step.
    pub fn set_step_apertures_with(&mut self, old: &mut Mask, inner: &[&CutGeometry]) {
        let (Some(cut), Some(old_cut)) = (self.cut.as_ref(), old.cut.as_ref()) else {
            return;
        };
        let lap = std::time::Instant::now();
        let n = inner.len() + 1;
        let w_end = 0.5 / n as f64;
        let w_in = 1.0 / n as f64;
        use rayon::prelude::*;
        let g = self.grid;
        let (nx, ny, nz) = (g.nx, g.ny, g.nz);
        let nxy = nx * ny;
        // P1-4: the changed set — cells touched by either build (this step's
        // apertures average the two geometries), dilated by one.
        let t_new = &cut.touched_cell;
        let t_old = &old_cut.touched_cell;
        let periodic = self.periodic_z;
        let changed: Vec<usize> = (0..g.cells())
            .into_par_iter()
            .filter(|&idx| {
                let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx);
                let t = |q: usize| t_new[q] || t_old[q];
                if t(idx) {
                    return true;
                }
                (i + 1 < nx && t(idx + 1))
                    || (i > 0 && t(idx - 1))
                    || (j + 1 < ny && t(idx + nx))
                    || (j > 0 && t(idx - nx))
                    || (k + 1 < nz && t(idx + nxy))
                    || (k > 0 && t(idx - nxy))
                    || (periodic && nz > 1 && k + 1 == nz && t(idx - (nz - 1) * nxy))
                    || (periodic && nz > 1 && k == 0 && t(idx + (nz - 1) * nxy))
            })
            .collect();
        let l_changed = lap.elapsed();
        // P1-5 (j): with the trapezoid alone and a previous step's arrays, a
        // face of no changed cell keeps its step aperture (its corners were
        // untouched in both builds, so αⁿ⁻¹ = αⁿ = αⁿ⁺¹): the old mask's
        // arrays move over and only the changed cells' faces and the
        // changed cells' activity are recomputed, with the same expressions.
        let prev = if inner.is_empty() {
            old.step_apertures.take().zip(old.step_open.take())
        } else {
            None
        };
        if let Some(((mut au, mut av, mut aw), (mut ou, mut ov, mut ow, mut active))) = prev {
            let mix = |a: &[f64], b: &[f64], f: usize| w_end * (a[f] + b[f]);
            for &idx in &changed {
                let (k, j, i) = (idx / nxy, (idx % nxy) / nx, idx % nx);
                for f in [g.uface(k, j, i), g.uface(k, j, i + 1)] {
                    au[f] = mix(&cut.a_u, &old_cut.a_u, f);
                    ou[f] = au[f] > 0.0;
                }
                for f in [g.vface(k, j, i), g.vface(k, j + 1, i)] {
                    av[f] = mix(&cut.a_v, &old_cut.a_v, f);
                    ov[f] = av[f] > 0.0;
                }
                for f in [g.wface(k, j, i), g.wface(k + 1, j, i)] {
                    aw[f] = mix(&cut.a_w, &old_cut.a_w, f);
                    ow[f] = aw[f] > 0.0;
                }
                active[idx] = self.cell_fluid[idx] || old.cell_fluid[idx];
            }
            self.step_open = Some((ou, ov, ow, active));
            self.step_apertures = Some((au, av, aw));
        } else {
            let avg = |pick: &dyn Fn(&CutGeometry) -> &[f64]| -> Vec<f64> {
                let a = pick(cut);
                let b = pick(old_cut);
                let mut out: Vec<f64> = a.iter().zip(b).map(|(x, y)| w_end * (x + y)).collect();
                for gi in inner {
                    for (o, v) in out.iter_mut().zip(pick(gi)) {
                        *o += w_in * v;
                    }
                }
                out
            };
            let au = avg(&|gg: &CutGeometry| &gg.a_u);
            let av = avg(&|gg: &CutGeometry| &gg.a_v);
            let aw = avg(&|gg: &CutGeometry| &gg.a_w);
            let open = |a: &[f64]| -> Vec<bool> { a.iter().map(|&x| x > 0.0).collect() };
            let active = self
                .cell_fluid
                .iter()
                .zip(&old.cell_fluid)
                .map(|(&n, &o)| n || o)
                .collect();
            self.step_open = Some((open(&au), open(&av), open(&aw), active));
            self.step_apertures = Some((au, av, aw));
        }
        let l_apert = lap.elapsed();
        self.compute_merging(Some(old));
        if profile_on() {
            eprintln!(
                "    step-aperture laps: changed set {:.0} ms, apertures {:.0} ms, merging {:.0} ms",
                l_changed.as_secs_f64() * 1e3,
                (l_apert - l_changed).as_secs_f64() * 1e3,
                (lap.elapsed() - l_apert).as_secs_f64() * 1e3
            );
        }
        self.changed_cells = Some(changed);
    }

    pub(super) fn lattice(&self) -> Lattice {
        Lattice {
            g: self.grid,
            periodic_z: self.periodic_z,
        }
    }

    /// Aperture of the face of component `c` at lattice `p` (1 without a
    /// cut geometry), `None` outside the grid.
    pub(super) fn aperture(&self, c: usize, p: [i64; 3]) -> Option<f64> {
        let f = self.lattice().face(c, p)?;
        Some(match c {
            0 => self.a_u(f),
            1 => self.a_v(f),
            _ => self.a_w(f),
        })
    }

    /// The momentum control volume of the unknown face of component `c` at
    /// `p`: its face apertures are the averages of the two adjacent cells'
    /// (the own-direction faces at the cell centres average the face's and
    /// its own-direction neighbours' apertures), its wall closes them, its
    /// wall distance is the face centre's signed distance moved to the
    /// fluid part's centre, `φ + ½h(1 − α)`, floored.
    pub(super) fn cv_geometry(&self, c: usize, p: [i64; 3]) -> CvGeometry {
        let g = self.grid;
        let h = [g.dx, g.dy, g.dz];
        let area = [g.dy * g.dz, g.dx * g.dz, g.dx * g.dy];
        let mut e = |d: usize| {
            let mut v = [0i64; 3];
            v[d] = 1;
            v
        };
        let add =
            |a: [i64; 3], b: [i64; 3], s: i64| [a[0] + s * b[0], a[1] + s * b[1], a[2] + s * b[2]];
        let ec = e(c);
        let cell_minus = add(p, ec, -1);
        let cell_plus = p;
        let alpha = self.aperture(c, p).unwrap_or(1.0);
        let mut ap = [[1.0; 2]; 3];
        for d in 0..3 {
            let ed = e(d);
            if d == c {
                let am = self.aperture(c, add(p, ec, -1)).unwrap_or(alpha);
                let apl = self.aperture(c, add(p, ec, 1)).unwrap_or(alpha);
                ap[d] = [0.5 * (am + alpha), 0.5 * (alpha + apl)];
            } else {
                let minus = 0.5
                    * (self.aperture(d, cell_minus).unwrap_or(1.0)
                        + self.aperture(d, cell_plus).unwrap_or(1.0));
                let plus = 0.5
                    * (self.aperture(d, add(cell_minus, ed, 1)).unwrap_or(1.0)
                        + self.aperture(d, add(cell_plus, ed, 1)).unwrap_or(1.0));
                ap[d] = [minus, plus];
            }
        }
        // S2-7: the sides' own apertures instead of the whole-face averages.
        if self.cv_sides_exact {
            if let Some(exact) = self
                .cut
                .as_ref()
                .and_then(|cut| self.exact_cv_sides(cut, c, p))
            {
                ap = exact;
            }
        }
        let mut wall = [0.0; 3];
        for d in 0..3 {
            wall[d] = -(ap[d][1] - ap[d][0]) * area[d];
        }
        let h_min = g.dx.min(g.dy).min(g.dz);
        let phi_face = self.cut.as_ref().map_or(h_min, |cut| {
            let f = self
                .lattice()
                .face(c, p)
                .expect("unknown face inside the grid");
            match c {
                0 => cut.d_u[f],
                1 => cut.d_v[f],
                _ => cut.d_w[f],
            }
        });
        // The open part's centroid sits ½h(1 − α) from the face centre IN THE
        // FACE PLANE: its wall distance gains that times the wall normal's
        // in-plane part (1 for a wall parallel to the face normal).
        let n_t = if self.wall_distance_oblique {
            let a_w = (wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]).sqrt();
            if a_w > 0.0 {
                let n_c = wall[c] / a_w;
                (1.0 - n_c * n_c).max(0.0).sqrt()
            } else {
                1.0
            }
        } else {
            1.0
        };
        let distance =
            (phi_face + 0.5 * h[c] * (1.0 - alpha) * n_t).max(self.distance_floor() * h_min);
        CvGeometry {
            alpha,
            ap,
            wall,
            distance,
        }
    }

    /// S2-7: the control volume's side apertures from the interpolant on the
    /// sides' OWN corners. An unknown face's control volume is the tile
    /// between the two adjacent cells' centres: across `c` its sides are
    /// two HALF faces (the far half of `cell_minus`'s face, the near half
    /// of `cell_plus`'s), in the own direction the two cells' centre
    /// planes. φ is linear along every edge, so the mid-edge values are
    /// exact for the interpolant; each half face / centre plane is a quad
    /// through the faces' own `quad_fraction`. The averages of whole-face
    /// apertures the default takes are wrong by O(1) wherever the wall
    /// crosses a side (the in-plane momentum residual on oblique walls).
    /// `None` at a domain side (the default stays).
    fn exact_cv_sides(&self, cut: &CutGeometry, c: usize, p: [i64; 3]) -> Option<[[f64; 2]; 3]> {
        let parts = self.exact_cv_side_parts(cut, c, p)?;
        let mut ap = [[1.0; 2]; 3];
        for d in 0..3 {
            for side in 0..2 {
                ap[d][side] = if d == c {
                    parts[d][side].0
                } else {
                    0.5 * (parts[d][side].0 + parts[d][side].1)
                };
            }
        }
        Some(ap)
    }

    /// The parts of [`Self::exact_cv_sides`]: across `c`, the (far half of
    /// `cell_minus`'s face, near half of `cell_plus`'s face) apertures per
    /// side; along `c`, the centre plane's aperture (twice).
    pub(super) fn exact_cv_side_parts(
        &self,
        cut: &CutGeometry,
        c: usize,
        p: [i64; 3],
    ) -> Option<[[(f64, f64); 2]; 3]> {
        let lat = self.lattice();
        let g = self.grid;
        let mut pm = p;
        pm[c] -= 1;
        let cells = [g.kji(lat.cell(pm)?), g.kji(lat.cell(p)?)];
        // The corner of cell (k, j, i) at unit offsets `o = [di, dj, dk]`.
        let corner = |cell: (usize, usize, usize), o: [usize; 3]| {
            cut.corner_phi(cell.0 + o[2], cell.1 + o[1], cell.2 + o[0])
        };
        // The value at a cell's corner or, with `half`, at the mid-point of
        // its edge along `c` (the interpolant's mean of the two corners).
        let value = |cell: (usize, usize, usize), mut o: [usize; 3], half: bool| -> f64 {
            if half {
                o[c] = 0;
                let a = corner(cell, o);
                o[c] = 1;
                0.5 * (a + corner(cell, o))
            } else {
                corner(cell, o)
            }
        };
        let mut ap = [[(1.0, 1.0); 2]; 3];
        for d in 0..3 {
            if d == c {
                // The two cells' centre planes across `c`: corners at the
                // mid-points of the cells' `c` edges.
                let (d1, d2) = ((c + 1) % 3, (c + 2) % 3);
                for (side, cell) in cells.iter().enumerate() {
                    let mid = |o1: usize, o2: usize| {
                        let mut o = [0; 3];
                        o[d1] = o1;
                        o[d2] = o2;
                        value(*cell, o, true)
                    };
                    let a = super::cut::quad_fraction(mid(0, 0), mid(1, 0), mid(0, 1), mid(1, 1));
                    ap[d][side] = (a, a);
                }
            } else {
                let e = 3 - c - d;
                for side in 0..2 {
                    // cell_minus's `d` face at `side`, its half nearer the
                    // unknown face (c from ½ to 1); cell_plus's, c from 0 to ½.
                    let half = |cell: (usize, usize, usize), far: bool| -> f64 {
                        let at = |oc: u8, oe: usize| {
                            let mut o = [0; 3];
                            o[d] = side;
                            o[e] = oe;
                            match oc {
                                0 => value(cell, o, false),
                                1 => value(cell, o, true),
                                _ => {
                                    o[c] = 1;
                                    value(cell, o, false)
                                }
                            }
                        };
                        if far {
                            super::cut::quad_fraction(at(1, 0), at(2, 0), at(1, 1), at(2, 1))
                        } else {
                            super::cut::quad_fraction(at(0, 0), at(1, 0), at(0, 1), at(1, 1))
                        }
                    };
                    ap[d][side] = (half(cells[0], true), half(cells[1], false));
                }
            }
        }
        Some(ap)
    }

    /// The surface velocity component `c` at the foot of the normal from
    /// the face centre `x`. With a cut geometry the signed distance and
    /// the normal come from the geometry's own corner values (the
    /// trilinear interpolant at the face centre, its gradient by central
    /// differences of the neighbouring face centres) — one call of the
    /// body's velocity per face instead of eight of its distance.
    pub fn surface_velocity_at(&self, body: &Body, x: [f64; 3], c: usize, t: f64) -> f64 {
        let g = self.grid;
        let (s, n) = match self.cut.as_ref() {
            Some(cut) => {
                let (s, n) = self.interpolant_distance_and_normal(cut, x);
                (s, n)
            }
            None => {
                let eps = 1e-6 * g.dx.min(g.dy).min(g.dz);
                let s = body.phi(x[0], x[1], x[2], t);
                let (n1, n2, n3) = body.normal(x[0], x[1], x[2], t, eps);
                (s, [n1, n2, n3])
            }
        };
        let v = body.surface_velocity(x[0] - s * n[0], x[1] - s * n[1], x[2] - s * n[2], t);
        [v.0, v.1, v.2][c]
    }

    /// φ and its unit gradient at `x` from the trilinear interpolant of the
    /// corner values (the cut geometry's own surface).
    pub(super) fn interpolant_distance_and_normal(
        &self,
        cut: &CutGeometry,
        x: [f64; 3],
    ) -> (f64, [f64; 3]) {
        let g = self.grid;
        let (nx, ny, nz) = (g.nx as i64, g.ny as i64, g.nz as i64);
        let h = [g.dx, g.dy, g.dz];
        let node = |i: i64, j: i64, k: i64| -> f64 {
            let i = i.clamp(0, nx);
            let j = j.clamp(0, ny);
            let k = k.clamp(0, nz);
            cut.phi[((k * (ny + 1) + j) * (nx + 1) + i) as usize]
        };
        let gx = x[0] / h[0];
        let gy = x[1] / h[1];
        let gz = x[2] / h[2];
        let (i0, j0, k0) = (gx.floor() as i64, gy.floor() as i64, gz.floor() as i64);
        let (fx, fy, fz) = (gx - i0 as f64, gy - j0 as f64, gz - k0 as f64);
        // Trilinear value and its partial derivatives.
        let c = |di: i64, dj: i64, dk: i64| node(i0 + di, j0 + dj, k0 + dk);
        let lerp = |a: f64, b: f64, f: f64| a + f * (b - a);
        let c00 = lerp(c(0, 0, 0), c(1, 0, 0), fx);
        let c10 = lerp(c(0, 1, 0), c(1, 1, 0), fx);
        let c01 = lerp(c(0, 0, 1), c(1, 0, 1), fx);
        let c11 = lerp(c(0, 1, 1), c(1, 1, 1), fx);
        let c0 = lerp(c00, c10, fy);
        let c1 = lerp(c01, c11, fy);
        let s = lerp(c0, c1, fz);
        let dx0 = lerp(c(1, 0, 0) - c(0, 0, 0), c(1, 1, 0) - c(0, 1, 0), fy);
        let dx1 = lerp(c(1, 0, 1) - c(0, 0, 1), c(1, 1, 1) - c(0, 1, 1), fy);
        let dphi_dx = lerp(dx0, dx1, fz) / h[0];
        let dy0 = lerp(c(0, 1, 0) - c(0, 0, 0), c(1, 1, 0) - c(1, 0, 0), fx);
        let dy1 = lerp(c(0, 1, 1) - c(0, 0, 1), c(1, 1, 1) - c(1, 0, 1), fx);
        let dphi_dy = lerp(dy0, dy1, fz) / h[1];
        let dz0 = lerp(c(0, 0, 1) - c(0, 0, 0), c(1, 0, 1) - c(1, 0, 0), fx);
        let dz1 = lerp(c(0, 1, 1) - c(0, 1, 0), c(1, 1, 1) - c(1, 1, 0), fx);
        let dphi_dz = lerp(dz0, dz1, fy) / h[2];
        let norm = (dphi_dx * dphi_dx + dphi_dy * dphi_dy + dphi_dz * dphi_dz).sqrt();
        if norm > 0.0 {
            (s, [dphi_dx / norm, dphi_dy / norm, dphi_dz / norm])
        } else {
            (s, [1.0, 0.0, 0.0])
        }
    }

    /// The volume fluxes of the surface velocity through every cell's wall
    /// into the body, `U_b·W_c` (zero for a body at rest; the porous
    /// manufactured surface's flux otherwise), made compatible: the net
    /// flux (the quadrature's defect on a closed surface — a rigid
    /// translation's is zero by closure) is redistributed over the wall
    /// cells by wall area, as the binary wall's ghost fluxes are. Returns
    /// the table and the correction (flux per unit wall area).
    pub fn wall_flux_table(&self, body: &Body, t: f64) -> (Vec<f64>, f64) {
        let mut table = vec![0.0; self.grid.cells()];
        let Some(cut) = self.cut.as_ref() else {
            return (table, 0.0);
        };
        let (mut net, mut area) = (0.0, 0.0);
        for idx in 0..table.len() {
            if !self.cell_fluid[idx] {
                continue;
            }
            let w = cut.wall[idx];
            let a = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
            if a == 0.0 {
                continue;
            }
            table[idx] = self.wall_flux(body, idx, t);
            net += table[idx];
            area += a;
        }
        let correction = if area > 0.0 { net / area } else { 0.0 };
        if correction != 0.0 {
            for (idx, w) in cut.wall.iter().enumerate() {
                let a = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
                table[idx] -= correction * a;
            }
        }
        (table, correction)
    }

    /// The moving rigid body's wall fluxes by the discrete geometric
    /// conservation law: `(V_c^{n+1} − V_c^n)/dt` per active cell (a dying
    /// cell's remaining volume leaves through its step-averaged apertures),
    /// the net (the cut geometry's closure defect) redistributed over the
    /// wall cells by wall area.
    pub fn gcl_flux_table(&self, old: &Mask, dt: f64) -> (Vec<f64>, f64) {
        let mut table = vec![0.0; self.grid.cells()];
        let (Some(cut), Some(old_cut)) = (self.cut.as_ref(), old.cut.as_ref()) else {
            return (table, 0.0);
        };
        let lap = std::time::Instant::now();
        let g = self.grid;
        let dv = g.dx * g.dy * g.dz;
        let (mut net, mut area) = (0.0, 0.0);
        // P1-5 (b): a cell untouched by this build has the old volume (its
        // entry is exactly 0 and adds nothing to `net`); a cell without a
        // wall adds nothing to `area`. The sums run over the changed cells
        // and the wall cells in ascending order — the serial loop's order
        // with its zero terms left out, bit for bit.
        let wall_cells: Vec<usize> = {
            use rayon::prelude::*;
            (0..table.len())
                .into_par_iter()
                .filter(|&idx| self.cell_active(idx) && cut.wall[idx] != [0.0; 3])
                .collect()
        };
        let cells: Vec<usize> = match self.changed_cells.as_deref() {
            Some(changed) => changed.to_vec(),
            None => (0..table.len()).collect(),
        };
        for &idx in &cells {
            if !self.cell_active(idx) {
                continue;
            }
            let entry = (cut.vol[idx] - old_cut.vol[idx]) * dv / dt;
            table[idx] = entry;
            net += entry;
        }
        for &idx in &wall_cells {
            let w = cut.wall[idx];
            area += (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
        }
        if std::env::var_os("RTX_E3_DEBUG").is_some() {
            let dead = (0..table.len())
                .filter(|&i| !self.cell_fluid[i] && old.cell_fluid[i])
                .count();
            let fresh = (0..table.len())
                .filter(|&i| self.cell_fluid[i] && !old.cell_fluid[i])
                .count();
            let (vn, vn1): (f64, f64) = (old_cut.vol.iter().sum(), cut.vol.iter().sum());
            eprintln!(
                "    gcl: dead {dead} fresh {fresh} net {net:.3e} area {area:.3e} ΣV old {vn:.6} new {vn1:.6} (Δ {:.3e})",
                vn1 - vn
            );
        }
        let l_sums = lap.elapsed();
        let correction = if area > 0.0 { net / area } else { 0.0 };
        if correction != 0.0 {
            // Every cell with a wall (the others subtract exactly 0).
            for idx in 0..table.len() {
                let w = cut.wall[idx];
                if w != [0.0; 3] {
                    let a = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
                    table[idx] -= correction * a;
                }
            }
        }
        if profile_on() {
            eprintln!(
                "    gcl laps: lists + sums {:.0} ms, correction {:.0} ms",
                l_sums.as_secs_f64() * 1e3,
                (lap.elapsed() - l_sums).as_secs_f64() * 1e3
            );
        }
        (table, correction)
    }

    /// The volume flux of the surface velocity through a cell's wall into
    /// the body, `U_b·W_c`, uncorrected. Zero without a cut geometry.
    pub fn wall_flux(&self, body: &Body, idx: usize, t: f64) -> f64 {
        let Some(cut) = self.cut.as_ref() else {
            return 0.0;
        };
        let w = cut.wall[idx];
        if w == [0.0; 3] {
            return 0.0;
        }
        let g = self.grid;
        let (k, j, i) = g.kji(idx);
        let x = [
            (i as f64 + 0.5) * g.dx,
            (j as f64 + 0.5) * g.dy,
            (k as f64 + 0.5) * g.dz,
        ];
        // S2-7b A1: the flux through the TRUE surface at the wall foot —
        // `A_w (u_b · n)` with n the body's own normal into the body.
        if self.wall_flux_true_normal {
            let (s, n) = self.interpolant_distance_and_normal(cut, x);
            let foot = [x[0] - s * n[0], x[1] - s * n[1], x[2] - s * n[2]];
            let v = body.surface_velocity(foot[0], foot[1], foot[2], t);
            let eps = 1e-6 * g.dx.min(g.dy).min(g.dz);
            let (n1, n2, n3) = body.normal(foot[0], foot[1], foot[2], t, eps);
            let a_w = (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
            // `normal` points out of the solid: into the body is its negative.
            return -a_w * (v.0 * n1 + v.1 * n2 + v.2 * n3);
        }
        let mut flux = 0.0;
        for (c, wc) in w.iter().enumerate() {
            flux += self.surface_velocity_at(body, x, c, t) * wc;
        }
        flux
    }

    /// The cut-cell load route split into its pressure and shear parts.
    pub fn cut_wall_force_parts(
        &self,
        body: &Body,
        f: &Field,
        mu: f64,
        t: f64,
    ) -> Option<([f64; 3], [f64; 3])> {
        self.cut_wall_force_parts_in(body, f, mu, t, None)
    }

    pub(super) fn cut_wall_force_parts_in(
        &self,
        body: &Body,
        f: &Field,
        mu: f64,
        t: f64,
        planes: Option<(usize, usize)>,
    ) -> Option<([f64; 3], [f64; 3])> {
        let cut = self.cut.as_ref()?;
        let g = self.grid;
        let (nx, ny, nz) = (g.nx, g.ny, g.nz);
        let (k0, k1) = planes.unwrap_or((0, nz));
        let mut pressure = [0.0; 3];
        let mut force = [0.0; 3];
        for (idx, w) in cut.wall.iter().enumerate() {
            let (k, _, i) = g.kji(idx);
            if self.cell_fluid[idx] && k >= k0 && k < k1 && in_load_window((i as f64 + 0.5) * g.dx)
            {
                for c in 0..3 {
                    pressure[c] += f.p[idx] * w[c];
                }
            }
        }
        let gw = self.gradient_weight_force(&f.p, Some((k0, k1)));
        for c in 0..3 {
            pressure[c] += gw[c];
        }
        let lat = self.lattice();
        let values: [&[f64]; 3] = [&f.u, &f.v, &f.w];
        let w_range = if self.periodic_z { 0..nz } else { 1..nz };
        for c in 0..3 {
            let (ir, jr, kr) = match c {
                0 => (1..nx, 0..ny, k0..k1),
                1 => (0..nx, 1..ny, k0..k1),
                _ => (0..nx, 0..ny, w_range.start.max(k0)..w_range.end.min(k1)),
            };
            for k in kr {
                for j in jr.clone() {
                    for i in ir.clone() {
                        let p = [i as i64, j as i64, k as i64];
                        let idx = lat.face(c, p).expect("face");
                        let kind = match c {
                            0 => self.u_kind[idx],
                            1 => self.v_kind[idx],
                            _ => self.w_kind[idx],
                        };
                        if kind != FaceKind::Fluid || !in_load_window(lat.face_position(c, p)[0]) {
                            continue;
                        }
                        let cv = self.cv_geometry(c, p);
                        let a_w = (cv.wall[0] * cv.wall[0]
                            + cv.wall[1] * cv.wall[1]
                            + cv.wall[2] * cv.wall[2])
                            .sqrt();
                        if a_w == 0.0 {
                            continue;
                        }
                        let x = lat.face_position(c, p);
                        let ub = self.surface_velocity_at(body, x, c, t);
                        let xi = if self.wall_advancing {
                            self.surface_normal_velocity_at(body, x, t) * cv.distance * self.density
                                / mu
                        } else {
                            0.0
                        };
                        let (c1, c2, nb) = self.wall_gradient(c, p, &cv, xi);
                        let un = nb.map_or(ub, |f| values[c][f]);
                        force[c] += mu * a_w * (c1 * (values[c][idx] - ub) + c2 * (un - ub));
                    }
                }
            }
        }
        Some((pressure, force))
    }
}