//! R8-c: the fluid–structure interface's LOAD side for a 3D flag — the
//! cut-cell wall's force as the list of its summands (each located), and
//! their consistent, conservative distribution onto a structured Hex20
//! plate (the R8-b structure's node layout).
//!
//! # The loads
//!
//! [`Mask::cut_wall_loads`] returns every summand of the operator load
//! route `cut_wall_force` — the per-cell pressure `p_c W_c` (at the cell
//! centre), the per-face implicit wall shear and the per-face wall exchange
//! (diffusive and convective; at the fluid face's position) — each with
//! its FOOT on the body's surface (`x − φ n`, two projection steps on the
//! body's φ). The sums are the route's own: the same loops, observed.
//! The route's total is returned with them; the summands' sum differs from
//! it by the summation order only (round-off).
//!
//! # The transfer
//!
//! [`HexPlate::transfer`] hands each load `(point, F)` to the Hex20 element
//! of the deformed plate that contains `point` (the isoparametric map
//! inverted by Newton to round-off; the nearest element's extrapolation
//! when no element contains it) and distributes it with the element's
//! shape functions, `f_a = N_a(ξ) F`: the consistent nodal load of a point
//! force (its virtual work). The serendipity functions are a partition of
//! unity and reproduce the coordinates (`Σ N_a x_a = x(ξ) = point` once
//! Newton has converged), so the total force AND the total moment about
//! any point are conserved to round-off, whatever the element. A foot on
//! the plate's faces loads that face's nodes only; a foot inside the solid
//! (the capsule's rounded span edges and tip lie inside the Hex box) also
//! loads the element's other nodes — reported as the interior share.
//! The motion side is the transpose: [`HexPlate::locate`] gives the same
//! weights for the velocity `Σ N_a v_a` at a fluid point (work-conjugate).
//!
//! Host only (phase 1): the loads are read from the host mirror of the
//! device state at the sample steps, as the load routes are today.

use super::body::Body;
use super::field::Field;
use super::plate::PlateSurface;
use super::wall::Mask;

/// Which summand of the operator load route a [`WallLoad`] is.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LoadKind {
    /// `p_c W_c` of a cut cell.
    Pressure,
    /// The implicit wall shear of a fluid face.
    Shear,
    /// The diffusive exchange with a prescribed neighbour face.
    ExchangeDiffusive,
    /// The convective exchange with a prescribed neighbour face.
    ExchangeConvective,
}

/// One summand of the cut-cell wall force (a force ON the body).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct WallLoad {
    pub kind: LoadKind,
    /// Where the operator evaluates it: the cell centre or the face position.
    pub x: [f64; 3],
    /// Its foot on the body's surface (the application point).
    pub foot: [f64; 3],
    /// The force on the body.
    pub f: [f64; 3],
}

impl Mask {
    /// The operator load route's summands (see the module doc) and the
    /// route's total as `cut_wall_force` computes it. `None` without a cut
    /// geometry, or with the S2-5 gradient weights on (their force is not
    /// decomposed).
    pub fn cut_wall_loads(
        &self,
        body: &Body,
        f: &Field,
        mu: f64,
        t: f64,
    ) -> Option<(Vec<WallLoad>, [f64; 3])> {
        if self.grad_weights.is_some() {
            return None;
        }
        // The loops visit every fluid cell and face; only the non-zero
        // summands are loads (the zero ones change no sum).
        let mut raw: Vec<(LoadKind, [f64; 3], [f64; 3])> = Vec::new();
        let mut keep = |k: LoadKind, x: [f64; 3], v: [f64; 3]| {
            if v != [0.0; 3] {
                raw.push((k, x, v));
            }
        };
        let (p, s) = self.cut_wall_force_parts_sink(body, f, mu, t, None, &mut keep)?;
        let (d, c) =
            self.cut_wall_exchange_parts_sink(body, f, mu, self.density, t, None, &mut keep)?;
        let xch = [d[0] + c[0], d[1] + c[1], d[2] + c[2]];
        let total = [
            p[0] + s[0] + xch[0],
            p[1] + s[1] + xch[1],
            p[2] + s[2] + xch[2],
        ];
        let g = self.grid;
        let eps = 1e-6 * g.dx.min(g.dy).min(g.dz);
        let foot = |x: [f64; 3]| {
            let mut q = x;
            for _ in 0..2 {
                let s = body.phi(q[0], q[1], q[2], t);
                let n = body.normal(q[0], q[1], q[2], t, eps);
                q = [q[0] - s * n.0, q[1] - s * n.1, q[2] - s * n.2];
            }
            q
        };
        use rayon::prelude::*;
        let loads = raw
            .par_iter()
            .map(|&(kind, x, f)| WallLoad {
                kind,
                x,
                foot: foot(x),
                f,
            })
            .collect();
        Some((loads, total))
    }
}

/// The Hex20 node order of R8-b's `Flag3d` (rtx-fea `analysis/flag3d.rs`):
/// the natural coordinates `(ξ, η, ζ)` ↔ lattice `(i, j, k)` = (length,
/// thickness, span); corners of `ζ = −1` then `ζ = +1` counter-clockwise
/// from `(−1, −1)`, the mid-edges of `ζ = −1`, of `ζ = +1`, then the four
/// span edges.
const HEX20: [[i8; 3]; 20] = [
    [-1, -1, -1],
    [1, -1, -1],
    [1, 1, -1],
    [-1, 1, -1],
    [-1, -1, 1],
    [1, -1, 1],
    [1, 1, 1],
    [-1, 1, 1],
    [0, -1, -1],
    [1, 0, -1],
    [0, 1, -1],
    [-1, 0, -1],
    [0, -1, 1],
    [1, 0, 1],
    [0, 1, 1],
    [-1, 0, 1],
    [-1, -1, 0],
    [1, -1, 0],
    [1, 1, 0],
    [-1, 1, 0],
];

/// The Hex20 serendipity shape functions and their natural derivatives.
#[must_use]
pub fn hex20(xi: [f64; 3]) -> ([f64; 20], [[f64; 3]; 20]) {
    let mut n = [0.0; 20];
    let mut d = [[0.0; 3]; 20];
    for (a, c) in HEX20.iter().enumerate() {
        let ca = [f64::from(c[0]), f64::from(c[1]), f64::from(c[2])];
        let lin = |m: usize| 1.0 + xi[m] * ca[m];
        if c.iter().all(|&v| v != 0) {
            let (p, q, r) = (lin(0), lin(1), lin(2));
            let s = xi[0] * ca[0] + xi[1] * ca[1] + xi[2] * ca[2];
            n[a] = 0.125 * p * q * r * (s - 2.0);
            d[a][0] = 0.125 * q * r * ca[0] * (s - 2.0 + p);
            d[a][1] = 0.125 * p * r * ca[1] * (s - 2.0 + q);
            d[a][2] = 0.125 * p * q * ca[2] * (s - 2.0 + r);
        } else {
            // The zero direction `m0`; the other two linear.
            let m0 = c.iter().position(|&v| v == 0).expect("mid-edge");
            let (m1, m2) = ((m0 + 1) % 3, (m0 + 2) % 3);
            let bub = 1.0 - xi[m0] * xi[m0];
            n[a] = 0.25 * bub * lin(m1) * lin(m2);
            d[a][m0] = 0.25 * (-2.0 * xi[m0]) * lin(m1) * lin(m2);
            d[a][m1] = 0.25 * bub * ca[m1] * lin(m2);
            d[a][m2] = 0.25 * bub * lin(m1) * ca[m2];
        }
    }
    (n, d)
}

/// A structured Hex20 plate on the `(2nx+1) × (2ny+1) × (2nz+1)`
/// serendipity lattice — `nx` elements along the length, `ny` through the
/// thickness, `nz` along the span — numbered exactly as R8-b's `Flag3d`:
/// lattice points with at most one odd index, scanned `i` (length)
/// outermost, then `j` (thickness), then `k` (span); node `n` is the
/// `n`-th such point. The elements scan `(ex, ey, ez)` the same way.
#[derive(Debug, Clone)]
pub struct HexPlate {
    pub nx: usize,
    pub ny: usize,
    pub nz: usize,
    dims: [usize; 3],
    lattice: Vec<Option<usize>>,
    points: Vec<[usize; 3]>,
    elements: Vec<[usize; 20]>,
}

/// A point located in the plate: its element, natural coordinates, the
/// nodes and weights `N_a(ξ)`, the Newton residual `|x(ξ) − p|` and how
/// far outside the element it is (`max |ξ| − 1`, ≤ 0 inside).
#[derive(Debug, Clone)]
pub struct Location {
    pub element: usize,
    pub xi: [f64; 3],
    pub nodes: [usize; 20],
    pub weights: [f64; 20],
    pub residual: f64,
    pub outside: f64,
}

/// The natural-coordinate distance outside an element from which a load
/// counts as extrapolated in [`PlateTransfer`] (0.02 of the half-element).
pub const FAR_OUTSIDE: f64 = 0.02;

/// The result of one load transfer (see [`HexPlate::transfer`]).
#[derive(Debug, Clone)]
pub struct PlateTransfer {
    /// The nodal forces, in the node numbering.
    pub nodal: Vec<[f64; 3]>,
    /// Σ F of the loads and Σ f_a of the nodes.
    pub force_in: [f64; 3],
    pub force_out: [f64; 3],
    /// The moments about `origin`: Σ (p − o) × F and Σ (x_a − o) × f_a.
    pub moment_in: [f64; 3],
    pub moment_out: [f64; 3],
    /// The largest Newton residual (m) and the largest outside-ness.
    pub max_residual: f64,
    pub max_outside: f64,
    /// Loads located outside every element by more than [`FAR_OUTSIDE`]
    /// in natural coordinates (a foot off the plate by more than a few %
    /// of an element — the surfaces' round-off-level mismatch is excluded),
    /// their Σ|F|, and the worst one's point.
    pub extrapolated: usize,
    pub extrapolated_load: f64,
    pub worst_point: [f64; 3],
    /// The share of Σ|f_a| on nodes off the wetted surface (the interior
    /// layers and the clamped root face).
    pub interior_share: f64,
}

fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
    [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    ]
}

fn norm(a: [f64; 3]) -> f64 {
    (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt()
}

impl HexPlate {
    /// The lattice and the elements of an `nx × ny × nz` plate.
    #[must_use]
    pub fn new(nx: usize, ny: usize, nz: usize) -> Self {
        assert!(
            nx > 0 && ny > 0 && nz > 0,
            "element counts must be positive"
        );
        let dims = [2 * nx + 1, 2 * ny + 1, 2 * nz + 1];
        let mut lattice = vec![None; dims[0] * dims[1] * dims[2]];
        let mut points = Vec::new();
        for i in 0..dims[0] {
            for j in 0..dims[1] {
                for k in 0..dims[2] {
                    if (i % 2) + (j % 2) + (k % 2) > 1 {
                        continue;
                    }
                    lattice[(i * dims[1] + j) * dims[2] + k] = Some(points.len());
                    points.push([i, j, k]);
                }
            }
        }
        let mut plate = Self {
            nx,
            ny,
            nz,
            dims,
            lattice,
            points,
            elements: Vec::new(),
        };
        for ex in 0..nx {
            for ey in 0..ny {
                for ez in 0..nz {
                    let (a, b, c) = (2 * ex, 2 * ey, 2 * ez);
                    let e = HEX20.map(|o| {
                        let at = |base: usize, v: i8| (base as i64 + 1 + i64::from(v)) as usize;
                        plate
                            .lattice_node(at(a, o[0]), at(b, o[1]), at(c, o[2]))
                            .expect("serendipity node")
                    });
                    plate.elements.push(e);
                }
            }
        }
        plate
    }

    /// The node at lattice `(i, j, k)`, if the point carries one.
    #[must_use]
    pub fn lattice_node(&self, i: usize, j: usize, k: usize) -> Option<usize> {
        if i >= self.dims[0] || j >= self.dims[1] || k >= self.dims[2] {
            return None;
        }
        self.lattice[(i * self.dims[1] + j) * self.dims[2] + k]
    }

    /// The lattice point of node `n`.
    #[must_use]
    pub fn lattice_of(&self, n: usize) -> [usize; 3] {
        self.points[n]
    }

    #[must_use]
    pub fn node_count(&self) -> usize {
        self.points.len()
    }

    #[must_use]
    pub fn elements(&self) -> &[[usize; 20]] {
        &self.elements
    }

    /// Whether node `n` lies on the wetted surface (the two faces, the tip,
    /// the span edges; not the clamped root face `i = 0` unless also on one
    /// of those).
    #[must_use]
    pub fn is_wetted(&self, n: usize) -> bool {
        let [i, j, k] = self.points[n];
        j == 0 || j == 2 * self.ny || i == 2 * self.nx || k == 0 || k == 2 * self.nz
    }

    /// Node positions from a placement `(s, η, ζ) → x` of the lattice's
    /// fractions: `s = i/(2nx) ∈ [0, 1]` along the length, `η = j/ny − 1 ∈
    /// [−1, 1]` through the thickness, `ζ = k/(2nz) ∈ [0, 1]` along the span.
    #[must_use]
    pub fn place<F: Fn(f64, f64, f64) -> [f64; 3]>(&self, f: F) -> Vec<[f64; 3]> {
        self.points
            .iter()
            .map(|&[i, j, k]| {
                f(
                    i as f64 / (2 * self.nx) as f64,
                    j as f64 / self.ny as f64 - 1.0,
                    k as f64 / (2 * self.nz) as f64,
                )
            })
            .collect()
    }

    /// Newton on the isoparametric map of element `e` for the point `p`.
    fn invert(&self, pos: &[[f64; 3]], e: usize, p: [f64; 3]) -> Location {
        let nodes = self.elements[e];
        let x: Vec<[f64; 3]> = nodes.iter().map(|&n| pos[n]).collect();
        let mut xi = [0.0f64; 3];
        let mut residual = f64::INFINITY;
        let mut weights = [0.0; 20];
        for it in 0..60 {
            let (n, d) = hex20(xi);
            let mut r = [-p[0], -p[1], -p[2]];
            let mut jac = [[0.0f64; 3]; 3];
            for a in 0..20 {
                for c in 0..3 {
                    r[c] += n[a] * x[a][c];
                    for m in 0..3 {
                        jac[c][m] += x[a][c] * d[a][m];
                    }
                }
            }
            weights = n;
            let rn = norm(r);
            // Converged: two more iterations past the first residual at
            // round-off level make the last one a no-op.
            if rn <= residual && rn < 1e-15 && it > 2 {
                residual = rn;
                break;
            }
            residual = rn;
            // δ = −J⁻¹ r by the adjugate.
            let det = jac[0][0] * (jac[1][1] * jac[2][2] - jac[1][2] * jac[2][1])
                - jac[0][1] * (jac[1][0] * jac[2][2] - jac[1][2] * jac[2][0])
                + jac[0][2] * (jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0]);
            if det == 0.0 || !det.is_finite() {
                break;
            }
            let inv = [
                [
                    jac[1][1] * jac[2][2] - jac[1][2] * jac[2][1],
                    jac[0][2] * jac[2][1] - jac[0][1] * jac[2][2],
                    jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1],
                ],
                [
                    jac[1][2] * jac[2][0] - jac[1][0] * jac[2][2],
                    jac[0][0] * jac[2][2] - jac[0][2] * jac[2][0],
                    jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2],
                ],
                [
                    jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0],
                    jac[0][1] * jac[2][0] - jac[0][0] * jac[2][1],
                    jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0],
                ],
            ];
            for m in 0..3 {
                let dm = (inv[m][0] * r[0] + inv[m][1] * r[1] + inv[m][2] * r[2]) / det;
                // Keep the iterate bounded (a far point in a distorted element).
                xi[m] = (xi[m] - dm).clamp(-4.0, 4.0);
            }
        }
        let outside = xi.iter().fold(f64::NEG_INFINITY, |m, v| m.max(v.abs())) - 1.0;
        Location {
            element: e,
            xi,
            nodes,
            weights,
            residual,
            outside,
        }
    }

    /// Locate `p` in the deformed plate `pos`: the containing element
    /// (Newton on the elements nearest by centroid), else the nearest
    /// element's extrapolation (the smallest outside-ness).
    #[must_use]
    pub fn locate(&self, pos: &[[f64; 3]], centroids: &[[f64; 3]], p: [f64; 3]) -> Location {
        let mut order: Vec<(f64, usize)> = centroids
            .iter()
            .enumerate()
            .map(|(e, c)| {
                let d = [c[0] - p[0], c[1] - p[1], c[2] - p[2]];
                (d[0] * d[0] + d[1] * d[1] + d[2] * d[2], e)
            })
            .collect();
        let take = 8.min(order.len());
        order.select_nth_unstable_by(take - 1, |a, b| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1)));
        order[..take].sort_by(|a, b| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1)));
        let mut best: Option<Location> = None;
        for &(_, e) in &order[..take] {
            let loc = self.invert(pos, e, p);
            if loc.outside <= 1e-9 && loc.residual < 1e-12 {
                return loc;
            }
            if best
                .as_ref()
                .is_none_or(|b| (loc.outside, loc.residual) < (b.outside, b.residual))
            {
                best = Some(loc);
            }
        }
        best.expect("an element")
    }

    /// The elements' corner centroids in the deformed plate.
    #[must_use]
    pub fn centroids(&self, pos: &[[f64; 3]]) -> Vec<[f64; 3]> {
        self.elements
            .iter()
            .map(|e| {
                let mut c = [0.0; 3];
                for &n in &e[..8] {
                    for m in 0..3 {
                        c[m] += 0.125 * pos[n][m];
                    }
                }
                c
            })
            .collect()
    }

    /// The consistent nodal load of point forces `(p, F)` on the deformed
    /// plate `pos` (see the module doc), with the conservation budget about
    /// `origin`.
    #[must_use]
    pub fn transfer(
        &self,
        pos: &[[f64; 3]],
        loads: &[([f64; 3], [f64; 3])],
        origin: [f64; 3],
    ) -> PlateTransfer {
        use rayon::prelude::*;
        assert_eq!(pos.len(), self.node_count(), "one position per node");
        let centroids = self.centroids(pos);
        let located: Vec<Location> = loads
            .par_iter()
            .map(|&(p, _)| self.locate(pos, &centroids, p))
            .collect();
        let mut nodal = vec![[0.0f64; 3]; self.node_count()];
        let (mut force_in, mut moment_in) = ([0.0f64; 3], [0.0f64; 3]);
        let (mut max_residual, mut max_outside) = (0.0f64, f64::NEG_INFINITY);
        let (mut extrapolated, mut extrapolated_load) = (0usize, 0.0f64);
        let mut worst_point = [0.0f64; 3];
        for (&(p, f), loc) in loads.iter().zip(&located) {
            for (&n, &w) in loc.nodes.iter().zip(&loc.weights) {
                for c in 0..3 {
                    nodal[n][c] += w * f[c];
                }
            }
            let m = cross([p[0] - origin[0], p[1] - origin[1], p[2] - origin[2]], f);
            for c in 0..3 {
                force_in[c] += f[c];
                moment_in[c] += m[c];
            }
            max_residual = max_residual.max(loc.residual);
            if loc.outside > max_outside {
                max_outside = loc.outside;
                worst_point = p;
            }
            if loc.outside > FAR_OUTSIDE {
                extrapolated += 1;
                extrapolated_load += norm(f);
            }
        }
        let (mut force_out, mut moment_out) = ([0.0f64; 3], [0.0f64; 3]);
        let (mut total_abs, mut interior_abs) = (0.0f64, 0.0f64);
        for (n, f) in nodal.iter().enumerate() {
            let x = pos[n];
            let m = cross([x[0] - origin[0], x[1] - origin[1], x[2] - origin[2]], *f);
            for c in 0..3 {
                force_out[c] += f[c];
                moment_out[c] += m[c];
            }
            let a = norm(*f);
            total_abs += a;
            if !self.is_wetted(n) || self.points[n][0] == 0 {
                interior_abs += a;
            }
        }
        PlateTransfer {
            nodal,
            force_in,
            force_out,
            moment_in,
            moment_out,
            max_residual,
            max_outside,
            extrapolated,
            extrapolated_load,
            worst_point,
            interior_share: if total_abs > 0.0 {
                interior_abs / total_abs
            } else {
                0.0
            },
        }
    }

    /// The mid-surface of the deformed plate for the fluid's body
    /// ([`PlateSurface`]): the lattice's middle layer `j = ny` (`ny` even:
    /// a node layer), the stations at the element boundaries `k` even (their
    /// z the nodes' at the root, `i = 0`), the points along the length every
    /// lattice step; the last point of each station pulled back along its
    /// last segment by `tip_inset` (the capsule's apex on the structure's
    /// tip: the flag test's `RTX_E3_FLAG_TIP_INSET`, one half-thickness).
    /// With nodal velocities, the matching station velocities.
    #[must_use]
    pub fn mid_surface(
        &self,
        pos: &[[f64; 3]],
        vel: Option<&[[f64; 3]]>,
        tip_inset: f64,
    ) -> PlateSurface {
        assert!(
            self.ny % 2 == 0,
            "the mid-surface is a node layer for an even thickness count"
        );
        let j = self.ny;
        let ns = self.dims[0];
        let mut z = Vec::new();
        let mut xy = Vec::new();
        let mut vv = Vec::new();
        for k in (0..self.dims[2]).step_by(2) {
            z.push(pos[self.lattice_node(0, j, k).expect("root node")][2]);
            let row0 = xy.len();
            for i in 0..ns {
                let n = self.lattice_node(i, j, k).expect("mid-layer node");
                xy.push([pos[n][0], pos[n][1]]);
                if let Some(v) = vel {
                    vv.push([v[n][0], v[n][1]]);
                }
            }
            if tip_inset > 0.0 {
                let [ax, ay] = xy[row0 + ns - 2];
                let [bx, by] = xy[row0 + ns - 1];
                let len = ((bx - ax).powi(2) + (by - ay).powi(2)).sqrt();
                let f = (1.0 - tip_inset / len).max(0.0);
                xy[row0 + ns - 1] = [ax + f * (bx - ax), ay + f * (by - ay)];
            }
        }
        PlateSurface { z, ns, xy, vel: vv }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn hex20_is_a_partition_of_unity_with_nodal_interpolation() {
        for (a, c) in HEX20.iter().enumerate() {
            let (n, _) = hex20([f64::from(c[0]), f64::from(c[1]), f64::from(c[2])]);
            for (b, v) in n.iter().enumerate() {
                assert!((v - if a == b { 1.0 } else { 0.0 }).abs() < 1e-15);
            }
        }
        let (n, d) = hex20([0.3, -0.7, 0.45]);
        assert!((n.iter().sum::<f64>() - 1.0).abs() < 1e-15);
        for m in 0..3 {
            assert!(d.iter().map(|v| v[m]).sum::<f64>().abs() < 1e-14);
        }
        // Derivatives against central differences.
        let h = 1e-6;
        for m in 0..3 {
            let mut p = [0.3, -0.7, 0.45];
            let mut q = p;
            p[m] += h;
            q[m] -= h;
            let (np, _) = hex20(p);
            let (nq, _) = hex20(q);
            for a in 0..20 {
                assert!(((np[a] - nq[a]) / (2.0 * h) - d[a][m]).abs() < 1e-8);
            }
        }
    }

    #[test]
    fn layout_matches_the_structure() {
        let p = HexPlate::new(35, 2, 10);
        // Serendipity count: (2nx+1)(2ny+1)(2nz+1) minus the points with ≥ 2 odd.
        let full = 71 * 5 * 21;
        let two_odd = 35 * 2 * 21 + 35 * 5 * 10 + 71 * 2 * 10 - 2 * 35 * 2 * 10;
        assert_eq!(p.node_count(), full - two_odd);
        assert_eq!(p.elements().len(), 35 * 2 * 10);
        // Node 0 is (0,0,0), the next along the span.
        assert_eq!(p.lattice_of(0), [0, 0, 0]);
        assert_eq!(p.lattice_of(1), [0, 0, 1]);
        assert_eq!(p.elements()[0][1], p.lattice_node(2, 0, 0).unwrap());
    }

    /// A bent, twisted plate: arbitrary point loads inside and just outside
    /// are distributed with the total force and moment conserved to
    /// round-off.
    #[test]
    fn transfer_conserves_force_and_moment() {
        let plate = HexPlate::new(35, 2, 10);
        let bend = |s: f64, zeta: f64| 0.06 * s * s * (1.0 + 0.5 * (2.0 * zeta - 1.0));
        let pos = plate.place(|s, eta, zeta| {
            let x = 0.25 + 0.35 * s;
            let y = 0.2 + bend(s, zeta) + 0.01 * eta;
            [x, y, 0.105 + 0.2 * zeta]
        });
        let mut loads = Vec::new();
        let mut state = 12345u64;
        let mut rnd = || {
            state = state
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            (state >> 11) as f64 / (1u64 << 53) as f64
        };
        for _ in 0..3000 {
            let (s, zeta, eta) = (rnd(), rnd(), 2.4 * rnd() - 1.2);
            let p = [
                0.25 + 0.35 * s,
                0.2 + bend(s, zeta) + 0.01 * eta,
                0.105 + 0.2 * zeta,
            ];
            loads.push((p, [rnd() - 0.5, 10.0 * (rnd() - 0.5), 0.1 * (rnd() - 0.5)]));
        }
        let tr = plate.transfer(&pos, &loads, [0.25, 0.2, 0.205]);
        let scale_f = loads.iter().map(|l| norm(l.1)).sum::<f64>();
        for c in 0..3 {
            assert!((tr.force_in[c] - tr.force_out[c]).abs() < 1e-13 * scale_f);
            assert!((tr.moment_in[c] - tr.moment_out[c]).abs() < 1e-13 * scale_f);
        }
        assert!(tr.max_residual < 1e-14);
        assert!(tr.extrapolated > 0, "the ±1.2 band reaches outside");
    }

    /// The mid-surface of the deformed Hex plate, thickened by the half
    /// thickness, passes through the plate's face nodes (φ ≈ 0 there): the
    /// motion side's geometry and the structure agree.
    #[test]
    fn mid_surface_passes_the_face_nodes() {
        use super::super::body::DeviceSdf;
        let plate = HexPlate::new(35, 2, 10);
        let half = 0.01;
        let bend = |s: f64, zeta: f64| 0.08 * s * s * (1.0 + 0.3 * (2.0 * zeta - 1.0));
        let pos = plate.place(|s, eta, zeta| {
            // Offsets along the 3D normal of the mid-surface y = w(x, z).
            let ds = 1e-7;
            let (x, y) = (0.25 + 0.35 * s, 0.2 + bend(s, zeta));
            let wx = (bend(s + ds, zeta) - bend(s - ds, zeta)) / (2.0 * ds) / 0.35;
            let wz = (bend(s, zeta + ds) - bend(s, zeta - ds)) / (2.0 * ds) / 0.2;
            let r = (1.0 + wx * wx + wz * wz).sqrt();
            [
                x - half * eta * wx / r,
                y + half * eta / r,
                0.105 + 0.2 * zeta - half * eta * wz / r,
            ]
        });
        let surf = plate.mid_surface(&pos, None, 0.0);
        surf.validate().unwrap();
        let sdf = DeviceSdf {
            cyl: [-10.0, -10.0, 0.05],
            cyl_cut: false,
            flag_cut: false,
            zc: 0.205,
            span: 0.2,
            r_edge: 0.0,
            half,
            fillet: 0.0,
            poly: Vec::new(),
            vel: Vec::new(),
            plate: Some(surf),
        };
        let mut worst = 0.0f64;
        for n in 0..plate.node_count() {
            let [i, j, _] = plate.lattice_of(n);
            if (j == 0 || j == 4) && i > 1 && i < 69 {
                let p = pos[n];
                worst = worst.max(sdf.phi_host(p[0], p[1], p[2]).abs());
            }
        }
        // Measured 4.5e-7 m (the mid-surface's 5 mm chords, the first-order
        // slope correction) against the 10 mm half-thickness.
        assert!(
            worst < 2e-6,
            "face nodes off the level set by {worst:.3e} m"
        );
    }
}