//! Shared harness for the Turek–Hron FSI2 tests: the benchmark geometry,
//! the flag mesh, the wetted-interface bookkeeping, the embedded fluid
//! configuration, and the load sampling (spike clamp + optional surface
//! smoothing). `turek_hron_fsi2.rs` runs the coupled march on it;
//! `fsi2_interface_noise.rs` measures the continuity of one coupling pass
//! on the same machinery.
//!
//! Everything here is code motion from the tenth-session FSI2 test —
//! the physics and defaults are unchanged unless a test says otherwise.

#![allow(dead_code)] // several test crates share this; each uses a subset

pub mod march;
pub mod overset;
pub mod overset_march;
pub mod replay;
pub mod rescue;

use std::cell::Cell;
use std::sync::{Arc, RwLock};

use nalgebra::Vector3;
use rtx_cfd::CfdConfig;
use rtx_cfd::solvers::incompressible::{
    AleBoundaries, ConvectionScheme, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
    FlowField, PoissonSolverKind, PolygonSdf, SideBoundary, polygon_interface_velocity,
    polygon_signed_distance,
};
use rtx_fea::assembly::dof_mapping::DofComponent;
use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
use rtx_fsi::{FluidFace, WettedSurface, smooth_tractions};

pub const L: f64 = 2.5;
pub const H: f64 = 0.41;
pub const RHO_F: f64 = 1000.0;
pub const NU_F: f64 = 1e-3;
pub const U_MEAN: f64 = 1.0;
pub const RHO_S: f64 = 10_000.0;
pub const E_S: f64 = 1.4e6;
pub const NU_S: f64 = 0.4;

pub const FLAG_X0: f64 = 0.25;
pub const FLAG_X1: f64 = 0.6;
pub const FLAG_Y0: f64 = 0.19;
pub const FLAG_Y1: f64 = 0.21;

/// The parameters that distinguish the self-excited Turek–Hron cases on
/// the shared geometry: mean inflow (Re = 100 U), solid density and
/// Young's modulus. Everything else — channel, cylinder, flag, fluid —
/// is common.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct BenchmarkCase {
    pub name: &'static str,
    pub u_mean: f64,
    pub rho_s: f64,
    pub e_s: f64,
    pub nu_s: f64,
    /// The rigid-flag CFD drag on this geometry at this Re (the fluid
    /// harness check before anything couples): CFD2 / CFD3 means.
    pub rigid_drag_reference: f64,
}

/// FSI2: Re 100, density ratio 10 — the heavy flag's resonant flapping.
pub const FSI2: BenchmarkCase = BenchmarkCase {
    name: "FSI2",
    u_mean: 1.0,
    rho_s: 10_000.0,
    e_s: 1.4e6,
    nu_s: 0.4,
    rigid_drag_reference: 136.7,
};

/// FSI3: Re 200, density ratio 1 (mu_s = 2e6 → E = 5.6e6) — the
/// added-mass regime.
pub const FSI3: BenchmarkCase = BenchmarkCase {
    name: "FSI3",
    u_mean: 2.0,
    rho_s: 1_000.0,
    e_s: 5.6e6,
    nu_s: 0.4,
    rigid_drag_reference: 439.45,
};

/// `RTX_{prefix}_UMEAN` over the case's benchmark inflow — TWIN-1's
/// sweep parameter (`docs/twin_composition_campaign.md` in omni-cortex)
/// — and `RTX_{prefix}_ES` over the case's benchmark Young's modulus —
/// TWIN-2's second parameter (`docs/twin2_stiffness_campaign.md`).
/// Unset, the case comes back unchanged: the same f64s flow and the
/// march is digit-identical by construction (and verified in vivo on
/// both committed defaults). When set, the march's printed rigid-drag
/// reference still names the BENCHMARK value, which only applies at the
/// case's own inflow — the override lines below keep logs honest.
pub fn case_from_env(prefix: &str, case: BenchmarkCase) -> BenchmarkCase {
    let u_mean = env_or(&format!("RTX_{prefix}_UMEAN"), case.u_mean);
    if u_mean.to_bits() != case.u_mean.to_bits() {
        println!(
            "  {} u_mean OVERRIDDEN to {u_mean} (benchmark {}; Re = {:.0})",
            case.name,
            case.u_mean,
            100.0 * u_mean
        );
    }
    let e_s = env_or(&format!("RTX_{prefix}_ES"), case.e_s);
    if e_s.to_bits() != case.e_s.to_bits() {
        println!(
            "  {} e_s OVERRIDDEN to {e_s:.4e} (benchmark {:.4e}; E/E0 = {:.4})",
            case.name,
            case.e_s,
            e_s / case.e_s
        );
    }
    BenchmarkCase {
        u_mean,
        e_s,
        ..case
    }
}

pub fn circle_sdf(x: f64, y: f64) -> f64 {
    ((x - 0.2).powi(2) + (y - 0.2).powi(2)).sqrt() - 0.05
}

/// The ramped parabolic inflow of the benchmark definition, for a mean
/// inflow `u_mean`.
pub fn inflow_for(u_mean: f64, y: f64, t: f64) -> f64 {
    let ramp = if t < 2.0 {
        0.5 * (1.0 - (std::f64::consts::PI * t / 2.0).cos())
    } else {
        1.0
    };
    ramp * 1.5 * u_mean * y * (H - y) / (0.5 * H).powi(2)
}

/// FSI2's inflow.
pub fn inflow(y: f64, t: f64) -> f64 {
    inflow_for(U_MEAN, y, t)
}

pub fn env_or(name: &str, default: f64) -> f64 {
    std::env::var(name)
        .map(|v| v.parse().expect(name))
        .unwrap_or(default)
}

/// The flag's Quad8 mesh (as in FSI1 and the CSM tests).
pub fn flag_mesh(nx: usize, ny: usize) -> Mesh {
    let mut mesh = Mesh::new(2).unwrap();
    let (lx, ly) = (2 * nx + 1, 2 * ny + 1);
    let mut grid = vec![vec![None; ly]; lx];
    for (i, column) in grid.iter_mut().enumerate() {
        for (j, slot) in column.iter_mut().enumerate() {
            if i % 2 == 1 && j % 2 == 1 {
                continue;
            }
            let x = FLAG_X0 + (FLAG_X1 - FLAG_X0) * i as f64 / (2 * nx) as f64;
            let y = FLAG_Y0 + (FLAG_Y1 - FLAG_Y0) * j as f64 / (2 * ny) as f64;
            *slot = Some(mesh.add_node(Node::new_2d(x, y)));
        }
    }
    for i in 0..nx {
        for j in 0..ny {
            let (a, b) = (2 * i, 2 * j);
            let nodes = vec![
                grid[a][b].unwrap(),
                grid[a + 2][b].unwrap(),
                grid[a + 2][b + 2].unwrap(),
                grid[a][b + 2].unwrap(),
                grid[a + 1][b].unwrap(),
                grid[a + 2][b + 1].unwrap(),
                grid[a + 1][b + 2].unwrap(),
                grid[a][b + 1].unwrap(),
            ];
            mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap())
                .unwrap();
        }
    }
    mesh
}

/// The wetted-interface bookkeeping (FSI1's, plus vertex velocities).
pub struct Interface {
    pub wetted: Vec<NodeId>,
    pub reference: Vec<(f64, f64)>,
    /// Ordered boundary walk: indices into `wetted` (`usize::MAX` marks
    /// the fixed anchor vertices inside the cylinder / at the clamp).
    pub walk: Vec<(usize, (f64, f64))>,
    /// The three wetted edges as indices into `wetted`: bottom (root →
    /// tip), tip (bottom corner → top corner, corners included), top
    /// (tip → root) — the P5 patch generator's input.
    pub bottom: Vec<usize>,
    pub tip: Vec<usize>,
    pub top: Vec<usize>,
}

impl Interface {
    pub fn build(mesh: &Mesh) -> Self {
        let eps = 1e-9;
        let on_bottom = |p: Vector3<f64>| (p.y - FLAG_Y0).abs() < eps;
        let on_top = |p: Vector3<f64>| (p.y - FLAG_Y1).abs() < eps;
        let on_tip = |p: Vector3<f64>| (p.x - FLAG_X1).abs() < eps;
        let clamped = |p: Vector3<f64>| (p.x - FLAG_X0).abs() < eps;

        let mut wetted: Vec<(NodeId, (f64, f64))> = mesh
            .nodes
            .iter()
            .filter(|(_, node)| {
                let p = node.position();
                (on_bottom(p) || on_top(p) || on_tip(p)) && !clamped(p)
            })
            .map(|(&id, node)| (id, (node.position().x, node.position().y)))
            .collect();
        wetted.sort_by_key(|(id, _)| *id);
        let index_of = |id: NodeId| wetted.iter().position(|(w, _)| *w == id).unwrap();

        let mut bottom: Vec<(NodeId, f64)> = mesh
            .nodes
            .iter()
            .filter(|(_, n)| on_bottom(n.position()) && !clamped(n.position()))
            .map(|(&id, n)| (id, n.position().x))
            .collect();
        bottom.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
        let mut tip: Vec<(NodeId, f64)> = mesh
            .nodes
            .iter()
            .filter(|(_, n)| {
                let p = n.position();
                on_tip(p) && !on_bottom(p) && !on_top(p)
            })
            .map(|(&id, n)| (id, n.position().y))
            .collect();
        tip.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());
        let mut top: Vec<(NodeId, f64)> = mesh
            .nodes
            .iter()
            .filter(|(_, n)| on_top(n.position()) && !clamped(n.position()))
            .map(|(&id, n)| (id, n.position().x))
            .collect();
        top.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap());

        let mut walk: Vec<(usize, (f64, f64))> = Vec::new();
        walk.push((usize::MAX, (0.22, FLAG_Y0)));
        walk.push((usize::MAX, (FLAG_X0, FLAG_Y0)));
        for (id, _) in &bottom {
            walk.push((index_of(*id), (0.0, 0.0)));
        }
        for (id, _) in &tip {
            walk.push((index_of(*id), (0.0, 0.0)));
        }
        for (id, _) in &top {
            walk.push((index_of(*id), (0.0, 0.0)));
        }
        walk.push((usize::MAX, (FLAG_X0, FLAG_Y1)));
        walk.push((usize::MAX, (0.22, FLAG_Y1)));

        let reference = wetted.iter().map(|(_, p)| *p).collect();
        let bottom_idx: Vec<usize> = bottom.iter().map(|(id, _)| index_of(*id)).collect();
        // The tip edge with its corners: the last bottom node, the tip's
        // interior nodes (sorted by y), the first top node.
        let mut tip_idx: Vec<usize> = vec![*bottom_idx.last().expect("bottom edge")];
        tip_idx.extend(tip.iter().map(|(id, _)| index_of(*id)));
        let top_idx: Vec<usize> = top.iter().map(|(id, _)| index_of(*id)).collect();
        tip_idx.push(*top_idx.first().expect("top edge"));
        Self {
            wetted: wetted.into_iter().map(|(id, _)| id).collect(),
            reference,
            walk,
            bottom: bottom_idx,
            tip: tip_idx,
            top: top_idx,
        }
    }

    /// Deformed polygon vertices for the interface vector `d`.
    pub fn polygon(&self, d: &[f64]) -> Vec<(f64, f64)> {
        self.walk
            .iter()
            .map(|&(k, anchor)| {
                if k == usize::MAX {
                    anchor
                } else {
                    let (x0, y0) = self.reference[k];
                    (x0 + d[2 * k], y0 + d[2 * k + 1])
                }
            })
            .collect()
    }

    /// Per-vertex velocities for the interface velocity vector `ddot`
    /// (anchors do not move).
    pub fn walk_velocities(&self, ddot: &[f64]) -> Vec<(f64, f64)> {
        self.walk
            .iter()
            .map(|&(k, _)| {
                if k == usize::MAX {
                    (0.0, 0.0)
                } else {
                    (ddot[2 * k], ddot[2 * k + 1])
                }
            })
            .collect()
    }

    /// The deformed wetted edges for the P5 patch generator, each edge
    /// starting / ending at the clamp line's fixed corner.
    pub fn edges(&self, d: &[f64]) -> rtx_cfd::mesh::patch_gen::FlagEdges {
        let at = |k: usize| -> [f64; 2] {
            let (x0, y0) = self.reference[k];
            [x0 + d[2 * k], y0 + d[2 * k + 1]]
        };
        let mut bottom = vec![[FLAG_X0, FLAG_Y0]];
        bottom.extend(self.bottom.iter().map(|&k| at(k)));
        let tip: Vec<[f64; 2]> = self.tip.iter().map(|&k| at(k)).collect();
        let mut top: Vec<[f64; 2]> = self.top.iter().map(|&k| at(k)).collect();
        top.push([FLAG_X0, FLAG_Y1]);
        rtx_cfd::mesh::patch_gen::FlagEdges { bottom, tip, top }
    }

    /// Deformed wetted node positions for the transfer.
    pub fn deformed_nodes(&self, d: &[f64]) -> Vec<Vector3<f64>> {
        self.reference
            .iter()
            .enumerate()
            .map(|(k, &(x0, y0))| Vector3::new(x0 + d[2 * k], y0 + d[2 * k + 1], 0.0))
            .collect()
    }
}

pub fn clamp_left(mesh: &Mesh) -> BoundaryConditionSet {
    let clamped: Vec<NodeId> = mesh
        .nodes
        .iter()
        .filter(|(_, node)| (node.position().x - FLAG_X0).abs() < 1e-9)
        .map(|(&id, _)| id)
        .collect();
    let mut set = BoundaryConditionSet::new();
    for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] {
        set.add_condition(BoundaryCondition::Dirichlet(DirichletBC {
            nodes: clamped.clone(),
            components: vec![component],
            condition_type: DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))),
            time_range: None,
            ramping_factor: 1.0,
            gradual_enforcement: false,
        }));
    }
    set
}

pub fn mid_amp(series: &[f64]) -> (f64, f64) {
    let max = series.iter().copied().fold(f64::MIN, f64::max);
    let min = series.iter().copied().fold(f64::MAX, f64::min);
    (0.5 * (max + min), 0.5 * (max - min))
}

/// Median of a sample buffer (sorts in place; empty buffers read 0).
pub fn median(samples: &mut [f64]) -> f64 {
    if samples.is_empty() {
        return 0.0;
    }
    samples.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let n = samples.len();
    if n % 2 == 1 {
        samples[n / 2]
    } else {
        0.5 * (samples[n / 2 - 1] + samples[n / 2])
    }
}

/// Frequency from linearly interpolated upward crossings of the mean.
pub fn crossing_frequency(times: &[f64], series: &[f64]) -> Option<f64> {
    let (mean, _) = mid_amp(series);
    let mut crossings: Vec<f64> = Vec::new();
    for k in 1..series.len() {
        let (a, b) = (series[k - 1] - mean, series[k] - mean);
        if a < 0.0 && b >= 0.0 {
            crossings.push(times[k - 1] + (a / (a - b)) * (times[k] - times[k - 1]));
        }
    }
    if crossings.len() < 3 {
        return None;
    }
    Some((crossings.len() - 1) as f64 / (crossings.last().unwrap() - crossings.first().unwrap()))
}

/// The fluid + interface machinery every FSI2 test shares: the embedded
/// solver configured for the benchmark channel, the deformable-geometry
/// lock, and the load sampling with its spike clamp and (optional)
/// surface smoothing.
pub struct Fsi2Harness {
    pub case: BenchmarkCase,
    pub mesh: Mesh,
    pub interface: Interface,
    pub a_node: NodeId,
    pub ny: usize,
    pub nx: usize,
    pub h: f64,
    pub mu: f64,
    pub dt_fluid: f64,
    /// Traction smoothing radius along the surface, in metres
    /// (0 disables). Set from `RTX_FSI2_SMOOTH` (in multiples of `h`).
    pub smooth_radius: f64,
    /// The deformable geometry AND its velocity, behind one lock: the
    /// fluid's per-step mask rebuild reads the polygon; the no-slip
    /// closure reads both.
    pub shared: Arc<RwLock<(PolygonSdf, Vec<(f64, f64)>)>>,
    pub spiked_total: Cell<usize>,
}

impl Fsi2Harness {
    /// Build the harness plus the configured solver and an at-rest field.
    pub fn build(
        ny: usize,
        flag_nx: usize,
        smooth_in_h: f64,
    ) -> (Self, EmbeddedPisoSolver, FlowField) {
        Self::build_case(FSI2, ny, flag_nx, smooth_in_h)
    }

    /// Build the harness for a benchmark case (FSI2 or FSI3 parameters
    /// on the shared geometry).
    pub fn build_case(
        case: BenchmarkCase,
        ny: usize,
        flag_nx: usize,
        smooth_in_h: f64,
    ) -> (Self, EmbeddedPisoSolver, FlowField) {
        let h = H / ny as f64;
        let nx = (L / h).round() as usize;
        let mu = RHO_F * NU_F;
        let u_mean = case.u_mean;
        let u_peak = 1.5 * 1.5 * u_mean;
        // The fluid's explicit limit; the coupling (and the flag's
        // Newmark) run at `subcycle` fluid steps per coupled step.
        let dt_fluid = 0.25 / (2.0 * u_peak / h + 4.0 * NU_F / (h * h));

        let mesh = flag_mesh(flag_nx, 2);
        let interface = Interface::build(&mesh);
        let a_node = mesh
            .nodes
            .iter()
            .find(|(_, n)| {
                (n.position().x - 0.6).abs() < 1e-9 && (n.position().y - 0.2).abs() < 1e-9
            })
            .map(|(&id, _)| id)
            .expect("point A");

        let zero_d = vec![0.0; 2 * interface.wetted.len()];
        // The polygon lives behind the lock as an INDEXED SDF
        // (bit-identical query; the brute-force walk was measured at
        // 51% of the fluid step, called for every mask cell and ghost).
        let shared = Arc::new(RwLock::new((
            PolygonSdf::new(interface.polygon(&zero_d)),
            interface.walk_velocities(&zero_d),
        )));
        let sdf_shared = shared.clone();
        let vel_shared = shared.clone();

        let config = CfdConfig::new()
            .with_density(RHO_F)
            .with_viscosity(mu)
            .with_reference_velocity(u_mean)
            .with_reference_length(0.1);
        let params = EmbeddedParameters {
            corrector_steps: 2,
            tolerance: 1e-7,
            boundaries: AleBoundaries {
                left: SideBoundary::Velocity,
                right: SideBoundary::PressureOutlet,
                bottom: SideBoundary::Velocity,
                top: SideBoundary::Velocity,
            },
            poisson_solver: PoissonSolverKind::Multigrid,
            poisson_precision: rtx_cfd::solvers::incompressible::MgPrecision::F64,
            poisson_smoother: rtx_cfd::solvers::incompressible::MgSmoother::Lexicographic,
            // TVD, deliberately: FSI2 marches in time and needs the
            // shedding physics upwind's numerical viscosity killed on
            // these grids (CFD3's finding). The limiter chatter that
            // defeats steady fixed points (FSI1's finding) is harmless
            // here — each step's fixed point is the interface
            // displacement of THAT step, not a steady load.
            convection_scheme: ConvectionScheme::TvdVanAlbada,
        };
        let mut solver = EmbeddedPisoSolver::new(config, params).unwrap();
        solver.set_boundary_velocity(move |x, y, t| {
            if x <= 0.0 {
                (inflow_for(u_mean, y, t), 0.0)
            } else {
                (0.0, 0.0)
            }
        });
        solver.set_moving_body(
            EmbeddedBody::from_sdf(move |x, y, _| {
                let geometry = sdf_shared.read().unwrap();
                circle_sdf(x, y).min(geometry.0.signed_distance(x, y))
            })
            .with_surface_velocity(move |x, y, _| {
                let geometry = vel_shared.read().unwrap();
                if circle_sdf(x, y) <= geometry.0.signed_distance(x, y) {
                    (0.0, 0.0)
                } else {
                    polygon_interface_velocity(geometry.0.vertices(), &geometry.1, x, y)
                }
            }),
        );

        // Start at rest; the ramp brings the inflow up from zero.
        let mut field = FlowField::new(nx, ny, h, h).unwrap();
        solver.initialize(&mut field).unwrap();

        let harness = Self {
            case,
            mesh,
            interface,
            a_node,
            ny,
            nx,
            h,
            mu,
            dt_fluid,
            smooth_radius: smooth_in_h * h,
            shared,
            spiked_total: Cell::new(0),
        };
        (harness, solver, field)
    }

    /// Publish an interface geometry (+ velocity) to the fluid.
    pub fn set_geometry(&self, d: &[f64], ddot: &[f64]) {
        let mut geometry = self.shared.write().unwrap();
        geometry.0 = PolygonSdf::new(self.interface.polygon(d));
        geometry.1 = self.interface.walk_velocities(ddot);
    }

    /// Surface drag and lift on cylinder + flag at the current geometry.
    pub fn measure_force(&self, solver: &EmbeddedPisoSolver, field: &FlowField) -> (f64, f64) {
        let mask = solver.mask().unwrap();
        let body = solver.body().unwrap();
        let vertices = self.shared.read().unwrap().0.vertices().to_vec();
        // Collect, then clamp, then integrate: the same 20x-median spike
        // clamp the coupling loads carry. Without it the REPORTED
        // drag/lift at large deformation are dominated by the rare wild
        // reconstructions (the s = 1 benchmark run printed +-4,000-scale
        // load swings against a +-78 reference while its displacements
        // matched the benchmark to 0.1%).
        let mut samples: Vec<(f64, f64, f64)> = Vec::new();
        let poly_probe = EmbeddedBody::polygon(vertices.clone());
        for s in poly_probe.surface_samples(0.5 * self.h) {
            if circle_sdf(s.x, s.y) < 1e-9 {
                continue;
            }
            if let Some((tx, ty)) = mask.traction_at(
                body, &field.u, &field.v, &field.p, self.mu, 0.0, s.x, s.y, s.nx, s.ny,
            ) {
                samples.push((tx, ty, s.ds));
            }
        }
        let circle_probe = EmbeddedBody::circle(0.2, 0.2, 0.05);
        for s in circle_probe.surface_samples(0.5 * self.h) {
            if polygon_signed_distance(&vertices, s.x, s.y) < 1e-9 {
                continue;
            }
            if let Some((tx, ty)) = mask.traction_at(
                body, &field.u, &field.v, &field.p, self.mu, 0.0, s.x, s.y, s.nx, s.ny,
            ) {
                samples.push((tx, ty, s.ds));
            }
        }
        let mut magnitudes: Vec<f64> = samples
            .iter()
            .map(|(tx, ty, _)| (tx * tx + ty * ty).sqrt())
            .collect();
        magnitudes.sort_by(|a, b| a.partial_cmp(b).unwrap());
        let median = magnitudes.get(magnitudes.len() / 2).copied().unwrap_or(0.0);
        let cap = 20.0 * median;
        let mut drag = 0.0;
        let mut lift = 0.0;
        for (tx, ty, ds) in samples {
            let norm = (tx * tx + ty * ty).sqrt();
            let scale = if median > 0.0 && norm > cap {
                cap / norm
            } else {
                1.0
            };
            drag += tx * scale * ds;
            lift += ty * scale * ds;
        }
        (drag, lift)
    }

    /// Tractions on the flag's wetted surface for a given geometry, from
    /// the solver's current field/mask; returns the transferred nodal
    /// forces, the conservation defect, and the samples dropped (probe
    /// failures plus spike rejections).
    pub fn sample_load(
        &self,
        solver: &EmbeddedPisoSolver,
        field: &FlowField,
        d: &[f64],
    ) -> (Vec<(NodeId, Vector3<f64>)>, f64, usize) {
        let vertices = self.interface.polygon(d);
        let poly_probe = EmbeddedBody::polygon(vertices);
        let mask = solver.mask().unwrap();
        let body = solver.body().unwrap();
        let mut faces = Vec::new();
        let mut tractions: Vec<Vector3<f64>> = Vec::new();
        let mut skipped = 0usize;
        for s in poly_probe.surface_samples(0.5 * self.h) {
            if circle_sdf(s.x, s.y) < 1e-9 {
                continue; // buried in the cylinder
            }
            match mask.traction_at(
                body, &field.u, &field.v, &field.p, self.mu, 0.0, s.x, s.y, s.nx, s.ny,
            ) {
                Some((tx, ty)) => {
                    faces.push(FluidFace {
                        centroid: Vector3::new(s.x, s.y, 0.0),
                        normal: Vector3::new(s.nx, s.ny, 0.0),
                        area: s.ds,
                    });
                    tractions.push(Vector3::new(tx, ty, 0.0));
                }
                None => skipped += 1,
            }
        }
        // Spike guard: a near-degenerate reconstruction can return a
        // finite but wild traction (the linear-fit condition sits just
        // above its truncation threshold at concave junctions). CLAMP
        // samples to 20x the median magnitude, keeping their direction —
        // the physical load varies smoothly along the surface — and COUNT
        // them: a non-zero count is a measurement of the pathology, not a
        // silent repair. Clamping, not dropping: a hard drop threshold
        // makes the coupling pass discontinuous in the candidate geometry
        // (a boundary sample flips in/out of the kept set between
        // subiterations, and the load jumps by the spike magnitude —
        // measured as a residual bouncing at the scale of the step
        // increment); the clamp is continuous.
        let mut magnitudes: Vec<f64> = tractions.iter().map(nalgebra::Vector3::norm).collect();
        magnitudes.sort_by(|a, b| a.partial_cmp(b).unwrap());
        let median = magnitudes.get(magnitudes.len() / 2).copied().unwrap_or(0.0);
        if median > 0.0 {
            let cap = 20.0 * median;
            for traction in &mut tractions {
                let norm = traction.norm();
                if norm > cap {
                    *traction *= cap / norm;
                    self.spiked_total.set(self.spiked_total.get() + 1);
                }
            }
        }
        // Surface smoothing (after the clamp: the clamp kills the wild
        // outliers, the smoothing spreads what remains over the stencil
        // the cell resolution can actually support — this is the
        // interface-noise-floor lever, measured by
        // `fsi2_interface_noise.rs`).
        if self.smooth_radius > 0.0 {
            tractions = smooth_tractions(&faces, &tractions, self.smooth_radius).unwrap();
        }
        let nodes_now = self.interface.deformed_nodes(d);
        let surface = WettedSurface::build(&faces, &nodes_now).expect("transfer build");
        let nodal = surface.transfer_load(&faces, &tractions).unwrap();
        let total_sampled: Vector3<f64> =
            faces.iter().zip(&tractions).map(|(f, t)| t * f.area).sum();
        let total_nodal: Vector3<f64> = nodal.iter().sum();
        let conservation = (total_nodal - total_sampled).norm() / total_sampled.norm().max(1e-30);
        (
            self.interface
                .wetted
                .iter()
                .zip(nodal)
                .map(|(&id, f)| (id, f))
                .collect(),
            conservation,
            skipped,
        )
    }

    /// Run `subcycle` fluid substeps from the current solver/field state,
    /// interpolating the interface geometry from `d_n` to `d_candidate`
    /// across the substeps with the candidate's constant interface
    /// velocity `(d_candidate - d_n) / dt` — one coupling pass's fluid
    /// half, exactly as the coupled march runs it.
    pub fn advance_subcycled(
        &self,
        solver: &mut EmbeddedPisoSolver,
        field: &mut FlowField,
        d_n: &[f64],
        d_candidate: &[f64],
        subcycle: usize,
        v_n: Option<&[f64]>,
    ) {
        self.advance_subcycled_with(
            solver,
            field,
            d_n,
            d_candidate,
            subcycle,
            v_n,
            self.dt_fluid,
        );
    }

    /// [`Self::advance_subcycled`] at an explicit fluid substep `dt_fluid`
    /// (the coupled step is `dt_fluid * subcycle`). The march's own path
    /// passes `self.dt_fluid` — same arithmetic, digit for digit; the
    /// coupling-level rescue passes `dt_fluid / n` for its substeps.
    #[allow(clippy::too_many_arguments)]
    pub fn advance_subcycled_with(
        &self,
        solver: &mut EmbeddedPisoSolver,
        field: &mut FlowField,
        d_n: &[f64],
        d_candidate: &[f64],
        subcycle: usize,
        v_n: Option<&[f64]>,
        dt_fluid: f64,
    ) {
        let dt = dt_fluid * subcycle as f64;
        let mean_velocity: Vec<f64> = d_candidate
            .iter()
            .zip(d_n)
            .map(|(new, old)| (new - old) / dt)
            .collect();
        for m in 1..=subcycle {
            let fraction = m as f64 / subcycle as f64;
            let (d_sub, ddot_sub): (Vec<f64>, Vec<f64>) = match v_n {
                // Constant velocity over the step: the geometry moves
                // linearly and the wall velocity JUMPS at the step
                // boundary — harmless for a heavy flag, but the
                // incompressible fluid answers a velocity jump with an
                // impulsive added-mass load ~ rho L dv / dt_fluid, which
                // at unit density ratio destroyed the flag in one step.
                None => (
                    d_n.iter()
                        .zip(d_candidate)
                        .map(|(old, new)| old + fraction * (new - old))
                        .collect(),
                    mean_velocity.clone(),
                ),
                // C^1 interface motion: constant acceleration across the
                // step from the previous end-of-step velocity to the
                // trapezoidal end velocity 2 dd/dt - v_n (Newmark
                // average acceleration's own kinematics), so the wall
                // velocity is continuous at the step boundary and the
                // impulse is gone. The end-of-substep velocity goes with
                // the end-of-substep geometry.
                Some(v_start) => {
                    let mut d_sub = Vec::with_capacity(d_n.len());
                    let mut ddot_sub = Vec::with_capacity(d_n.len());
                    for k in 0..d_n.len() {
                        let v_end = 2.0 * mean_velocity[k] - v_start[k];
                        let accel = (v_end - v_start[k]) / dt;
                        let tau = fraction * dt;
                        d_sub.push(d_n[k] + v_start[k] * tau + 0.5 * accel * tau * tau);
                        ddot_sub.push(v_start[k] + accel * tau);
                    }
                    (d_sub, ddot_sub)
                }
            };
            self.set_geometry(&d_sub, &ddot_sub);
            futures::executor::block_on(solver.advance(field, dt_fluid)).unwrap();
        }
    }
}