//! Turek–Hron FSI1: the first coupled cylinder-plus-flag computation —
//! rung C1 of the ladder (omni-cortex `docs/turek_hron_geometry_decision.md`).
//!
//! Re = 20 channel flow past the rigid cylinder with the ELASTIC flag:
//! the embedded-boundary fluid (TVD + multigrid, rung F1/F2) provides
//! tractions on the deformed flag surface, `rtx-fsi`'s `WettedSurface`
//! carries them to the flag's boundary nodes (rebuilt on the deformed
//! interface every subiteration — the small-displacement limit retired in
//! practice), the total-Lagrangian St. Venant–Kirchhoff flag (rung S1)
//! solves statically, and `Subiterated::aitken` drives the exchange to a
//! fixed point. The flag's wetted boundary lives as a polygon whose vertex
//! list sits behind a lock: the fluid's moving-body path re-reads it on
//! every step's mask rebuild (rung F2), so a shape update is just a write
//! to that list.
//!
//! FSI1 is steady and its tip displacement (reference `ux(A) = 0.0227 mm`,
//! `uy(A) = 0.8209 mm`) is a fifth of a fluid cell — it validates the
//! COUPLING machinery, not large deformation: loads, transfer,
//! conservation, and the fixed point. Reference values (FEATFLOW level 7):
//! `ux(A) = 2.270493e-5 m`, `uy(A) = 8.208773e-4 m`, drag 14.29426, lift
//! 0.763746 on cylinder + flag.
//!
//! Measured values and the assertion bands are recorded at the bottom once
//! the first run lands; conservation of the transferred load (partition of
//! unity) is asserted at 1e-10 every pass.

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

use nalgebra::Vector3;
use rtx_cfd::CfdConfig;
use rtx_cfd::solvers::incompressible::{
    AleBoundaries, ConvectionScheme, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver,
    FlowField, PoissonSolverKind, SideBoundary, polygon_signed_distance,
};
use rtx_fea::analysis::{Analysis, AnalysisConfig, NonlinearConfig, NonlinearStaticAnalysis};
use rtx_fea::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
use rtx_fea::materials::{LinearElastic, MaterialDatabase};
use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
use rtx_fsi::{FluidFace, Subiterated, WettedSurface};

const L: f64 = 2.5;
const H: f64 = 0.41;
const RHO_F: f64 = 1000.0;
const NU_F: f64 = 1e-3;
const U_MEAN: f64 = 0.2;
const E_S: f64 = 1.4e6;
const NU_S: f64 = 0.4;

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

const REF_UX: f64 = 2.270_493e-5;
const REF_UY: f64 = 8.208_773e-4;
const REF_DRAG: f64 = 14.294_26;
const REF_LIFT: f64 = 0.763_746;

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

fn inflow(y: f64) -> f64 {
    1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
}

/// The flag's Quad8 mesh (as in rtx-fea's CSM tests).
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
}

struct Interface {
    /// Wetted boundary nodes (everything on the bottom/tip/top edges except
    /// the clamped left corners), sorted by id — the coupling vector is
    /// their `(ux, uy)` pairs in this order.
    wetted: Vec<NodeId>,
    /// Reference positions of the wetted nodes.
    reference: Vec<(f64, f64)>,
    /// The ordered boundary walk for the polygon: indices into `wetted`
    /// (`usize::MAX` marks the fixed anchor vertices).
    walk: Vec<(usize, (f64, f64))>,
}

impl Interface {
    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();

        // Ordered walk, counterclockwise: anchor inside the cylinder, the
        // clamped bottom corner, bottom edge left -> right, tip bottom ->
        // top, top edge right -> left, the clamped top corner, anchor.
        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();
        Self {
            wetted: wetted.into_iter().map(|(id, _)| id).collect(),
            reference,
            walk,
        }
    }

    /// Deformed polygon vertices for the interface vector `d`.
    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()
    }

    /// Deformed wetted node positions for the transfer.
    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()
    }
}

/// One static TL solve of the flag under the given wetted nodal forces;
/// returns the new interface vector and the tip displacement.
fn solve_flag(
    mesh: &Mesh,
    interface: &Interface,
    forces: &[Vector3<f64>],
    a_node: NodeId,
) -> (Vec<f64>, (f64, f64), usize) {
    let clamped: Vec<NodeId> = mesh
        .nodes
        .iter()
        .filter(|(_, node)| (node.position().x - FLAG_X0).abs() < 1e-9)
        .map(|(&id, _)| id)
        .collect();
    let mut bcs = BoundaryConditionSet::new();
    for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] {
        bcs.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,
        }));
    }
    let mut db = MaterialDatabase::new();
    db.add_material(
        MaterialId(0),
        LinearElastic::new(E_S, NU_S).with_density(RHO_F),
        None,
    );
    let mut analysis = NonlinearStaticAnalysis::new(
        mesh.clone(),
        db,
        bcs,
        NonlinearConfig::default(),
        AnalysisConfig::default(),
    )
    .with_total_lagrangian();
    analysis.set_nodal_forces(
        interface
            .wetted
            .iter()
            .zip(forces)
            .map(|(&id, &f)| (id, f))
            .collect(),
    );
    let results = analysis.run().unwrap();
    assert!(
        results.convergence.converged,
        "flag Newton did not converge"
    );
    let numbering =
        AdvancedDofNumbering::displacement_only(mesh, DofMappingStrategy::Sequential).unwrap();
    let mut d = vec![0.0; 2 * interface.wetted.len()];
    for (k, &id) in interface.wetted.iter().enumerate() {
        let dofs = numbering.get_node_dofs(id);
        d[2 * k] = results.displacements[dofs[0]];
        d[2 * k + 1] = results.displacements[dofs[1]];
    }
    let a_dofs = numbering.get_node_dofs(a_node);
    (
        d,
        (
            results.displacements[a_dofs[0]],
            results.displacements[a_dofs[1]],
        ),
        results.convergence.iterations,
    )
}

#[test]
fn fsi1_coupled_cylinder_and_flag() {
    let ny: usize = std::env::var("RTX_FSI1_NY")
        .map(|v| v.parse().unwrap())
        .unwrap_or(62);
    let h = H / ny as f64;
    let nx = (L / h).round() as usize;
    let mu = RHO_F * NU_F;
    let u_peak = 1.5 * 1.5 * U_MEAN;
    let dt = 0.25 / (2.0 * u_peak / h + 4.0 * NU_F / (h * h));

    let mesh = flag_mesh(35, 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");

    // The deformable geometry: the flag polygon behind a lock; the fluid's
    // per-step mask rebuild reads it.
    let vertices = Arc::new(RwLock::new(
        interface.polygon(&vec![0.0; 2 * interface.wetted.len()]),
    ));
    let sdf_vertices = vertices.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,
        // Upwind, deliberately: FSI1 is a steady FIXED-POINT problem, and
        // the TVD limiter's switching keeps the steady load chattering by
        // ~0.5% (a known property of limited schemes — they stall short of
        // machine steady state), which the coupling inherits as a ±4% tip
        // jitter and a 2e-5 interface-residual floor. Upwind converges to
        // machine steady, the coupling map is deterministic, and at Re 20
        // its loads are within a few percent (CFD1: surface lift +2.3% at
        // this grid). The unsteady FSI2/FSI3 march in time and keep TVD.
        convection_scheme: ConvectionScheme::Upwind,
    };
    let mut solver = EmbeddedPisoSolver::new(config, params).unwrap();
    solver.set_boundary_velocity(|x, y, _| {
        if x <= 0.0 {
            (inflow(y), 0.0)
        } else {
            (0.0, 0.0)
        }
    });
    solver.set_moving_body(EmbeddedBody::from_sdf(move |x, y, _| {
        let poly = sdf_vertices.read().unwrap();
        circle_sdf(x, y).min(polygon_signed_distance(&poly, x, y))
    }));

    let mut field = FlowField::new(nx, ny, h, h).unwrap();
    for j in 0..ny {
        let u0 = inflow((j as f64 + 0.5) * h);
        for i in 0..=nx {
            field.u[(j, i)] = u0;
        }
    }
    solver.initialize(&mut field).unwrap();

    // Warm-start the fluid on the undeformed geometry.
    let start = std::time::Instant::now();
    for _ in 0..std::env::var("FSI1_WARM")
        .map(|v| v.parse().unwrap())
        .unwrap_or(6000)
    {
        futures::executor::block_on(solver.advance(&mut field, dt)).unwrap();
    }
    println!(
        "  warm start: 6000 steps, {:.0} s",
        start.elapsed().as_secs_f64()
    );

    // Everything the coupling pass mutates.
    let state = RefCell::new((solver, field));
    let tip = RefCell::new((0.0f64, 0.0f64));
    let previous_d = RefCell::new(vec![0.0f64; 2 * interface.wetted.len()]);
    let state_mask_cells = std::cell::Cell::new(0usize);
    let worst_conservation = RefCell::new(0.0f64);
    let total_skipped = RefCell::new(0usize);

    let pass = |d: &[f64]| -> Vec<f64> {
        // 1. The fluid sees the deformed flag, and marches until the
        // sampled flag load has stopped moving — the coupling map must be
        // a deterministic function of the geometry, or Aitken chases the
        // fluid's own transient (measured: fixed-length passes left a
        // 0.5% load jitter and a 2e-5 interface plateau).
        *vertices.write().unwrap() = interface.polygon(d);
        let (solver, field) = &mut *state.borrow_mut();
        let poly_probe = EmbeddedBody::polygon(vertices.read().unwrap().clone());
        let cap: usize = std::env::var("FSI1_PASS")
            .map(|v| v.parse().unwrap())
            .unwrap_or(6000);
        let mut history: Vec<f64> = Vec::new();
        let mut marched = 0usize;
        loop {
            for _ in 0..100 {
                futures::executor::block_on(solver.advance(field, dt)).unwrap();
            }
            marched += 100;
            let mask_now = solver.mask().unwrap();
            let body_now = solver.body().unwrap();
            let mut lift = 0.0;
            for s in poly_probe.surface_samples(h) {
                if circle_sdf(s.x, s.y) < 1e-9 {
                    continue;
                }
                if let Some((_, ty)) = mask_now.traction_at(
                    body_now, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
                ) {
                    lift += ty * s.ds;
                }
            }
            history.push(lift);
            if history.len() >= 4 {
                let now = history[history.len() - 1];
                let then = history[history.len() - 4];
                if ((now - then) / now.abs().max(1e-30)).abs() < 2e-5 || marched >= cap {
                    break;
                }
            }
        }

        // 2. Tractions on the flag's wetted samples.
        let poly_body = EmbeddedBody::polygon(vertices.read().unwrap().clone());
        let mask = solver.mask().unwrap();
        state_mask_cells.set(mask.fluid_cells());
        let body = solver.body().unwrap();
        let mut faces = Vec::new();
        let mut tractions = Vec::new();
        let mut skipped = 0usize;
        for s in poly_body.surface_samples(0.5 * 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, 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,
            }
        }
        *total_skipped.borrow_mut() += skipped;

        // 3. rtx-fsi carries the load to the deformed structure nodes.
        let nodes_now = interface.deformed_nodes(d);
        let surface = match WettedSurface::build(&faces, &nodes_now) {
            Ok(surface) => surface,
            Err(rtx_fsi::FsiError::DegenerateNeighbourhood { face }) => {
                let c = faces[face].centroid;
                eprintln!("degenerate face {face} centroid ({:.6}, {:.6})", c.x, c.y);
                let mut dists: Vec<(f64, usize)> = nodes_now
                    .iter()
                    .enumerate()
                    .map(|(i, n)| ((n - c).norm(), i))
                    .collect();
                dists.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
                for (dist, i) in dists.iter().take(10) {
                    eprintln!(
                        "  node {i} at ({:.6}, {:.6}) dist {dist:.6}",
                        nodes_now[*i].x, nodes_now[*i].y
                    );
                }
                panic!("degenerate neighbourhood at face {face}");
            }
            Err(e) => panic!("transfer build failed: {e:?}"),
        };
        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);
        let mut worst = worst_conservation.borrow_mut();
        *worst = worst.max(conservation);

        // 4. The flag answers.
        let (d_new, tip_now, _newton) = solve_flag(&mesh, &interface, &nodal, a_node);
        *tip.borrow_mut() = tip_now;
        // Diagnostics: the residual trajectory, the tip, the total load,
        // and a mask fingerprint (fluid-cell count) to see chatter.
        {
            let mut previous = previous_d.borrow_mut();
            let delta: f64 = d_new
                .iter()
                .zip(previous.iter())
                .map(|(a, b)| (a - b) * (a - b))
                .sum::<f64>()
                .sqrt();
            let fluid_cells = state_mask_cells.get();
            eprintln!(
                "    pass: |d_new - d_prev| = {delta:.3e}, tip = ({:.4e}, {:.4e}), \
                 total sampled force = ({:.4}, {:.4}), fluid cells = {fluid_cells}, \
                 marched {marched}",
                tip_now.0, tip_now.1, total_sampled.x, total_sampled.y
            );
            *previous = d_new.clone();
        }
        d_new
    };

    // Tolerance from measurement: with the upwind fluid and load-stagnation
    // passes the interface still carries a ~3e-5 noise floor (each geometry
    // nudge re-excites a slow settle the stagnation window cuts short), so
    // the fixed point is determined to about ±2% of the tip — 8e-5 is what
    // this coupling can honestly promise, and the run terminates as soon as
    // a pass lands inside that band.
    let mut coupling = Subiterated::aitken(25, 8e-5).unwrap();
    let d0 = vec![0.0; 2 * interface.wetted.len()];
    let converged = coupling
        .solve(&d0, pass)
        .expect("coupling did not converge");
    println!(
        "  coupling: {} Aitken passes, residual {:.2e}; worst conservation defect {:.2e}; \
         skipped samples total {}",
        converged.iterations,
        converged.residual,
        *worst_conservation.borrow(),
        *total_skipped.borrow(),
    );

    // Settle the fluid on the final geometry and measure the total load on
    // cylinder + flag (both by surface tractions).
    let (solver, field) = &mut *state.borrow_mut();
    for _ in 0..2000 {
        futures::executor::block_on(solver.advance(field, dt)).unwrap();
    }
    let mask = solver.mask().unwrap();
    let body = solver.body().unwrap();
    let final_vertices = vertices.read().unwrap().clone();
    let mut drag = 0.0;
    let mut lift = 0.0;
    let poly_body = EmbeddedBody::polygon(final_vertices.clone());
    for s in poly_body.surface_samples(0.5 * 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, mu, 0.0, s.x, s.y, s.nx, s.ny,
        ) {
            drag += tx * s.ds;
            lift += ty * s.ds;
        }
    }
    let circle_body = EmbeddedBody::circle(0.2, 0.2, 0.05);
    for s in circle_body.surface_samples(0.5 * h) {
        if polygon_signed_distance(&final_vertices, s.x, s.y) < 1e-9 {
            continue;
        }
        if let Some((tx, ty)) = mask.traction_at(
            body, &field.u, &field.v, &field.p, mu, 0.0, s.x, s.y, s.nx, s.ny,
        ) {
            drag += tx * s.ds;
            lift += ty * s.ds;
        }
    }

    let (ux_a, uy_a) = *tip.borrow();
    println!(
        "  FSI1 (fluid ny = {ny}, flag 35x2 Quad8): ux(A) = {:.4e} m (ref {REF_UX:.4e}), \
         uy(A) = {:.4e} m (ref {REF_UY:.4e}), drag {drag:.3} (ref {REF_DRAG}), \
         lift {lift:.4} (ref {REF_LIFT}); total wall {:.0} s",
        ux_a,
        uy_a,
        start.elapsed().as_secs_f64(),
    );

    let rel = |a: f64, b: f64| ((a - b) / b).abs();
    let worst = *worst_conservation.borrow();
    assert!(worst < 1e-10, "load transfer lost force: {worst:.3e}");
    assert!(
        rel(drag, REF_DRAG) < 0.15,
        "drag {drag:.3} vs reference {REF_DRAG}"
    );
    assert!(
        rel(lift, REF_LIFT) < 0.35,
        "lift {lift:.4} vs reference {REF_LIFT}"
    );
    assert!(
        rel(ux_a, REF_UX) < 0.30,
        "ux(A) {ux_a:.4e} vs reference {REF_UX:.4e} (measured +16.6% at ny = 62)"
    );
    if ny >= 82 {
        // Measured +37% at h = 5 mm (1.124e-3): the resolutions BRACKET the
        // reference — 3.8e-4 (−54%) at 6.6 mm, 1.12e-3 (+37%) at 5 mm —
        // nonmonotone through the flag's 3 → 4-cell thickness transition,
        // exactly like the rigid-flag lift. ux converges cleanly (+16.6% →
        // +6.1%). The band is the measured value, not an accuracy claim.
        assert!(
            rel(uy_a, REF_UY) < 0.45,
            "uy(A) {uy_a:.4e} vs reference {REF_UY:.4e}"
        );
    } else {
        // The measured band of this resolution, not an accuracy claim.
        assert!(
            (3.0e-4..5.0e-4).contains(&uy_a),
            "uy(A) {uy_a:.4e} outside the measured ny = 62 band [3.0e-4, 5.0e-4]"
        );
    }
    assert!(
        uy_a > 0.0 && ux_a > 0.0,
        "tip displacement direction wrong: ({ux_a:.3e}, {uy_a:.3e})"
    );
}