rtx-cfd + rtx-fea: embedded-boundary PISO and total-Lagrangian SVK — the first two Turek–Hron rungs

The Turek–Hron geometry decision (omni-cortex
docs/turek_hron_geometry_decision.md) chose an embedded boundary on the
fixed Cartesian MAC grid over body-fitted unstructured ALE; this commit
builds the first rung on each side of the ladder, verified MMS-first.

rtx-cfd — solvers::incompressible::{embedded, embedded_body}:
EmbeddedPisoSolver is the fixed-grid PISO predictor/projection with
per-side domain boundaries (ALE's SideBoundary semantics, so the channel
has an outlet), a (x, y, t) boundary-velocity function, and an optional
EmbeddedBody (signed distance + surface velocity; circle / rectangle /
union). EmbeddedMask classifies cells (fluid iff phi > 0 at the centre)
and faces (fluid iff both cells fluid; ghost within 1.5 h; solid deeper);
the predictor updates fluid faces only, the projection enforces continuity
on fluid cells with zero coefficient across prescribed faces, ghost faces
are re-imposed after each projection from a boundary-intercept
least-squares linear fit (exact for linear fields), the net ghost mass flux
is removed uniformly so a Neumann projection stays compatible, and loads
come by two routes: surface-stress reconstruction (full viscous traction)
and a control-volume momentum balance.

Verified (tests/embedded_mms.rs, tests/turek_hron_cfd.rs):
- no body, closed box: bit-identical to PisoSolver over 200 steps;
- embedded off-centre circle MMS 16/32/64: velocity orders 0.92, 0.97
  (plain PISO 0.85, 0.91), pressure 0.96, 0.90, max |div u| <= 9e-8 on
  every fluid cell, compatibility correction 6e-4 -> 3e-5; force on the
  circle vs the exact surface integral: surface route 0.52 -> 0.29 -> 0.15,
  control-volume route 0.61 -> 0.30 -> 0.15 (both first order, two
  unrelated readings of the same solution);
- Turek–Hron CFD1 (Re 20, h = 10 mm, flag two cells thick), settled to
  four digits: surface drag 15.71 / lift 0.94, control-volume drag 15.62 /
  lift 1.08 vs reference 14.29 / 1.119 — the drag routes agree to 0.6%,
  both +9.5%. A coarse first number; the refinement study waits on a
  multigrid projection (SOR: 0.1 s/step at 250x41 in the test profile).

Fourteenth defect of the campaign: the fixed-grid PISO predictor zeroes
the transverse convective face velocity on its domain sides (exact for
walls); carried into a solver with an outlet it dropped the OUTGOING
momentum flux through the outlet side of the v control volumes, the last
column accumulated, and CFD1 went NaN at t ~ 4 s. Found by printing where
max |u| lived (x = 2.5) after halving dt changed nothing. Fluxes now come
from the stored boundary faces on every side.

rtx-fea — elements::total_lagrangian + NonlinearStaticAnalysis::
with_total_lagrangian(): Green–Lagrange strain, second Piola–Kirchhoff
stress from a St. Venant–Kirchhoff law on the material's Lamé parameters
(plane strain in 2-D), B_L of the current deformation, material plus
geometric tangent; dead-load body force per reference volume.

Verified (tests/total_lagrangian_svk.rs):
- zero displacement: the plane-strain stiffness to 1e-13;
- tangent = d f_int/du by central differences at 20% random displacement
  (Quad4, Quad8, Hex8): relative < 1e-7, symmetric to 1e-12;
- a 34-degree rigid rotation produces no internal force; the small-strain
  routine does (negative control);
- manufactured finite-strain solution, body force by FD of the exact
  P = F S: Quad4 orders 1.95, 1.98; Quad8 2.93, 3.03, 3.02 (an 8%
  amplitude, Green–Lagrange strain to -0.25 near SVK's compressive limit
  E = -1/3, broke Newton on fine meshes — the material, not the code; 3%
  is clean);
- Turek–Hron CSM1 at 70x4 Quad8: u(A) = (-7.060, -65.43) mm vs
  (-7.188, -66.10), 1.0% / 1.8%, converging from below (35x2: -65.14);
  CSM2: (-0.4604, -16.79) vs (-0.4690, -16.97), 1.1% / 1.8%.

rtx-cfd 293 -> 301 green (5 unit + 3 integration), rtx-fea 559 -> 564.

Co-Authored-By: Claude Fable 5 <[email protected]>

crates/specialized/rtx-cfd/src/solvers/incompressible/embedded.rs

@@ -0,0 +1,700 @@
+//! PISO on the fixed uniform staggered grid with an embedded body.
+//!
+//! This is the fixed-grid PISO scheme (`piso.rs`: explicit momentum
+//! predictor, SOR pressure-correction projections) with three additions:
+//!
+//! - **per-side domain boundaries** ([`SideBoundary`]: prescribed velocity,
+//!   slip wall, pressure outlet), with exactly the ALE solver's semantics
+//!   — the channel of the Turek–Hron benchmark needs an outlet and the
+//!   fixed-grid PISO had none;
+//! - a **time-dependent boundary-velocity function** `(x, y, t)` supplying
+//!   the normal data and the tangential wall values, and a source
+//!   `(x, y, t)`, as on the ALE solver;
+//! - an optional **embedded body** ([`EmbeddedBody`]) classified onto the
+//!   grid by [`EmbeddedMask`]: the momentum predictor updates only fluid
+//!   faces, the projection enforces continuity only on fluid cells with
+//!   zero coefficient across every prescribed face, and after each
+//!   projection the ghost faces are re-imposed from the corrected fluid
+//!   field (see `embedded_body.rs` for the reconstruction and the
+//!   compatibility correction).
+//!
+//! With no body, velocity on every side and no normal flow through the
+//! sides (a closed box), `advance` is the fixed-grid PISO step to the last
+//! bit — `tests/embedded_mms.rs` pins that degeneracy before it measures
+//! anything else. With through-flow the two differ by design: the
+//! fixed-grid PISO zeroes the transverse convective face velocity on its
+//! domain sides (exact for walls), this solver takes it from the stored
+//! boundary faces, which is what an inlet or outlet needs — dropping the
+//! outgoing flux at an outlet let the last column accumulate momentum and
+//! the Turek–Hron channel blew up at t ≈ 4 s.
+//!
+//! # Boundary history
+//!
+//! As the ALE solver learned (the twelfth defect of the campaign), the
+//! start-of-step boundary faces are *not* re-stamped from the boundary
+//! function: the previous step's end-of-step application is the material
+//! history the explicit predictor differentiates. Only the first step
+//! stamps `t = 0` data, via [`EmbeddedPisoSolver::initialize`] or lazily.
+
+use super::ale::{AleBoundaries, SideBoundary};
+use super::embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind};
+use super::{FlowField, SolverResult};
+use crate::{CfdConfig, CfdError, CfdResult};
+
+type VelocityFn = Box<dyn Fn(f64, f64, f64) -> (f64, f64) + Send + Sync>;
+type SourceFn = Box<dyn Fn(f64, f64, f64) -> (f64, f64) + Send + Sync>;
+
+/// Parameters for the embedded-boundary PISO solver.
+#[derive(Debug, Clone)]
+pub struct EmbeddedParameters {
+    /// Projection passes per step.
+    pub corrector_steps: usize,
+    /// Convergence tolerance on the normalised mass imbalance after
+    /// correction.
+    pub tolerance: f64,
+    /// Boundary type per domain side (all prescribed velocity by default).
+    /// The struct is the ALE solver's; the semantics are identical.
+    pub boundaries: AleBoundaries,
+}
+
+impl Default for EmbeddedParameters {
+    fn default() -> Self {
+        Self {
+            corrector_steps: 2,
+            tolerance: 1e-6,
+            boundaries: AleBoundaries::default(),
+        }
+    }
+}
+
+/// Result of one embedded PISO step.
+#[derive(Debug, Clone)]
+pub struct EmbeddedResult {
+    /// Base solver result information.
+    pub solver_result: SolverResult,
+    /// Number of projection passes performed.
+    pub corrector_steps_performed: usize,
+    /// The per-face compatibility correction applied to the ghost faces at
+    /// the end of the step (velocity units); zero without a body.
+    pub ghost_correction: f64,
+}
+
+/// The embedded-boundary PISO solver. See the module docs.
+pub struct EmbeddedPisoSolver {
+    config: CfdConfig,
+    parameters: EmbeddedParameters,
+    momentum_source: Option<SourceFn>,
+    boundary_velocity: Option<VelocityFn>,
+    body: Option<EmbeddedBody>,
+    mask: Option<EmbeddedMask>,
+    time: f64,
+    initialized: bool,
+}
+
+impl EmbeddedPisoSolver {
+    /// Create the solver.
+    pub fn new(config: CfdConfig, parameters: EmbeddedParameters) -> CfdResult<Self> {
+        config.validate()?;
+        Ok(Self {
+            config,
+            parameters,
+            momentum_source: None,
+            boundary_velocity: None,
+            body: None,
+            mask: None,
+            time: 0.0,
+            initialized: false,
+        })
+    }
+
+    /// Volumetric momentum source `(x, y, t) -> (f_x, f_y)` per unit volume.
+    pub fn set_momentum_source<F>(&mut self, f: F)
+    where
+        F: Fn(f64, f64, f64) -> (f64, f64) + Send + Sync + 'static,
+    {
+        self.momentum_source = Some(Box::new(f));
+    }
+
+    /// Boundary velocity `(x, y, t) -> (u, v)` on the domain sides: normal
+    /// component prescribed on `Velocity` and `SlipWall` sides, tangential
+    /// component the no-slip value on `Velocity` sides.
+    pub fn set_boundary_velocity<F>(&mut self, f: F)
+    where
+        F: Fn(f64, f64, f64) -> (f64, f64) + Send + Sync + 'static,
+    {
+        self.boundary_velocity = Some(Box::new(f));
+    }
+
+    /// Embed a body. The mask is built on the first step (the body is
+    /// treated as fixed in shape and position for now — moving bodies
+    /// arrive with the next rung).
+    pub fn set_body(&mut self, body: EmbeddedBody) {
+        self.body = Some(body);
+        self.mask = None;
+    }
+
+    /// The body, if any.
+    pub fn body(&self) -> Option<&EmbeddedBody> {
+        self.body.as_ref()
+    }
+
+    /// The mask, once built (after `initialize` or the first step).
+    pub fn mask(&self) -> Option<&EmbeddedMask> {
+        self.mask.as_ref()
+    }
+
+    /// Accumulated time.
+    pub fn time(&self) -> f64 {
+        self.time
+    }
+
+    /// Reset the accumulated time.
+    pub fn set_time(&mut self, t: f64) {
+        self.time = t;
+    }
+
+    /// Solver configuration.
+    pub fn config(&self) -> &CfdConfig {
+        &self.config
+    }
+
+    /// Solver parameters.
+    pub fn parameters(&self) -> &EmbeddedParameters {
+        &self.parameters
+    }
+
+    fn boundary(&self, x: f64, y: f64, t: f64) -> (f64, f64) {
+        self.boundary_velocity
+            .as_ref()
+            .map_or((0.0, 0.0), |f| f(x, y, t))
+    }
+
+    /// Build the mask (if a body is set) and stamp the `t = time` boundary
+    /// data and ghost values onto the field. Called lazily by the first
+    /// `advance`; call it explicitly to inspect the mask or to run
+    /// diagnostics on the initial field.
+    pub fn initialize(&mut self, field: &mut FlowField) -> CfdResult<()> {
+        let (nx, ny, dx, dy) = field.grid_info();
+        if let Some(body) = &self.body {
+            if self.mask.is_none() {
+                self.mask = Some(EmbeddedMask::build(body, nx, ny, dx, dy, self.time)?);
+            }
+        }
+        let t = self.time;
+        self.apply_boundary_normals(field, t);
+        if let (Some(body), Some(mask)) = (&self.body, &self.mask) {
+            mask.impose(body, &mut field.u, &mut field.v, t);
+        }
+        self.initialized = true;
+        Ok(())
+    }
+
+    /// Write the prescribed normal velocities onto the boundary faces at
+    /// time `t`. Outlet faces are unknowns and are left alone.
+    fn apply_boundary_normals(&self, field: &mut FlowField, t: f64) {
+        let (nx, ny, dx, dy) = field.grid_info();
+        let b = self.parameters.boundaries;
+        let outlet = SideBoundary::PressureOutlet;
+        for j in 0..ny {
+            let y = (j as f64 + 0.5) * dy;
+            if b.left != outlet {
+                field.u[(j, 0)] = self.boundary(0.0, y, t).0;
+            }
+            if b.right != outlet {
+                field.u[(j, nx)] = self.boundary(nx as f64 * dx, y, t).0;
+            }
+        }
+        for i in 0..nx {
+            let x = (i as f64 + 0.5) * dx;
+            if b.bottom != outlet {
+                field.v[(0, i)] = self.boundary(x, 0.0, t).1;
+            }
+            if b.top != outlet {
+                field.v[(ny, i)] = self.boundary(x, ny as f64 * dy, t).1;
+            }
+        }
+    }
+
+    #[inline]
+    fn u_is_fluid(&self, j: usize, i: usize) -> bool {
+        self.mask
+            .as_ref()
+            .is_none_or(|m| m.u_kind(j, i) == FaceKind::Fluid)
+    }
+
+    #[inline]
+    fn v_is_fluid(&self, j: usize, i: usize) -> bool {
+        self.mask
+            .as_ref()
+            .is_none_or(|m| m.v_kind(j, i) == FaceKind::Fluid)
+    }
+
+    #[inline]
+    fn cell_is_fluid(&self, j: usize, i: usize) -> bool {
+        self.mask.as_ref().is_none_or(|m| m.is_fluid_cell(j, i))
+    }
+
+    fn upwind(face_velocity: f64, upstream: f64, downstream: f64) -> f64 {
+        if face_velocity >= 0.0 {
+            upstream
+        } else {
+            downstream
+        }
+    }
+
+    /// Explicit momentum predictor on the fluid faces, expression for
+    /// expression the fixed-grid PISO's (so the no-body case is identical
+    /// to the bit), plus the slip-wall / outlet arms of the ALE solver on
+    /// the domain sides. Non-fluid faces keep their prescribed values.
+    #[allow(clippy::too_many_lines)]
+    fn momentum_predictor(&self, field: &mut FlowField, dt: f64, t_old: f64) -> CfdResult<()> {
+        let (nx, ny, dx, dy) = field.grid_info();
+        let rho = self.config.density;
+        let nu = self.config.viscosity / rho;
+        let b = self.parameters.boundaries;
+        let velocity = SideBoundary::Velocity;
+
+        for j in 0..ny {
+            for i in 1..nx {
+                if !self.u_is_fluid(j, i) {
+                    continue;
+                }
+                let uo = &field.u_old;
+                let vo = &field.v_old;
+                let u_p = uo[(j, i)];
+
+                let ue_face = 0.5 * (uo[(j, i)] + uo[(j, i + 1)]);
+                let uw_face = 0.5 * (uo[(j, i - 1)] + uo[(j, i)]);
+
+                let south_is_wall = j == 0;
+                let north_is_wall = j + 1 == ny;
+
+                // Transverse face velocities from the stored v faces — on a
+                // domain side these are the prescribed boundary normals
+                // (zero on a wall, the outflow on an outlet). The fixed-grid
+                // PISO zeroes them on its walls, which is the same number on
+                // a wall and wrong on an outlet: the outgoing mass flux must
+                // carry momentum out, or the last row accumulates it.
+                let vn_face = 0.5 * (vo[(j + 1, i - 1)] + vo[(j + 1, i)]);
+                let vs_face = 0.5 * (vo[(j, i - 1)] + vo[(j, i)]);
+
+                // Upwind value across a domain side: the boundary function's
+                // tangential value on a Velocity side, the interior value
+                // otherwise (zero-gradient).
+                let beyond_north = if b.top == velocity {
+                    self.boundary(i as f64 * dx, ny as f64 * dy, t_old).0
+                } else {
+                    u_p
+                };
+                let beyond_south = if b.bottom == velocity {
+                    self.boundary(i as f64 * dx, 0.0, t_old).0
+                } else {
+                    u_p
+                };
+
+                let conv_x = (ue_face * Self::upwind(ue_face, uo[(j, i)], uo[(j, i + 1)])
+                    - uw_face * Self::upwind(uw_face, uo[(j, i - 1)], uo[(j, i)]))
+                    / dx;
+                let conv_y = (vn_face
+                    * if north_is_wall {
+                        Self::upwind(vn_face, u_p, beyond_north)
+                    } else {
+                        Self::upwind(vn_face, uo[(j, i)], uo[(j + 1, i)])
+                    }
+                    - vs_face
+                        * if south_is_wall {
+                            Self::upwind(vs_face, beyond_south, u_p)
+                        } else {
+                            Self::upwind(vs_face, uo[(j - 1, i)], uo[(j, i)])
+                        })
+                    / dy;
+
+                let diff_x = nu * (uo[(j, i + 1)] - 2.0 * u_p + uo[(j, i - 1)]) / (dx * dx);
+
+                // Wall-adjacent diffusive fluxes act over half a cell on a
+                // Velocity side; a slip wall or outlet carries no shear.
+                let flux_north = if north_is_wall {
+                    if b.top == velocity {
+                        let u_wall = self.boundary(i as f64 * dx, ny as f64 * dy, t_old).0;
+                        nu * (u_wall - u_p) / (0.5 * dy)
+                    } else {
+                        0.0
+                    }
+                } else {
+                    nu * (uo[(j + 1, i)] - u_p) / dy
+                };
+                let flux_south = if south_is_wall {
+                    if b.bottom == velocity {
+                        let u_wall = self.boundary(i as f64 * dx, 0.0, t_old).0;
+                        nu * (u_p - u_wall) / (0.5 * dy)
+                    } else {
+                        0.0
+                    }
+                } else {
+                    nu * (u_p - uo[(j - 1, i)]) / dy
+                };
+                let diff_y = (flux_north - flux_south) / dy;
+
+                let pressure_gradient = -(field.p[(j, i)] - field.p[(j, i - 1)]) / (rho * dx);
+
+                let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
+                    f(i as f64 * dx, (j as f64 + 0.5) * dy, t_old).0 / rho
+                });
+
+                field.u[(j, i)] = u_p
+                    + dt * (-conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force);
+            }
+        }
+
+        for j in 1..ny {
+            for i in 0..nx {
+                if !self.v_is_fluid(j, i) {
+                    continue;
+                }
+                let uo = &field.u_old;
+                let vo = &field.v_old;
+                let v_p = vo[(j, i)];
+
+                let vn_face = 0.5 * (vo[(j, i)] + vo[(j + 1, i)]);
+                let vs_face = 0.5 * (vo[(j - 1, i)] + vo[(j, i)]);
+
+                let west_is_wall = i == 0;
+                let east_is_wall = i + 1 == nx;
+
+                let ue_face = 0.5 * (uo[(j - 1, i + 1)] + uo[(j, i + 1)]);
+                let uw_face = 0.5 * (uo[(j - 1, i)] + uo[(j, i)]);
+                let beyond_east = if b.right == velocity {
+                    self.boundary(nx as f64 * dx, j as f64 * dy, t_old).1
+                } else {
+                    v_p
+                };
+                let beyond_west = if b.left == velocity {
+                    self.boundary(0.0, j as f64 * dy, t_old).1
+                } else {
+                    v_p
+                };
+
+                let conv_y = (vn_face * Self::upwind(vn_face, vo[(j, i)], vo[(j + 1, i)])
+                    - vs_face * Self::upwind(vs_face, vo[(j - 1, i)], vo[(j, i)]))
+                    / dy;
+                let conv_x = (ue_face
+                    * if east_is_wall {
+                        Self::upwind(ue_face, v_p, beyond_east)
+                    } else {
+                        Self::upwind(ue_face, vo[(j, i)], vo[(j, i + 1)])
+                    }
+                    - uw_face
+                        * if west_is_wall {
+                            Self::upwind(uw_face, beyond_west, v_p)
+                        } else {
+                            Self::upwind(uw_face, vo[(j, i - 1)], vo[(j, i)])
+                        })
+                    / dx;
+
+                let diff_y = nu * (vo[(j + 1, i)] - 2.0 * v_p + vo[(j - 1, i)]) / (dy * dy);
+
+                let flux_east = if east_is_wall {
+                    if b.right == velocity {
+                        let v_wall = self.boundary(nx as f64 * dx, j as f64 * dy, t_old).1;
+                        nu * (v_wall - v_p) / (0.5 * dx)
+                    } else {
+                        0.0
+                    }
+                } else {
+                    nu * (vo[(j, i + 1)] - v_p) / dx
+                };
+                let flux_west = if west_is_wall {
+                    if b.left == velocity {
+                        let v_wall = self.boundary(0.0, j as f64 * dy, t_old).1;
+                        nu * (v_p - v_wall) / (0.5 * dx)
+                    } else {
+                        0.0
+                    }
+                } else {
+                    nu * (v_p - vo[(j, i - 1)]) / dx
+                };
+                let diff_x = (flux_east - flux_west) / dx;
+
+                let pressure_gradient = -(field.p[(j, i)] - field.p[(j - 1, i)]) / (rho * dy);
+
+                let body_force = self.momentum_source.as_ref().map_or(0.0, |f| {
+                    f((i as f64 + 0.5) * dx, j as f64 * dy, t_old).1 / rho
+                });
+
+                field.v[(j, i)] = v_p
+                    + dt * (-conv_x - conv_y + diff_x + diff_y + pressure_gradient + body_force);
+            }
+        }
+
+        // Outlet faces: zero-gradient predictor value, corrected by the
+        // projection.
+        if b.left == SideBoundary::PressureOutlet {
+            for j in 0..ny {
+                field.u[(j, 0)] = field.u[(j, 1)];
+            }
+        }
+        if b.right == SideBoundary::PressureOutlet {
+            for j in 0..ny {
+                field.u[(j, nx)] = field.u[(j, nx - 1)];
+            }
+        }
+        if b.bottom == SideBoundary::PressureOutlet {
+            for i in 0..nx {
+                field.v[(0, i)] = field.v[(1, i)];
+            }
+        }
+        if b.top == SideBoundary::PressureOutlet {
+            for i in 0..nx {
+                field.v[(ny, i)] = field.v[(ny - 1, i)];
+            }
+        }
+
+        field.copy_to_starred();
+        Ok(())
+    }
+
+    /// One projection on the fluid cells: the fixed-grid PISO's, with a
+    /// zero coefficient across every prescribed face (domain Velocity /
+    /// SlipWall sides and every non-fluid interior face), a Dirichlet `p' =
+    /// 0` half a cell beyond an outlet face, and the anchor on the first
+    /// fluid cell when no outlet exists. Returns the normalised mass
+    /// imbalance of the corrected field over the fluid cells.
+    #[allow(clippy::too_many_lines)]
+    fn project(&self, field: &mut FlowField, dt: f64) -> CfdResult<f64> {
+        let (nx, ny, dx, dy) = field.grid_info();
+        let rho = self.config.density;
+        let b = self.parameters.boundaries;
+        let outlet = SideBoundary::PressureOutlet;
+        let any_outlet = [b.left, b.right, b.bottom, b.top].contains(&outlet);
+        let anchor = self.mask.as_ref().map_or((1, 1), EmbeddedMask::anchor);
+
+        field.p_prime.fill(0.0);
+
+        let mut source_scale = 0.0;
+        for j in 0..ny {
+            for i in 0..nx {
+                if !self.cell_is_fluid(j, i) {
+                    field.sp[(j, i)] = 0.0;
+                    continue;
+                }
+                let divergence_flux = rho
+                    * ((field.u_star[(j, i + 1)] - field.u_star[(j, i)]) * dy
+                        + (field.v_star[(j + 1, i)] - field.v_star[(j, i)]) * dx);
+                field.sp[(j, i)] = -divergence_flux;
+                source_scale += divergence_flux.abs();
+            }
+        }
+
+        let ae_interior = dt * dy / dx;
+        let an_interior = dt * dx / dy;
+        // Outlet face: p' = 0 half a cell away.
+        let ae_outlet = dt * dy / (0.5 * dx);
+        let an_outlet = dt * dx / (0.5 * dy);
+
+        let reference_flux = rho * self.config.reference_velocity * self.config.reference_length;
+        let inner_stop =
+            (1e-2 * source_scale).max(0.1 * self.parameters.tolerance * reference_flux) + 1e-14;
+        let omega = 2.0 / (1.0 + (std::f64::consts::PI / nx.max(ny) as f64).sin());
+        for _sweep in 0..2000 {
+            let mut residual = 0.0;
+            for j in 0..ny {
+                for i in 0..nx {
+                    if !self.cell_is_fluid(j, i) {
+                        continue;
+                    }
+                    if !any_outlet && (j, i) == anchor {
+                        field.p_prime[(j, i)] = 0.0;
+                        continue;
+                    }
+
+                    // A coefficient is zero exactly when the face is
+                    // prescribed: a domain side with velocity data, or a
+                    // non-fluid interior face.
+                    let ae = if i + 1 == nx {
+                        if b.right == outlet { ae_outlet } else { 0.0 }
+                    } else if self.u_is_fluid(j, i + 1) {
+                        ae_interior
+                    } else {
+                        0.0
+                    };
+                    let aw = if i == 0 {
+                        if b.left == outlet { ae_outlet } else { 0.0 }
+                    } else if self.u_is_fluid(j, i) {
+                        ae_interior
+                    } else {
+                        0.0
+                    };
+                    let an = if j + 1 == ny {
+                        if b.top == outlet { an_outlet } else { 0.0 }
+                    } else if self.v_is_fluid(j + 1, i) {
+                        an_interior
+                    } else {
+                        0.0
+                    };
+                    let as_ = if j == 0 {
+                        if b.bottom == outlet { an_outlet } else { 0.0 }
+                    } else if self.v_is_fluid(j, i) {
+                        an_interior
+                    } else {
+                        0.0
+                    };
+                    let ap = ae + aw + an + as_;
+                    if ap == 0.0 {
+                        // An isolated fluid cell enclosed by prescribed
+                        // faces has no equation; leave p' = 0 there.
+                        continue;
+                    }
+
+                    let east = if i + 1 < nx {
+                        ae * field.p_prime[(j, i + 1)]
+                    } else {
+                        0.0
+                    };
+                    let west = if i > 0 {
+                        aw * field.p_prime[(j, i - 1)]
+                    } else {
+                        0.0
+                    };
+                    let north = if j + 1 < ny {
+                        an * field.p_prime[(j + 1, i)]
+                    } else {
+                        0.0
+                    };
+                    let south = if j > 0 {
+                        as_ * field.p_prime[(j - 1, i)]
+                    } else {
+                        0.0
+                    };
+
+                    let rhs = field.sp[(j, i)] + east + west + north + south;
+                    let p_old = field.p_prime[(j, i)];
+                    residual += (rhs - ap * p_old).abs();
+                    field.p_prime[(j, i)] = (1.0 - omega) * p_old + omega * rhs / ap;
+                }
+            }
+            if residual < inner_stop {
+                break;
+            }
+        }
+
+        // Correct exactly the faces the equations treated as correctable:
+        // fluid interior faces, and outlet faces against p' = 0 outside.
+        for j in 0..ny {
+            for i in 1..nx {
+                if self.u_is_fluid(j, i) {
+                    let dp_dx = (field.p_prime[(j, i)] - field.p_prime[(j, i - 1)]) / dx;
+                    field.u[(j, i)] = field.u_star[(j, i)] - (dt / rho) * dp_dx;
+                }
+            }
+            if b.left == outlet {
+                let dp_dx = (field.p_prime[(j, 0)] - 0.0) / (0.5 * dx);
+                field.u[(j, 0)] = field.u_star[(j, 0)] - (dt / rho) * dp_dx;
+            }
+            if b.right == outlet {
+                let dp_dx = (0.0 - field.p_prime[(j, nx - 1)]) / (0.5 * dx);
+                field.u[(j, nx)] = field.u_star[(j, nx)] - (dt / rho) * dp_dx;
+            }
+        }
+        for i in 0..nx {
+            for j in 1..ny {
+                if self.v_is_fluid(j, i) {
+                    let dp_dy = (field.p_prime[(j, i)] - field.p_prime[(j - 1, i)]) / dy;
+                    field.v[(j, i)] = field.v_star[(j, i)] - (dt / rho) * dp_dy;
+                }
+            }
+            if b.bottom == outlet {
+                let dp_dy = (field.p_prime[(0, i)] - 0.0) / (0.5 * dy);
+                field.v[(0, i)] = field.v_star[(0, i)] - (dt / rho) * dp_dy;
+            }
+            if b.top == outlet {
+                let dp_dy = (0.0 - field.p_prime[(ny - 1, i)]) / (0.5 * dy);
+                field.v[(ny, i)] = field.v_star[(ny, i)] - (dt / rho) * dp_dy;
+            }
+        }
+
+        for j in 0..ny {
+            for i in 0..nx {
+                if self.cell_is_fluid(j, i) {
+                    field.p[(j, i)] += field.p_prime[(j, i)];
+                }
+            }
+        }
+
+        let mut mass_imbalance = 0.0;
+        for j in 0..ny {
+            for i in 0..nx {
+                if !self.cell_is_fluid(j, i) {
+                    continue;
+                }
+                let divergence_flux = rho
+                    * ((field.u[(j, i + 1)] - field.u[(j, i)]) * dy
+                        + (field.v[(j + 1, i)] - field.v[(j, i)]) * dx);
+                mass_imbalance += divergence_flux.abs();
+            }
+        }
+
+        Ok(if reference_flux > 0.0 {
+            mass_imbalance / reference_flux
+        } else {
+            mass_imbalance
+        })
+    }
+
+    /// Advance one step of `dt`.
+    pub async fn advance(&mut self, field: &mut FlowField, dt: f64) -> CfdResult<EmbeddedResult> {
+        if dt <= 0.0 || !dt.is_finite() {
+            return Err(CfdError::invalid_parameter(format!(
+                "time step must be positive and finite, got {dt}"
+            )));
+        }
+        if !self.initialized {
+            self.initialize(field)?;
+        }
+        let start_time = std::time::Instant::now();
+        let t_old = self.time;
+        let t_new = t_old + dt;
+
+        // The start-of-step state is whatever the previous step left on the
+        // boundary and ghost faces — no re-stamping (see module docs).
+        field.update_old_values();
+        self.momentum_predictor(field, dt, t_old)?;
+        // Boundary data for the new interval; the predictor's `u` holds the
+        // old boundary values until now.
+        self.apply_boundary_normals(field, t_new);
+        field.copy_to_starred();
+
+        let mut residual_history = Vec::new();
+        let mut total_corrector_steps = 0;
+        let mut final_residual = f64::INFINITY;
+        for _corrector in 0..self.parameters.corrector_steps.max(1) {
+            let mass_residual = self.project(field, dt)?;
+            residual_history.push(mass_residual);
+            final_residual = mass_residual;
+            total_corrector_steps += 1;
+            if mass_residual < self.parameters.tolerance {
+                break;
+            }
+            field.copy_to_starred();
+        }
+
+        // Ghost faces follow the corrected fluid field; they are the
+        // stencil and flux data of the next step.
+        let ghost_correction = match (&self.body, &self.mask) {
+            (Some(body), Some(mask)) => mask.impose(body, &mut field.u, &mut field.v, t_new),
+            _ => 0.0,
+        };
+
+        self.time = t_new;
+        Ok(EmbeddedResult {
+            solver_result: SolverResult {
+                converged: final_residual < self.parameters.tolerance,
+                iterations: total_corrector_steps,
+                final_residual,
+                residual_history,
+                solve_time: start_time.elapsed(),
+            },
+            corrector_steps_performed: total_corrector_steps,
+            ghost_correction,
+        })
+    }
+}

crates/specialized/rtx-cfd/src/solvers/incompressible/mod.rs

@@ -13,6 +13,10 @@ use crate::{CfdConfig, CfdError, CfdResult};
 pub mod ale;
 /// Boundary conditions
 pub mod boundary_conditions;
+/// PISO on the fixed grid with an embedded body
+pub mod embedded;
+/// Embedded-body geometry, classification and loads
+pub mod embedded_body;
 /// Flow field data structures
 pub mod flow_field;
 /// PISO algorithm implementation
@@ -33,6 +37,8 @@ pub use ale::{
 pub use boundary_conditions::{
     BoundaryCondition, BoundaryConditions, BoundaryLocation, BoundaryType,
 };
+pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult};
+pub use embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample};
 pub use flow_field::FlowField;
 pub use piso::{PisoParameters, PisoResult, PisoSolver};
 #[cfg(feature = "cuda")]

crates/specialized/rtx-cfd/tests/embedded_mms.rs

@@ -0,0 +1,471 @@
+//! Code verification of the embedded-boundary PISO solver.
+//!
+//! Two claims, in the order they must be established:
+//!
+//! 1. **With no body it IS the fixed-grid PISO** — the same manufactured
+//!    problem marched by `PisoSolver` and by `EmbeddedPisoSolver` must
+//!    produce bit-identical fields, because the predictor and projection
+//!    are the same expressions and the only additions are masked out.
+//!
+//! 2. **With an embedded circle the manufactured solution is recovered at
+//!    the discretisation's order**, the field is divergence-free on every
+//!    fluid cell, the compatibility correction shrinks with the mesh, and
+//!    the force on the circle by *both* load routes — surface-stress
+//!    reconstruction and a control-volume momentum balance — converges to
+//!    the exact surface integral of the manufactured stress.
+//!
+//! The manufactured field, source and grid convention are those of
+//! `tests/mms_piso.rs` (`u = sin(pi x) cos(pi y)`, `v = -cos(pi x) sin(pi y)`,
+//! `p = sin(pi x) sin(pi y)`); the circle (centre (0.6, 0.45), r = 0.2 — off-centre, so the exact force is not zero by symmetry)
+//! carries the exact field as its surface velocity, so the embedded wall is
+//! a Dirichlet boundary on a curve that cuts the grid arbitrarily — which
+//! is exactly what the ghost reconstruction has to get right.
+
+use rtx_cfd::solvers::incompressible::{
+    BoundaryConditions, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver, FaceKind, FlowField,
+    IncompressibleSolver, PisoParameters, PisoSolver,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+use std::f64::consts::PI;
+
+const RHO: f64 = 1.0;
+const MU: f64 = 0.05;
+const CX: f64 = 0.6;
+const CY: f64 = 0.45;
+const R: f64 = 0.2;
+
+fn u_exact(x: f64, y: f64) -> f64 {
+    (PI * x).sin() * (PI * y).cos()
+}
+
+fn v_exact(x: f64, y: f64) -> f64 {
+    -(PI * x).cos() * (PI * y).sin()
+}
+
+fn p_exact(x: f64, y: f64) -> f64 {
+    (PI * x).sin() * (PI * y).sin()
+}
+
+fn source(x: f64, y: f64) -> (f64, f64) {
+    let fx = RHO * 0.5 * PI * (2.0 * PI * x).sin()
+        + 2.0 * PI * PI * MU * u_exact(x, y)
+        + PI * (PI * x).cos() * (PI * y).sin();
+    let fy = RHO * 0.5 * PI * (2.0 * PI * y).sin()
+        + 2.0 * PI * PI * MU * v_exact(x, y)
+        + PI * (PI * x).sin() * (PI * y).cos();
+    (fx, fy)
+}
+
+/// Exact force on the circle from the manufactured stress,
+/// `F = oint (-p I + mu (grad u + grad u^T)) n ds`, by fine quadrature, and
+/// the momentum flux through the circle `M = oint rho u (u . n) ds` — the
+/// manufactured surface velocity has a normal component, so the "body" is
+/// porous. A control-volume balance around it therefore measures `F - M`
+/// (the surface-stress route measures `F`); for a rigid no-slip body `M`
+/// is zero and the two routes measure the same thing.
+fn exact_force_and_flux() -> ((f64, f64), (f64, f64)) {
+    let n = 20_000;
+    let (mut fx, mut fy) = (0.0, 0.0);
+    let (mut mx, mut my) = (0.0, 0.0);
+    for k in 0..n {
+        let theta = (k as f64 + 0.5) * 2.0 * PI / n as f64;
+        let (s, c) = theta.sin_cos();
+        let (x, y) = (CX + R * c, CY + R * s);
+        let ux = PI * (PI * x).cos() * (PI * y).cos();
+        let uy = -PI * (PI * x).sin() * (PI * y).sin();
+        let vx = PI * (PI * x).sin() * (PI * y).sin();
+        let vy = -PI * (PI * x).cos() * (PI * y).cos();
+        let p = p_exact(x, y);
+        let sxx = -p + 2.0 * MU * ux;
+        let syy = -p + 2.0 * MU * vy;
+        let sxy = MU * (uy + vx);
+        let ds = 2.0 * PI * R / n as f64;
+        fx += (sxx * c + sxy * s) * ds;
+        fy += (sxy * c + syy * s) * ds;
+        let (u, v) = (u_exact(x, y), v_exact(x, y));
+        let un = u * c + v * s;
+        mx += RHO * u * un * ds;
+        my += RHO * v * un * ds;
+    }
+    ((fx, fy), (mx, my))
+}
+
+/// The manufactured field on the box boundary with the normal components
+/// snapped to their exact analytic zero: `sin(pi)` evaluates to 1.2e-16,
+/// and a 1e-16 through-flow is enough to separate the two solvers at the
+/// last bit (the embedded solver carries boundary fluxes faithfully).
+fn boundary_exact(x: f64, y: f64) -> (f64, f64) {
+    let u = if x <= 0.0 || x >= 1.0 {
+        0.0
+    } else {
+        u_exact(x, y)
+    };
+    let v = if y <= 0.0 || y >= 1.0 {
+        0.0
+    } else {
+        v_exact(x, y)
+    };
+    (u, v)
+}
+
+fn config() -> CfdConfig {
+    CfdConfig::new()
+        .with_density(RHO)
+        .with_viscosity(MU)
+        .with_reference_velocity(1.0)
+        .with_reference_length(1.0)
+}
+
+fn initial_field(n: usize) -> CfdResult<FlowField> {
+    let dx = 1.0 / n as f64;
+    let mut field = FlowField::new(n, n, dx, dx)?;
+    for j in 0..n {
+        let y = (j as f64 + 0.5) * dx;
+        field.u[(j, 0)] = boundary_exact(0.0, y).0;
+        field.u[(j, n)] = boundary_exact(1.0, y).0;
+    }
+    for i in 0..n {
+        let x = (i as f64 + 0.5) * dx;
+        field.v[(0, i)] = boundary_exact(x, 0.0).1;
+        field.v[(n, i)] = boundary_exact(x, 1.0).1;
+    }
+    Ok(field)
+}
+
+fn time_step(n: usize) -> f64 {
+    let dx = 1.0 / n as f64;
+    let nu = MU / RHO;
+    0.4 * (dx * dx / (4.0 * nu)).min(dx)
+}
+
+/// Claim 1: no body, velocity on every side — the two solvers must agree
+/// to the bit over a couple of hundred steps from the same start.
+#[tokio::test]
+async fn without_a_body_the_embedded_solver_is_piso_to_the_bit() -> CfdResult<()> {
+    let n = 16;
+    let dt = time_step(n);
+
+    let mut piso = PisoSolver::new(
+        config(),
+        PisoParameters {
+            corrector_steps: 2,
+            time_step: dt,
+            tolerance: 1e-8,
+        },
+    )?;
+    piso.set_momentum_source(source);
+    piso.set_wall_velocity(boundary_exact);
+
+    let mut embedded = EmbeddedPisoSolver::new(
+        config(),
+        EmbeddedParameters {
+            corrector_steps: 2,
+            tolerance: 1e-8,
+            ..EmbeddedParameters::default()
+        },
+    )?;
+    embedded.set_momentum_source(|x, y, _| source(x, y));
+    embedded.set_boundary_velocity(|x, y, _| boundary_exact(x, y));
+
+    let mut a = initial_field(n)?;
+    let mut b = initial_field(n)?;
+    let empty = BoundaryConditions::new();
+    for _ in 0..200 {
+        piso.solve_time_step(&mut a, &empty, dt).await?;
+        embedded.advance(&mut b, dt).await?;
+    }
+    let mut max_diff: f64 = 0.0;
+    for (x, y) in a.u.iter().zip(b.u.iter()) {
+        max_diff = max_diff.max((x - y).abs());
+    }
+    for (x, y) in a.v.iter().zip(b.v.iter()) {
+        max_diff = max_diff.max((x - y).abs());
+    }
+    for (x, y) in a.p.iter().zip(b.p.iter()) {
+        max_diff = max_diff.max((x - y).abs());
+    }
+    assert!(
+        max_diff == 0.0,
+        "embedded solver without a body differs from PISO by {max_diff:.3e}"
+    );
+    Ok(())
+}
+
+struct Measurement {
+    l2_velocity: f64,
+    l2_pressure: f64,
+    max_div: f64,
+    ghost_correction: f64,
+    force_surface: (f64, f64),
+    skipped_samples: usize,
+    force_cv: (f64, f64),
+}
+
+/// March the manufactured problem with the embedded circle to steady
+/// state on an `n` by `n` grid and measure everything.
+async fn measure(n: usize) -> CfdResult<Measurement> {
+    let dx = 1.0 / n as f64;
+    let dt = time_step(n);
+
+    let mut solver = EmbeddedPisoSolver::new(
+        config(),
+        EmbeddedParameters {
+            corrector_steps: 2,
+            tolerance: 1e-8,
+            ..EmbeddedParameters::default()
+        },
+    )?;
+    solver.set_momentum_source(|x, y, _| source(x, y));
+    solver.set_boundary_velocity(|x, y, _| boundary_exact(x, y));
+    solver.set_body(
+        EmbeddedBody::circle(CX, CY, R)
+            .with_surface_velocity(|x, y, _| (u_exact(x, y), v_exact(x, y))),
+    );
+
+    let mut field = initial_field(n)?;
+    solver.initialize(&mut field)?;
+
+    let mut steady_residual = f64::INFINITY;
+    let mut last_correction = 0.0;
+    for _step in 0..200_000 {
+        let u_before = field.u.clone();
+        let v_before = field.v.clone();
+        let result = solver.advance(&mut field, dt).await?;
+        last_correction = result.ghost_correction;
+
+        let mut max_change: f64 = 0.0;
+        for (a, b) in field.u.iter().zip(u_before.iter()) {
+            max_change = max_change.max((a - b).abs());
+        }
+        for (a, b) in field.v.iter().zip(v_before.iter()) {
+            max_change = max_change.max((a - b).abs());
+        }
+        steady_residual = max_change / dt;
+        if steady_residual < 1e-6 {
+            break;
+        }
+    }
+    assert!(
+        steady_residual < 1e-6,
+        "embedded PISO did not reach a steady state at n = {n}: |du/dt| = {steady_residual:.3e}"
+    );
+
+    let mask = solver.mask().expect("mask built");
+
+    // Velocity error over the fluid faces, pressure error over the fluid
+    // cells (mean-shifted: the level is arbitrary), divergence on every
+    // fluid cell.
+    let mut squared = 0.0;
+    let mut volume = 0.0;
+    for j in 0..n {
+        for i in 1..n {
+            if mask.u_kind(j, i) == FaceKind::Fluid {
+                let e = field.u[(j, i)] - u_exact(i as f64 * dx, (j as f64 + 0.5) * dx);
+                squared += e * e * dx * dx;
+                volume += dx * dx;
+            }
+        }
+    }
+    for j in 1..n {
+        for i in 0..n {
+            if mask.v_kind(j, i) == FaceKind::Fluid {
+                let e = field.v[(j, i)] - v_exact((i as f64 + 0.5) * dx, j as f64 * dx);
+                squared += e * e * dx * dx;
+                volume += dx * dx;
+            }
+        }
+    }
+    let l2_velocity = (squared / volume).sqrt();
+
+    let mut diff_sum = 0.0;
+    let mut cells = 0usize;
+    for j in 0..n {
+        for i in 0..n {
+            if mask.is_fluid_cell(j, i) {
+                diff_sum += field.p[(j, i)] - p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dx);
+                cells += 1;
+            }
+        }
+    }
+    let shift = diff_sum / cells as f64;
+    let mut p_sq = 0.0;
+    let mut max_div: f64 = 0.0;
+    for j in 0..n {
+        for i in 0..n {
+            if mask.is_fluid_cell(j, i) {
+                let e =
+                    field.p[(j, i)] - shift - p_exact((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dx);
+                p_sq += e * e;
+                let div = (field.u[(j, i + 1)] - field.u[(j, i)]) / dx
+                    + (field.v[(j + 1, i)] - field.v[(j, i)]) / dx;
+                max_div = max_div.max(div.abs());
+            }
+        }
+    }
+    let l2_pressure = (p_sq / cells as f64).sqrt();
+
+    let body = solver.body().expect("body set");
+    let surface = mask.surface_force(
+        body,
+        &field.u,
+        &field.v,
+        &field.p,
+        MU,
+        solver.time(),
+        0.5 * dx,
+    );
+    // Control volume: the middle three quarters of the box, whole cells.
+    let i0 = n / 8;
+    let i1 = n - n / 8;
+    let src = |x: f64, y: f64| source(x, y);
+    let force_cv = mask.control_volume_force(
+        &field.u,
+        &field.v,
+        &field.p,
+        &field.u_old,
+        &field.v_old,
+        dt,
+        RHO,
+        MU,
+        Some(&src),
+        (i0, i1, i0, i1),
+    );
+
+    Ok(Measurement {
+        l2_velocity,
+        l2_pressure,
+        max_div,
+        ghost_correction: last_correction.abs(),
+        force_surface: (surface.fx, surface.fy),
+        skipped_samples: surface.skipped,
+        force_cv,
+    })
+}
+
+/// Claim 2. Measured (16 -> 32 -> 64): L2 velocity 1.607e-2, 8.489e-3,
+/// 4.341e-3 — orders 0.92, 0.97 (PISO without a body: 0.85, 0.91);
+/// L2 pressure orders 0.96, 0.90; max |div u| <= 9e-8 on every fluid cell;
+/// compatibility correction 6.2e-4, 1.1e-5, 2.8e-5; force error relative
+/// to the exact |F|: surface route 0.52, 0.29, 0.15, control-volume route
+/// 0.61, 0.30, 0.15 — both first order, both from the same solution, by
+/// two unrelated readings of it. The structure of what must hold was fixed
+/// before the numbers were known:
+/// - the velocity error falls at the scheme's order (first-order upwind:
+///   approaching 1; the ghost treatment must not drag it below),
+/// - every fluid cell is divergence-free,
+/// - the compatibility correction shrinks with the mesh,
+/// - both force routes converge to the exact force.
+#[tokio::test]
+async fn embedded_circle_recovers_the_manufactured_solution() -> CfdResult<()> {
+    let resolutions = [16usize, 32, 64];
+    let mut measurements = Vec::new();
+    for &n in &resolutions {
+        measurements.push(measure(n).await?);
+    }
+    let ((fx_exact, fy_exact), (mx, my)) = exact_force_and_flux();
+    let f_scale = (fx_exact * fx_exact + fy_exact * fy_exact).sqrt();
+    let (fx_cv_exact, fy_cv_exact) = (fx_exact - mx, fy_exact - my);
+    println!(
+        "  exact force on the circle: ({fx_exact:.6e}, {fy_exact:.6e}); momentum flux through it \
+         ({mx:.6e}, {my:.6e}); the control-volume route measures ({fx_cv_exact:.6e}, {fy_cv_exact:.6e})"
+    );
+
+    let errors: Vec<f64> = measurements.iter().map(|m| m.l2_velocity).collect();
+    let p_errors: Vec<f64> = measurements.iter().map(|m| m.l2_pressure).collect();
+    let rates: Vec<f64> = errors.windows(2).map(|w| (w[0] / w[1]).log2()).collect();
+    let p_rates: Vec<f64> = p_errors.windows(2).map(|w| (w[0] / w[1]).log2()).collect();
+
+    let mut surface_errors = Vec::new();
+    let mut cv_errors = Vec::new();
+    for (k, (m, &n)) in measurements.iter().zip(&resolutions).enumerate() {
+        let rate = if k == 0 {
+            String::from("    -")
+        } else {
+            format!("{:5.2}", rates[k - 1])
+        };
+        let p_rate = if k == 0 {
+            String::from("    -")
+        } else {
+            format!("{:5.2}", p_rates[k - 1])
+        };
+        let surface = (m.force_surface.0 - fx_exact).hypot(m.force_surface.1 - fy_exact) / f_scale;
+        let cv = (m.force_cv.0 - fx_cv_exact).hypot(m.force_cv.1 - fy_cv_exact) / f_scale;
+        println!(
+            "  n = {n:3}   L2 u {:.4e} (order {rate})   L2 p {:.4e} (order {p_rate})   \
+             max div {:.2e}   ghost corr {:.2e}   F_surface ({:.5e}, {:.5e}) rel {:.3e} skipped {}   \
+             F_cv ({:.5e}, {:.5e}) rel {:.3e}",
+            m.l2_velocity,
+            m.l2_pressure,
+            m.max_div,
+            m.ghost_correction,
+            m.force_surface.0,
+            m.force_surface.1,
+            surface,
+            m.skipped_samples,
+            m.force_cv.0,
+            m.force_cv.1,
+            cv
+        );
+        surface_errors.push(surface);
+        cv_errors.push(cv);
+    }
+
+    assert!(
+        errors.windows(2).all(|w| w[1] < w[0]),
+        "velocity error must fall under refinement: {errors:?}"
+    );
+    for (k, &rate) in rates.iter().enumerate() {
+        assert!(
+            rate > 0.75,
+            "refinement {} -> {}: velocity order {rate:.3} below what first-order upwind delivers \
+             without a body (0.85, 0.91) — the embedded treatment is polluting the order. \
+             Errors {errors:?}",
+            resolutions[k],
+            resolutions[k + 1]
+        );
+        assert!(
+            rate < 2.3,
+            "velocity order {rate:.3} above the scheme's — suspect the measure"
+        );
+    }
+    assert!(
+        p_errors.windows(2).all(|w| w[1] < w[0]),
+        "pressure error must fall under refinement: {p_errors:?}"
+    );
+    for m in &measurements {
+        assert!(
+            m.max_div < 1e-5,
+            "a fluid cell is not divergence-free: max |div u| = {:.3e}",
+            m.max_div
+        );
+    }
+    let corrections: Vec<f64> = measurements.iter().map(|m| m.ghost_correction).collect();
+    assert!(
+        corrections.last().unwrap() < corrections.first().unwrap(),
+        "the compatibility correction must shrink with the mesh: {corrections:?}"
+    );
+    for (m, &n) in measurements.iter().zip(&resolutions) {
+        assert!(
+            m.skipped_samples == 0,
+            "surface-force reconstruction skipped {} samples at n = {n}",
+            m.skipped_samples
+        );
+    }
+    // Both routes read a first-order-accurate solution, so their errors fall
+    // at first order: measured surface 0.52 / 0.29 / 0.15 at 16 / 32 / 64
+    // (halving each refinement). Monotone convergence is the claim; a
+    // sub-5% load needs a finer grid than this suite runs.
+    assert!(
+        surface_errors.windows(2).all(|w| w[1] < w[0]) && cv_errors.windows(2).all(|w| w[1] < w[0]),
+        "both force routes must converge toward the exact force: surface {surface_errors:?}, \
+         control volume {cv_errors:?}"
+    );
+    assert!(
+        *surface_errors.last().unwrap() < 0.2 && *cv_errors.last().unwrap() < 0.2,
+        "at n = 64 both routes must be within 20% of the exact force: surface {:.3e}, \
+         control volume {:.3e}",
+        surface_errors.last().unwrap(),
+        cv_errors.last().unwrap()
+    );
+    Ok(())
+}

crates/specialized/rtx-cfd/tests/turek_hron_cfd.rs

@@ -0,0 +1,258 @@
+//! Turek–Hron CFD1: steady laminar flow (Re = 20) past the rigid cylinder
+//! with the rigid flag attached, on the embedded-boundary PISO solver.
+//!
+//! Geometry and parameters from the FEATFLOW benchmark definition (sourced
+//! 2026-08-20, see omni-cortex `docs/turek_hron_geometry_decision.md`):
+//! channel `[0, 2.5] x [0, 0.41]`, cylinder centre (0.2, 0.2) radius 0.05,
+//! flag `[0.25, 0.6] x [0.19, 0.21]`, `rho = 1000`, `nu = 1e-3`, parabolic
+//! inflow with mean `U = 0.2` (max 0.3), no-slip walls, outlet at the right.
+//! Reference (level 6): **drag 14.2929, lift 1.11905** on cylinder + flag.
+//!
+//! The flag is modelled as `[0.20, 0.6] x [0.19, 0.21]`: its left 5 cm lie
+//! inside the cylinder, which removes the two 1 mm fluid wedges the literal
+//! corners (0.25, 0.19 ± 0.01) would leave between bar and circle — below
+//! the grid scale here, and filled in the benchmark's own meshes.
+//!
+//! This is the first quantitative claim of the embedded solver against an
+//! external reference. The in-suite resolution is bounded by the dev
+//! profile's speed (the SOR projection: ~0.1 s/step at h = 10 mm, an hour
+//! per run at 5 mm); the assertion is correspondingly the measured band,
+//! not an accuracy claim, with the two load routes required to agree with
+//! each other as well. Measured at h = 10 mm (flag two cells thick), fully
+//! settled (drag stagnant to four digits): surface route drag 15.71, lift
+//! 0.936 (3 junction samples skipped); control-volume route drag 15.62,
+//! lift 1.079 — drag routes agree to 0.6%, both +9.5% on the reference.
+//! The refinement study that turns this into a claim waits on a multigrid
+//! Poisson solver (falsifier 4 of the geometry decision).
+
+use rtx_cfd::solvers::incompressible::{
+    AleBoundaries, EmbeddedBody, EmbeddedParameters, EmbeddedPisoSolver, FlowField, SideBoundary,
+};
+use rtx_cfd::{CfdConfig, CfdResult};
+
+const L: f64 = 2.5;
+const H: f64 = 0.41;
+const RHO: f64 = 1000.0;
+const NU: f64 = 1e-3;
+const U_MEAN: f64 = 0.2;
+const REF_DRAG: f64 = 14.2929;
+const REF_LIFT: f64 = 1.11905;
+
+fn inflow(y: f64) -> f64 {
+    1.5 * U_MEAN * y * (H - y) / (0.5 * H).powi(2)
+}
+
+fn body() -> EmbeddedBody {
+    EmbeddedBody::union(
+        EmbeddedBody::circle(0.2, 0.2, 0.05),
+        EmbeddedBody::rectangle(0.20, 0.19, 0.6, 0.21),
+    )
+}
+
+struct Cfd1 {
+    drag_surface: f64,
+    lift_surface: f64,
+    skipped: usize,
+    drag_cv: f64,
+    lift_cv: f64,
+    steps: usize,
+    seconds: f64,
+}
+
+/// March CFD1 to a steady state on a grid of `ny` cells across the channel.
+/// Steady means the control-volume drag has stopped moving: its relative
+/// change over the last 200 steps below `1e-4`, after at least one flow-
+/// through time — "the answer stopped moving", not a residual.
+async fn run_cfd1(ny: usize) -> CfdResult<Cfd1> {
+    let h = H / ny as f64;
+    let nx = (L / h).round() as usize;
+    let mu = RHO * NU;
+    // Explicit predictor: the COMBINED criterion — convective Courant numbers
+    // in both directions plus the diffusion number must stay below one —
+    // with the local peak velocity taken as 1.5x the inflow peak for the
+    // 24% blockage. (Taking 0.4 of the smaller single limit, as the MMS
+    // tests do, went NaN here at t ~ 7 s: both limits are active at once.)
+    let u_peak = 1.5 * 1.5 * U_MEAN;
+    let dt = 0.25 / (2.0 * u_peak / h + 4.0 * NU / (h * h));
+
+    let config = CfdConfig::new()
+        .with_density(RHO)
+        .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,
+        },
+    };
+    let mut solver = EmbeddedPisoSolver::new(config, params)?;
+    solver.set_boundary_velocity(|x, y, _| {
+        if x <= 0.0 {
+            (inflow(y), 0.0)
+        } else {
+            (0.0, 0.0)
+        }
+    });
+    solver.set_body(body());
+
+    let mut field = FlowField::new(nx, ny, h, h)?;
+    // Start from the inflow profile everywhere (the body's faces are
+    // overwritten by the mask at initialisation).
+    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)?;
+
+    // Control volume for the momentum balance: whole cells, in the fluid
+    // on its boundary, enclosing cylinder and flag.
+    let cv = (
+        (0.10 / h).round() as usize,
+        (0.75 / h).round() as usize,
+        (0.05 / h).round() as usize,
+        (0.36 / h).round() as usize,
+    );
+    let cv_force =
+        |field: &FlowField, mask: &rtx_cfd::solvers::incompressible::EmbeddedMask, dt: f64| {
+            mask.control_volume_force(
+                &field.u,
+                &field.v,
+                &field.p,
+                &field.u_old,
+                &field.v_old,
+                dt,
+                RHO,
+                mu,
+                None,
+                cv,
+            )
+        };
+
+    let start = std::time::Instant::now();
+    let flow_through = L / U_MEAN;
+    let min_steps = (flow_through / dt).ceil() as usize;
+    let mut history: Vec<f64> = Vec::new();
+    let mut steps = 0;
+    loop {
+        let result = solver.advance(&mut field, dt).await?;
+        steps += 1;
+        if steps % 50 == 0 {
+            let (fx, _) = cv_force(&field, solver.mask().unwrap(), dt);
+            history.push(fx);
+            // Diagnostics: where is the velocity largest, did the projection
+            // converge, how big was the ghost correction.
+            let (mut umax, mut at) = (0.0f64, (0usize, 0usize));
+            for j in 0..ny {
+                for i in 0..=nx {
+                    let a = field.u[(j, i)].abs();
+                    if a > umax {
+                        umax = a;
+                        at = (j, i);
+                    }
+                }
+            }
+            if steps % 250 == 0 || umax > 3.0 * 1.5 * U_MEAN || !umax.is_finite() {
+                println!(
+                    "    ny = {ny}: step {steps} t = {:.2} s  drag_cv = {fx:.4}  max|u| = {umax:.4} at (x={:.3}, y={:.3})  \
+                     projection: converged {} residual {:.2e} passes {}  ghost corr {:.2e}  [{:.0} s wall]",
+                    solver.time(),
+                    at.1 as f64 * h,
+                    (at.0 as f64 + 0.5) * h,
+                    result.solver_result.converged,
+                    result.solver_result.final_residual,
+                    result.corrector_steps_performed,
+                    result.ghost_correction,
+                    start.elapsed().as_secs_f64()
+                );
+            }
+            assert!(
+                umax.is_finite(),
+                "velocity became non-finite at step {steps}"
+            );
+            if steps >= min_steps && history.len() > 4 {
+                let now = history[history.len() - 1];
+                let then = history[history.len() - 5];
+                if ((now - then) / now).abs() < 1e-4 {
+                    break;
+                }
+            }
+        }
+        assert!(steps < 400_000, "CFD1 at ny = {ny} did not settle");
+    }
+    let seconds = start.elapsed().as_secs_f64();
+
+    let mask = solver.mask().unwrap();
+    let surface = mask.surface_force(
+        solver.body().unwrap(),
+        &field.u,
+        &field.v,
+        &field.p,
+        mu,
+        solver.time(),
+        0.5 * h,
+    );
+    let (drag_cv, lift_cv) = cv_force(&field, mask, dt);
+    Ok(Cfd1 {
+        drag_surface: surface.fx,
+        lift_surface: surface.fy,
+        skipped: surface.skipped,
+        drag_cv,
+        lift_cv,
+        steps,
+        seconds,
+    })
+}
+
+#[tokio::test]
+async fn cfd1_drag_and_lift_against_the_featflow_reference() -> CfdResult<()> {
+    // One resolution until the projection has a multigrid solver: at the
+    // SOR cost a 62-cell run takes over an hour in the test profile.
+    let resolutions = [41usize];
+    let mut results = Vec::new();
+    for &ny in &resolutions {
+        let r = run_cfd1(ny).await?;
+        println!(
+            "  ny = {ny:3} (h = {:.4})   surface: drag {:.4} lift {:.4} (skipped {})   \
+             control volume: drag {:.4} lift {:.4}   [{} steps, {:.0} s]   reference drag {REF_DRAG} lift {REF_LIFT}",
+            H / ny as f64,
+            r.drag_surface,
+            r.lift_surface,
+            r.skipped,
+            r.drag_cv,
+            r.lift_cv,
+            r.steps,
+            r.seconds
+        );
+        results.push(r);
+    }
+    let fine = results.last().unwrap();
+    let rel = |a: f64, b: f64| ((a - b) / b).abs();
+    // Both routes within the measured band of the reference (9.5% at this
+    // grid) and of each other.
+    assert!(
+        rel(fine.drag_surface, REF_DRAG) < 0.12 && rel(fine.drag_cv, REF_DRAG) < 0.12,
+        "drag: surface {:.4}, control volume {:.4}, reference {REF_DRAG}",
+        fine.drag_surface,
+        fine.drag_cv
+    );
+    assert!(
+        rel(fine.drag_surface, fine.drag_cv) < 0.05,
+        "the two drag routes disagree: surface {:.4} vs control volume {:.4}",
+        fine.drag_surface,
+        fine.drag_cv
+    );
+    assert!(
+        rel(fine.lift_surface, REF_LIFT) < 0.25 && rel(fine.lift_cv, REF_LIFT) < 0.25,
+        "lift: surface {:.4}, control volume {:.4}, reference {REF_LIFT}",
+        fine.lift_surface,
+        fine.lift_cv
+    );
+    Ok(())
+}

crates/specialized/rtx-fea/src/analysis/nonlinear_analysis.rs

@@ -28,6 +28,7 @@ use crate::analysis::{LoadSteppingStrategy, NonlinearSolverType};
 use crate::assembly::SparseMatrix;
 use crate::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
 use crate::boundary::{BoundaryCondition, BoundaryConditionSet};
+use crate::elements::total_lagrangian::{self, saint_venant_kirchhoff};
 use crate::elements::{ElementMatrixComputer, StandardFiniteElement};
 use crate::error::{AnalysisError, FeaResult};
 use crate::materials::{MaterialDatabase, reduced_constitutive};
@@ -48,6 +49,11 @@ pub struct NonlinearStaticAnalysis {
     /// (`∫ N_i f dV`) into the external force vector.
     #[allow(clippy::type_complexity)]
     body_force: Option<Box<dyn Fn(Vector3<f64>) -> Vector3<f64> + Send + Sync>>,
+    /// Use the total-Lagrangian finite-deformation formulation with a
+    /// St. Venant–Kirchhoff law built from each material's Lamé parameters
+    /// (plane strain in 2-D), instead of the small-strain path through
+    /// `reduced_constitutive`.
+    total_lagrangian: bool,
 }
 
 impl std::fmt::Debug for NonlinearStaticAnalysis {
@@ -77,9 +83,22 @@ impl NonlinearStaticAnalysis {
             progress: 0.0,
             complete: false,
             body_force: None,
+            total_lagrangian: false,
         }
     }
 
+    /// Switch to the total-Lagrangian finite-deformation formulation
+    /// (`elements::total_lagrangian`): Green–Lagrange strain, second
+    /// Piola–Kirchhoff stress from a St. Venant–Kirchhoff law with the
+    /// material's Lamé parameters, consistent material + geometric tangent.
+    /// In 2-D this is plane strain. The body force is then a dead load per
+    /// unit reference volume.
+    #[must_use]
+    pub fn with_total_lagrangian(mut self) -> Self {
+        self.total_lagrangian = true;
+        self
+    }
+
     /// Set a body force per unit volume. See the module docs.
     pub fn set_body_force<F>(&mut self, f: F)
     where
@@ -111,8 +130,6 @@ impl NonlinearStaticAnalysis {
                         element.material_id.0
                     ))
                 })?;
-            let constitutive = reduced_constitutive(material, dim)?;
-
             let node_coords: Vec<Vector3<f64>> = element
                 .nodes
                 .iter()
@@ -130,13 +147,26 @@ impl NonlinearStaticAnalysis {
                 element_displacement[local] = solution[dof];
             }
 
-            let (f_int, k_t) = ElementMatrixComputer::compute_internal_force_and_tangent(
+            let (f_int, k_t) = if self.total_lagrangian {
+                let (lambda, mu) = material.properties().lame_parameters();
+                let constitutive = saint_venant_kirchhoff(lambda, mu, dim);
+                total_lagrangian::internal_force_and_tangent(
                     &fe,
                     &node_coords,
                     &element_displacement,
                     constitutive.as_ref(),
                     None,
-            )?;
+                )?
+            } else {
+                let constitutive = reduced_constitutive(material, dim)?;
+                ElementMatrixComputer::compute_internal_force_and_tangent(
+                    &fe,
+                    &node_coords,
+                    &element_displacement,
+                    constitutive.as_ref(),
+                    None,
+                )?
+            };
 
             for (local_row, &dof_row) in dofs.iter().enumerate() {
                 let Some(free_row) = free_index[dof_row] else {

crates/specialized/rtx-fea/src/elements/mod.rs

@@ -15,6 +15,8 @@ pub mod isoparametric;
 pub mod jacobian;
 pub mod quadrature;
 pub mod shape_functions;
+/// Total-Lagrangian finite-deformation element routines
+pub mod total_lagrangian;
 
 use crate::error::{ElementError, FeaResult};
 use crate::mesh::ElementType;

crates/specialized/rtx-fea/src/elements/total_lagrangian.rs

@@ -0,0 +1,196 @@
+//! Total-Lagrangian internal force and consistent tangent for finite
+//! deformation, and the St. Venant–Kirchhoff constitutive closure.
+//!
+//! Everything is integrated on the *reference* configuration: the element's
+//! node coordinates are the undeformed ones, shape-function derivatives are
+//! `dN/dX`, and the stress measure is the second Piola–Kirchhoff tensor `S`
+//! work-conjugate to the Green–Lagrange strain `E = (FᵀF − I)/2`.
+//!
+//! Per quadrature point, with `H = ∂u/∂X`, `F = I + H`:
+//!
+//! ```text
+//!   f_int = ∫ B_Lᵀ S dV,     K_T = ∫ B_Lᵀ C B_L dV + ∫ (∇N_a · S ∇N_b) I dV
+//! ```
+//!
+//! where `B_L` is the linearised strain–displacement operator of the
+//! *current* deformation (`δE = sym(Fᵀ δH)`), `C = ∂S/∂E` the material
+//! tangent, and the second term the geometric (initial-stress) stiffness.
+//! Voigt convention is the element library's: 2-D `[E11, E22, 2E12]`, 3-D
+//! `[E11, E22, E33, 2E23, 2E13, 2E12]` — engineering shear, so `B_Lᵀ S` is
+//! the virtual work `S11 δE11 + S22 δE22 + 2 S12 δE12` with no extra factor.
+//!
+//! In 2-D the formulation is **plane strain** (`F33 = 1`, `E33 = 0`), which
+//! is what a 2-D St. Venant–Kirchhoff continuum means; the Turek–Hron CSM
+//! and FSI benchmarks are defined this way.
+//!
+//! What the tests pin: at zero displacement the tangent equals the
+//! small-strain plane-strain stiffness; the tangent is the exact derivative
+//! of the internal force (finite differences); a finite rigid rotation
+//! produces no internal force (the small-strain routine fails this — its
+//! negative control); manufactured solutions at the element orders; and
+//! the Turek–Hron CSM1/CSM2 deflections.
+
+use super::FiniteElement;
+use crate::error::{ElementError, FeaResult};
+use nalgebra::{DMatrix, DVector, Vector3};
+
+/// Constitutive closure on Green–Lagrange strain in Voigt form, returning
+/// the second Piola–Kirchhoff stress (Voigt) and the tangent `∂S/∂E`.
+pub type HyperelasticClosure =
+    Box<dyn Fn(&DVector<f64>) -> FeaResult<(DVector<f64>, DMatrix<f64>)>>;
+
+/// St. Venant–Kirchhoff: `S = λ tr(E) I + 2μ E`, with the constant tangent
+/// `C = λ 1⊗1 + 2μ I` (shear diagonal `μ` in engineering-shear Voigt).
+/// `spatial_dim` 2 is plane strain.
+pub fn saint_venant_kirchhoff(lambda: f64, mu: f64, spatial_dim: usize) -> HyperelasticClosure {
+    let n = if spatial_dim == 2 { 3 } else { 6 };
+    let normal = spatial_dim;
+    let mut c = DMatrix::zeros(n, n);
+    for i in 0..normal {
+        for j in 0..normal {
+            c[(i, j)] = lambda;
+        }
+        c[(i, i)] += 2.0 * mu;
+    }
+    for i in normal..n {
+        c[(i, i)] = mu;
+    }
+    Box::new(move |e: &DVector<f64>| {
+        if e.len() != n {
+            return Err(ElementError::MatrixComputationFailed {
+                reason: format!("strain has {} components, expected {n}", e.len()),
+            }
+            .into());
+        }
+        Ok((&c * e, c.clone()))
+    })
+}
+
+/// Total-Lagrangian element internal force and consistent tangent at the
+/// given element displacement (interleaved per node `[u0, v0, (w0), u1, …]`).
+pub fn internal_force_and_tangent(
+    element: &dyn FiniteElement,
+    node_coords: &[Vector3<f64>],
+    element_displacement: &DVector<f64>,
+    constitutive: &dyn Fn(&DVector<f64>) -> FeaResult<(DVector<f64>, DMatrix<f64>)>,
+    quadrature_order: Option<usize>,
+) -> FeaResult<(DVector<f64>, DMatrix<f64>)> {
+    let quad_rule = element.quadrature_rule(quadrature_order)?;
+    let dim = element.spatial_dimension();
+    let num_nodes = element.num_nodes();
+    let total_dofs = num_nodes * dim;
+    if element_displacement.len() != total_dofs {
+        return Err(ElementError::MatrixComputationFailed {
+            reason: format!(
+                "element displacement has {} entries, element has {total_dofs} DOFs",
+                element_displacement.len()
+            ),
+        }
+        .into());
+    }
+    if dim != 2 && dim != 3 {
+        return Err(ElementError::MatrixComputationFailed {
+            reason: format!("total-Lagrangian formulation needs 2-D or 3-D, got {dim}"),
+        }
+        .into());
+    }
+    let n_strain = if dim == 2 { 3 } else { 6 };
+    // Voigt index pairs: normals first, then shears in the library's order.
+    let pairs: Vec<(usize, usize)> = if dim == 2 {
+        vec![(0, 0), (1, 1), (0, 1)]
+    } else {
+        vec![(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]
+    };
+
+    let mut internal_force = DVector::zeros(total_dofs);
+    let mut tangent = DMatrix::zeros(total_dofs, total_dofs);
+
+    for point in &quad_rule.points {
+        let shape_eval = element.shape_functions(&point.coords)?;
+        let jacobian_eval = element.jacobian(&point.coords, node_coords)?;
+        if !jacobian_eval.is_valid() {
+            return Err(ElementError::JacobianSingular {
+                det: jacobian_eval.determinant,
+            }
+            .into());
+        }
+        // dN_a/dX_j, num_nodes x dim.
+        let g = jacobian_eval.transform_derivatives(&shape_eval.derivatives)?;
+
+        // Deformation gradient F = I + sum_a u_a ⊗ ∇N_a.
+        let mut f: DMatrix<f64> = DMatrix::identity(dim, dim);
+        for a in 0..num_nodes {
+            for i in 0..dim {
+                let ui = element_displacement[a * dim + i];
+                for j in 0..dim {
+                    f[(i, j)] += ui * g[(a, j)];
+                }
+            }
+        }
+        // Green–Lagrange strain, Voigt with engineering shear.
+        let c_right = f.transpose() * &f;
+        let mut e_voigt: DVector<f64> = DVector::zeros(n_strain);
+        for (k, &(i, j)) in pairs.iter().enumerate() {
+            let e_ij = 0.5 * (c_right[(i, j)] - if i == j { 1.0 } else { 0.0 });
+            e_voigt[k] = if i == j { e_ij } else { 2.0 * e_ij };
+        }
+        let (s_voigt, c_mat) = constitutive(&e_voigt)?;
+        if s_voigt.len() != n_strain || c_mat.nrows() != n_strain || c_mat.ncols() != n_strain {
+            return Err(ElementError::MatrixComputationFailed {
+                reason: format!(
+                    "constitutive closure returned stress of {} and tangent {}x{}, \
+                     expected {n_strain} components",
+                    s_voigt.len(),
+                    c_mat.nrows(),
+                    c_mat.ncols()
+                ),
+            }
+            .into());
+        }
+        // PK2 as a tensor.
+        let mut s: DMatrix<f64> = DMatrix::zeros(dim, dim);
+        for (k, &(i, j)) in pairs.iter().enumerate() {
+            s[(i, j)] = s_voigt[k];
+            s[(j, i)] = s_voigt[k];
+        }
+
+        // B_L: row k (pair i,j), column (a, m):
+        //   i == j : F[m][i] g[a][i]
+        //   i != j : F[m][i] g[a][j] + F[m][j] g[a][i]
+        let mut b_l: DMatrix<f64> = DMatrix::zeros(n_strain, total_dofs);
+        for (k, &(i, j)) in pairs.iter().enumerate() {
+            for a in 0..num_nodes {
+                for m in 0..dim {
+                    let col = a * dim + m;
+                    b_l[(k, col)] = if i == j {
+                        f[(m, i)] * g[(a, i)]
+                    } else {
+                        f[(m, i)] * g[(a, j)] + f[(m, j)] * g[(a, i)]
+                    };
+                }
+            }
+        }
+
+        let scale = jacobian_eval.determinant().abs() * point.weight;
+        internal_force += b_l.transpose() * &s_voigt * scale;
+        tangent += b_l.transpose() * &c_mat * &b_l * scale;
+
+        // Geometric stiffness: (∇N_a · S ∇N_b) on each diagonal block.
+        for a in 0..num_nodes {
+            for b in 0..num_nodes {
+                let mut gsg = 0.0;
+                for i in 0..dim {
+                    for j in 0..dim {
+                        gsg += g[(a, i)] * s[(i, j)] * g[(b, j)];
+                    }
+                }
+                let value = gsg * scale;
+                for m in 0..dim {
+                    tangent[(a * dim + m, b * dim + m)] += value;
+                }
+            }
+        }
+    }
+
+    Ok((internal_force, tangent))
+}

crates/specialized/rtx-fea/tests/total_lagrangian_svk.rs

@@ -0,0 +1,553 @@
+//! Verification of the total-Lagrangian St. Venant–Kirchhoff path
+//! (`elements::total_lagrangian`, `NonlinearStaticAnalysis::with_total_lagrangian`).
+//!
+//! In the order the claims must be established:
+//!
+//! 1. At zero displacement the TL tangent is the small-strain plane-strain
+//!    stiffness, to rounding.
+//! 2. The tangent is the derivative of the internal force — finite
+//!    differences at a finite random displacement, for Quad4, Quad8, Hex8.
+//!    A plausible-but-wrong geometric stiffness fails this and nothing else.
+//! 3. A finite rigid rotation produces no internal force (Green–Lagrange
+//!    strain is objective); the small-strain routine does produce one — the
+//!    negative control that shows the test has teeth.
+//! 4. Manufactured solutions at finite strain recover the element orders
+//!    (Quad4 ~2, Quad8 ~3 in L2), with the body force obtained by finite
+//!    differences of the exact first Piola–Kirchhoff stress.
+//! 5. Turek–Hron CSM1 and CSM2 (flag under gravity, clamped at the cylinder):
+//!    reference `u_x(A) = −7.18777e-3, u_y(A) = −66.1029e-3` (CSM1) and
+//!    `−0.469006e-3, −16.9740e-3` (CSM2), FEATFLOW tables.
+
+use nalgebra::{DMatrix, DVector, Vector3};
+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::elements::total_lagrangian::{internal_force_and_tangent, saint_venant_kirchhoff};
+use rtx_fea::elements::{ElementMatrixComputer, FiniteElement, StandardFiniteElement};
+use rtx_fea::materials::{LinearElastic, MaterialDatabase};
+use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+use std::f64::consts::PI;
+
+const E_MOD: f64 = 1.4e6;
+const NU: f64 = 0.4;
+
+fn lame() -> (f64, f64) {
+    let mu = E_MOD / (2.0 * (1.0 + NU));
+    let lambda = E_MOD * NU / ((1.0 + NU) * (1.0 - 2.0 * NU));
+    (lambda, mu)
+}
+
+fn plane_strain_d() -> DMatrix<f64> {
+    let (lambda, mu) = lame();
+    let mut d = DMatrix::zeros(3, 3);
+    d[(0, 0)] = lambda + 2.0 * mu;
+    d[(1, 1)] = lambda + 2.0 * mu;
+    d[(0, 1)] = lambda;
+    d[(1, 0)] = lambda;
+    d[(2, 2)] = mu;
+    d
+}
+
+fn quad4_coords() -> Vec<Vector3<f64>> {
+    // A deliberately non-rectangular quadrilateral.
+    vec![
+        Vector3::new(0.0, 0.0, 0.0),
+        Vector3::new(1.1, 0.1, 0.0),
+        Vector3::new(0.9, 1.0, 0.0),
+        Vector3::new(-0.1, 0.8, 0.0),
+    ]
+}
+
+fn quad8_coords() -> Vec<Vector3<f64>> {
+    let c = quad4_coords();
+    let mid = |a: usize, b: usize| 0.5 * (c[a] + c[b]);
+    vec![
+        c[0],
+        c[1],
+        c[2],
+        c[3],
+        mid(0, 1),
+        mid(1, 2),
+        mid(2, 3),
+        mid(3, 0),
+    ]
+}
+
+fn hex8_coords() -> Vec<Vector3<f64>> {
+    vec![
+        Vector3::new(0.0, 0.0, 0.0),
+        Vector3::new(1.0, 0.0, 0.1),
+        Vector3::new(1.1, 1.0, 0.0),
+        Vector3::new(0.0, 0.9, 0.0),
+        Vector3::new(0.1, 0.0, 1.0),
+        Vector3::new(1.0, 0.1, 1.0),
+        Vector3::new(1.0, 1.0, 1.1),
+        Vector3::new(0.0, 1.0, 0.9),
+    ]
+}
+
+/// Deterministic pseudo-random displacement of amplitude `amp`.
+fn pseudo_random(n: usize, amp: f64, seed: u64) -> DVector<f64> {
+    let mut x = seed;
+    DVector::from_iterator(
+        n,
+        (0..n).map(|_| {
+            x = x
+                .wrapping_mul(6364136223846793005)
+                .wrapping_add(1442695040888963407);
+            amp * (((x >> 33) as f64) / (1u64 << 31) as f64 - 1.0)
+        }),
+    )
+}
+
+/// 1. Zero displacement: TL tangent == small-strain plane-strain stiffness.
+#[test]
+fn at_zero_displacement_the_tangent_is_the_plane_strain_stiffness() {
+    let coords = quad4_coords();
+    let fe = StandardFiniteElement::new(ElementType::Quad4, coords.clone());
+    let (lambda, mu) = lame();
+    let svk = saint_venant_kirchhoff(lambda, mu, 2);
+    let u0 = DVector::zeros(8);
+    let (f_int, k_tl) = internal_force_and_tangent(&fe, &coords, &u0, svk.as_ref(), None).unwrap();
+    let d = plane_strain_d();
+    let linear = move |strain: &DVector<f64>| Ok((&d * strain, d.clone()));
+    let (_, k_lin) =
+        ElementMatrixComputer::compute_internal_force_and_tangent(&fe, &coords, &u0, &linear, None)
+            .unwrap();
+    assert!(
+        f_int.norm() == 0.0,
+        "internal force at zero displacement: {}",
+        f_int.norm()
+    );
+    let diff = (&k_tl - &k_lin).norm() / k_lin.norm();
+    assert!(
+        diff < 1e-13,
+        "TL tangent at u = 0 differs from the plane-strain stiffness by {diff:.3e}"
+    );
+}
+
+/// 2. The tangent is the derivative of the internal force (central
+/// differences at a finite displacement; step chosen so the FD error is far
+/// below the tolerance for a smooth polynomial residual).
+#[test]
+fn tangent_is_the_derivative_of_the_internal_force() {
+    let (lambda, mu) = lame();
+    for (name, element_type, coords, dim) in [
+        ("Quad4", ElementType::Quad4, quad4_coords(), 2usize),
+        ("Quad8", ElementType::Quad8, quad8_coords(), 2),
+        ("Hex8", ElementType::Hex8, hex8_coords(), 3),
+    ] {
+        let fe = StandardFiniteElement::new(element_type, coords.clone());
+        let svk = saint_venant_kirchhoff(lambda, mu, dim);
+        let n = coords.len() * dim;
+        // 20% strain-level displacement: genuinely finite.
+        let u = pseudo_random(n, 0.2, 7);
+        let (_, k) = internal_force_and_tangent(&fe, &coords, &u, svk.as_ref(), None).unwrap();
+        let eps = 1e-6;
+        let mut k_fd = DMatrix::zeros(n, n);
+        for j in 0..n {
+            let mut up = u.clone();
+            let mut um = u.clone();
+            up[j] += eps;
+            um[j] -= eps;
+            let (fp, _) =
+                internal_force_and_tangent(&fe, &coords, &up, svk.as_ref(), None).unwrap();
+            let (fm, _) =
+                internal_force_and_tangent(&fe, &coords, &um, svk.as_ref(), None).unwrap();
+            let col = (fp - fm) / (2.0 * eps);
+            k_fd.set_column(j, &col);
+        }
+        let rel = (&k - &k_fd).norm() / k.norm();
+        assert!(
+            rel < 1e-7,
+            "{name}: tangent vs finite-difference derivative of f_int: relative {rel:.3e}"
+        );
+        let asym = (&k - k.transpose()).norm() / k.norm();
+        assert!(asym < 1e-12, "{name}: tangent not symmetric: {asym:.3e}");
+    }
+}
+
+/// 3. Finite rigid rotation: no internal force in the TL routine; the
+/// small-strain routine sees strain and pushes back (negative control).
+#[test]
+fn rigid_rotation_produces_no_internal_force_and_the_small_strain_routine_fails_this() {
+    let coords = quad8_coords();
+    let fe = StandardFiniteElement::new(ElementType::Quad8, coords.clone());
+    let theta: f64 = 0.6; // 34 degrees
+    let (s, c) = theta.sin_cos();
+    let mut u = DVector::zeros(16);
+    for (a, x) in coords.iter().enumerate() {
+        u[2 * a] = c * x.x - s * x.y - x.x;
+        u[2 * a + 1] = s * x.x + c * x.y - x.y;
+    }
+    let (lambda, mu) = lame();
+    let svk = saint_venant_kirchhoff(lambda, mu, 2);
+    let (f_tl, _) = internal_force_and_tangent(&fe, &coords, &u, svk.as_ref(), None).unwrap();
+    let d = plane_strain_d();
+    let linear = move |strain: &DVector<f64>| Ok((&d * strain, d.clone()));
+    let (f_small, _) =
+        ElementMatrixComputer::compute_internal_force_and_tangent(&fe, &coords, &u, &linear, None)
+            .unwrap();
+    let scale = E_MOD * u.norm();
+    assert!(
+        f_tl.norm() < 1e-10 * scale,
+        "TL internal force under a rigid rotation: {:.3e} (scale {scale:.3e})",
+        f_tl.norm()
+    );
+    assert!(
+        f_small.norm() > 1e-2 * scale,
+        "the small-strain routine should see a rigid rotation as strain; got {:.3e}",
+        f_small.norm()
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Manufactured solution at finite strain
+// ---------------------------------------------------------------------------
+
+const AMP: f64 = 0.03;
+
+fn u_exact(p: Vector3<f64>) -> Vector3<f64> {
+    let (x, y) = (p.x, p.y);
+    Vector3::new(
+        AMP * (PI * x).sin() * (PI * y).sin(),
+        AMP * (PI * x).sin() * (PI * y).sin() * 0.5 + AMP * 0.3 * x * y,
+        0.0,
+    )
+}
+
+fn grad_u(p: Vector3<f64>) -> DMatrix<f64> {
+    let (x, y) = (p.x, p.y);
+    let sx = (PI * x).sin();
+    let cx = (PI * x).cos();
+    let sy = (PI * y).sin();
+    let cy = (PI * y).cos();
+    let mut h = DMatrix::zeros(2, 2);
+    h[(0, 0)] = AMP * PI * cx * sy;
+    h[(0, 1)] = AMP * PI * sx * cy;
+    h[(1, 0)] = 0.5 * AMP * PI * cx * sy + AMP * 0.3 * y;
+    h[(1, 1)] = 0.5 * AMP * PI * sx * cy + AMP * 0.3 * x;
+    h
+}
+
+/// First Piola–Kirchhoff stress `P = F S` of the manufactured field.
+fn piola(p: Vector3<f64>) -> DMatrix<f64> {
+    let (lambda, mu) = lame();
+    let f = DMatrix::identity(2, 2) + grad_u(p);
+    let e = 0.5 * (f.transpose() * &f - DMatrix::identity(2, 2));
+    let s = lambda * e.trace() * DMatrix::identity(2, 2) + 2.0 * mu * &e;
+    f * s
+}
+
+/// Body force per unit reference volume `b = −Div P`, by central
+/// differences of the analytic `P` (step 1e-6: FD error ~1e-12 relative).
+fn body_force(p: Vector3<f64>) -> Vector3<f64> {
+    let eps = 1e-6;
+    let mut div = Vector3::zeros();
+    for j in 0..2 {
+        let mut dp = Vector3::zeros();
+        dp[j] = eps;
+        let plus = piola(p + dp);
+        let minus = piola(p - dp);
+        for i in 0..2 {
+            div[i] += (plus[(i, j)] - minus[(i, j)]) / (2.0 * eps);
+        }
+    }
+    -div
+}
+
+fn on_unit_square_boundary(p: Vector3<f64>) -> bool {
+    (0..2).any(|d| p[d].abs() < 1e-12 || (p[d] - 1.0).abs() < 1e-12)
+}
+
+fn quad4_mesh(n: usize) -> Mesh {
+    let mut mesh = Mesh::new(2).unwrap();
+    let mut grid = vec![vec![NodeId(0); n + 1]; n + 1];
+    for (i, column) in grid.iter_mut().enumerate() {
+        for (j, slot) in column.iter_mut().enumerate() {
+            *slot = mesh.add_node(Node::new_2d(i as f64 / n as f64, j as f64 / n as f64));
+        }
+    }
+    for i in 0..n {
+        for j in 0..n {
+            let nodes = vec![
+                grid[i][j],
+                grid[i + 1][j],
+                grid[i + 1][j + 1],
+                grid[i][j + 1],
+            ];
+            mesh.add_element(Element::new(ElementType::Quad4, nodes, MaterialId(0)).unwrap())
+                .unwrap();
+        }
+    }
+    mesh
+}
+
+/// `nx` by `ny` Quad8 mesh of `[x0, x1] x [y0, y1]` (serendipity lattice with
+/// cell centres left out), corners counter-clockwise then mid-edges.
+fn quad8_rect_mesh(x0: f64, x1: f64, y0: f64, y1: f64, 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 = x0 + (x1 - x0) * i as f64 / (2 * nx) as f64;
+            let y = y0 + (y1 - 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
+}
+
+fn materials() -> MaterialDatabase {
+    let mut db = MaterialDatabase::new();
+    db.add_material(MaterialId(0), LinearElastic::new(E_MOD, NU), None);
+    db
+}
+
+fn dirichlet(
+    nodes: Vec<NodeId>,
+    component: DofComponent,
+    f: impl Fn(Vector3<f64>) -> f64 + Send + Sync + 'static,
+) -> BoundaryCondition {
+    BoundaryCondition::Dirichlet(DirichletBC {
+        nodes,
+        components: vec![component],
+        condition_type: DirichletType::Spatial(SpatialFunction(Box::new(move |p| f(*p)))),
+        time_range: None,
+        ramping_factor: 1.0,
+        gradual_enforcement: false,
+    })
+}
+
+fn exact_boundary_conditions(mesh: &Mesh) -> BoundaryConditionSet {
+    let boundary: Vec<NodeId> = mesh
+        .nodes
+        .iter()
+        .filter(|(_, node)| on_unit_square_boundary(node.position()))
+        .map(|(&id, _)| id)
+        .collect();
+    let mut set = BoundaryConditionSet::new();
+    set.add_condition(dirichlet(
+        boundary.clone(),
+        DofComponent::DisplacementX,
+        |p| u_exact(p).x,
+    ));
+    set.add_condition(dirichlet(boundary, DofComponent::DisplacementY, |p| {
+        u_exact(p).y
+    }));
+    set
+}
+
+/// Quadrature-integrated L2 error against the manufactured field.
+fn l2_error(mesh: &Mesh, dof_numbering: &AdvancedDofNumbering, solution: &DVector<f64>) -> f64 {
+    let mut squared = 0.0;
+    for element in mesh.elements.values() {
+        let coords: Vec<Vector3<f64>> = element
+            .nodes
+            .iter()
+            .map(|id| mesh.get_node(*id).unwrap().position())
+            .collect();
+        let fe = StandardFiniteElement::new(element.element_type, coords.clone());
+        let rule = fe.quadrature_rule(Some(4)).unwrap();
+        let dofs: Vec<usize> = element
+            .nodes
+            .iter()
+            .flat_map(|node| dof_numbering.get_node_dofs(*node))
+            .collect();
+        for point in &rule.points {
+            let shape = fe.shape_functions(&point.coords).unwrap();
+            let jac = fe.jacobian(&point.coords, &coords).unwrap();
+            let physical = fe.map_to_physical(&point.coords, &coords).unwrap();
+            let mut uh = Vector3::zeros();
+            for a in 0..element.nodes.len() {
+                let n = shape.value(a).unwrap();
+                uh.x += n * solution[dofs[2 * a]];
+                uh.y += n * solution[dofs[2 * a + 1]];
+            }
+            let exact = u_exact(physical.coords);
+            squared += (uh - exact).norm_squared() * point.weight * jac.determinant().abs();
+        }
+    }
+    squared.sqrt()
+}
+
+fn solve_mms(mesh: Mesh) -> f64 {
+    let dof_numbering =
+        AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
+    let bcs = exact_boundary_conditions(&mesh);
+    let config = NonlinearConfig {
+        max_load_steps: 5,
+        ..NonlinearConfig::default()
+    };
+    let mut analysis = NonlinearStaticAnalysis::new(
+        mesh.clone(),
+        materials(),
+        bcs,
+        config,
+        AnalysisConfig::default(),
+    )
+    .with_total_lagrangian();
+    analysis.set_body_force(body_force);
+    let results = analysis.run().unwrap();
+    assert!(
+        results.convergence.converged,
+        "Newton did not converge on the manufactured problem"
+    );
+    l2_error(&mesh, &dof_numbering, &results.displacements)
+}
+
+fn observed_order(errors: &[f64]) -> Vec<f64> {
+    errors.windows(2).map(|w| (w[0] / w[1]).log2()).collect()
+}
+
+/// 4. Manufactured finite-strain solution: Quad4 at order 2, Quad8 at 3.
+#[test]
+fn manufactured_finite_strain_solution_converges_at_the_element_orders() {
+    // Sanity on the manufactured data: max |grad u| ~ AMP*pi ~ 0.25 — finite.
+    let quad4_errors: Vec<f64> = [4usize, 8, 16]
+        .iter()
+        .map(|&n| solve_mms(quad4_mesh(n)))
+        .collect();
+    let quad8_errors: Vec<f64> = [2usize, 4, 8, 16]
+        .iter()
+        .map(|&n| solve_mms(quad8_rect_mesh(0.0, 1.0, 0.0, 1.0, n, n)))
+        .collect();
+    let o4 = observed_order(&quad4_errors);
+    let o8 = observed_order(&quad8_errors);
+    println!("  Quad4 L2 errors {quad4_errors:?} orders {o4:?}");
+    println!("  Quad8 L2 errors {quad8_errors:?} orders {o8:?}");
+    assert!(o4.last().unwrap() > &1.8, "Quad4 order {o4:?}");
+    assert!(o8.last().unwrap() > &2.7, "Quad8 order {o8:?}");
+}
+
+// ---------------------------------------------------------------------------
+// Turek–Hron CSM1 / CSM2
+// ---------------------------------------------------------------------------
+
+struct Csm {
+    ux_a: f64,
+    uy_a: f64,
+    iterations: usize,
+}
+
+/// The flag `[0.25, 0.6] x [0.19, 0.21]`, clamped at `x = 0.25`, under
+/// gravity `g = 2` downward, density 1000, plane-strain SVK with the given
+/// shear modulus and `nu = 0.4`. Returns the displacement of
+/// `A = (0.6, 0.2)`.
+fn run_csm(mu_s: f64, nx: usize, ny: usize, load_steps: usize) -> Csm {
+    let e_mod = 2.0 * mu_s * (1.0 + NU);
+    let mesh = quad8_rect_mesh(0.25, 0.6, 0.19, 0.21, nx, ny);
+    let clamped: Vec<NodeId> = mesh
+        .nodes
+        .iter()
+        .filter(|(_, node)| (node.position().x - 0.25).abs() < 1e-12)
+        .map(|(&id, _)| id)
+        .collect();
+    let point_a = mesh
+        .nodes
+        .iter()
+        .find(|(_, node)| {
+            (node.position().x - 0.6).abs() < 1e-12 && (node.position().y - 0.2).abs() < 1e-12
+        })
+        .map(|(&id, _)| id)
+        .expect("point A (0.6, 0.2) must be a mesh node");
+    let mut bcs = BoundaryConditionSet::new();
+    bcs.add_condition(dirichlet(
+        clamped.clone(),
+        DofComponent::DisplacementX,
+        |_| 0.0,
+    ));
+    bcs.add_condition(dirichlet(clamped, DofComponent::DisplacementY, |_| 0.0));
+
+    let mut db = MaterialDatabase::new();
+    db.add_material(MaterialId(0), LinearElastic::new(e_mod, NU), None);
+    let dof_numbering =
+        AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
+    let config = NonlinearConfig {
+        max_load_steps: load_steps,
+        ..NonlinearConfig::default()
+    };
+    let mut analysis =
+        NonlinearStaticAnalysis::new(mesh, db, bcs, config, AnalysisConfig::default())
+            .with_total_lagrangian();
+    analysis.set_body_force(|_| Vector3::new(0.0, -1000.0 * 2.0, 0.0));
+    let results = analysis.run().unwrap();
+    assert!(results.convergence.converged, "CSM Newton did not converge");
+    let dofs = dof_numbering.get_node_dofs(point_a);
+    Csm {
+        ux_a: results.displacements[dofs[0]],
+        uy_a: results.displacements[dofs[1]],
+        iterations: results.convergence.iterations,
+    }
+}
+
+#[test]
+fn turek_hron_csm1_and_csm2_deflections() {
+    // (mu_s, reference ux, reference uy) in metres.
+    let cases = [
+        ("CSM1", 0.5e6, -7.18777e-3, -66.1029e-3),
+        ("CSM2", 2.0e6, -0.469006e-3, -16.9740e-3),
+    ];
+    // Measured (Quad8, plane-strain SVK, 5 load steps): CSM1 35x2
+    // (-7.006e-3, -65.14e-3), 70x4 (-7.060e-3, -65.43e-3) — converging
+    // from below onto the reference (-7.188e-3, -66.10e-3), 1.0% short in
+    // u_y and 1.8% in u_x at 5 mm elements. `TL_FINE=1` runs 140x8 as well
+    // (too slow for the suite) for the convergence record.
+    let fine_mesh = if std::env::var("TL_FINE").is_ok() {
+        (140, 8)
+    } else {
+        (70, 4)
+    };
+    for (name, mu_s, ref_ux, ref_uy) in cases {
+        let coarse = run_csm(mu_s, 35, 2, 5);
+        let fine = run_csm(mu_s, fine_mesh.0, fine_mesh.1, 5);
+        println!(
+            "  {name}: 35x2 Quad8 u(A) = ({:.5e}, {:.5e}); {}x{} Quad8 u(A) = ({:.5e}, {:.5e}) \
+             [{} Newton iterations]; reference ({ref_ux:.5e}, {ref_uy:.5e})",
+            coarse.ux_a,
+            coarse.uy_a,
+            fine_mesh.0,
+            fine_mesh.1,
+            fine.ux_a,
+            fine.uy_a,
+            fine.iterations
+        );
+        let rel = |a: f64, b: f64| ((a - b) / b).abs();
+        assert!(
+            rel(fine.uy_a, ref_uy) < 0.015,
+            "{name}: u_y(A) = {:.5e} vs reference {ref_uy:.5e}",
+            fine.uy_a
+        );
+        assert!(
+            rel(fine.ux_a, ref_ux) < 0.03,
+            "{name}: u_x(A) = {:.5e} vs reference {ref_ux:.5e}",
+            fine.ux_a
+        );
+        assert!(
+            rel(fine.uy_a, ref_uy) <= rel(coarse.uy_a, ref_uy) + 1e-4,
+            "{name}: refinement moved away from the reference"
+        );
+    }
+}