rtx-fsi: partitioned fluid-structure coupling

rtx-cfd (18,715 lines) and rtx-fea (36,576 lines) both exist and nothing
connects them -- rtx-fea is commented out of rtx-cfd's dependencies. This
is the coupling layer, and it is the piece Prof. Charbel Farhat's 2026
Guggenheim Medal citation is actually about.

It depends on NEITHER solver. The properties that make a partitioned
coupling correct -- conservation of force, moment and interface work --
are statements about the transfer operators alone, so they can be
validated now, on solvers whose canonical-benchmark validation is still
outstanding. Adapters to the concrete solvers belong above this.

TRANSFER (transfer.rs). Weights satisfy two constraints:
  sum(w_i) = 1            partition of unity  -> force conserved
  sum(w_i x_i) = x_face   linear reproduction -> MOMENT conserved

The second is the one that gets skipped. Inverse-distance weighting
satisfies the first and generally violates the second, conserving force
while corrupting moment -- which shows up as slow spurious rotation rather
than as an obvious error. Underdetermined for >4 nodes, so it takes the
minimum-norm solution w = A^T (A A^T)^+ b.

That is a PSEUDO-inverse, and not for defensiveness. A wetted surface is a
surface, so its nodes are usually planar, and for a planar patch the z
constraint row is an affine multiple of the ones row -- A A^T is genuinely
rank-deficient. The constraint is redundant there, not unsatisfiable. An
ordinary inverse rejects the most ordinary interface there is; I found
this because my first test fixture was collinear and the code correctly
refused it. Constraints are then verified against the weights actually
obtained, since a pseudo-inverse returns a least-squares answer whether or
not the system was consistent.

Motion transfer uses the TRANSPOSE of the load operator, which makes
interface work conserved identically: (Hf).v = f.(H^T v). Any other
pairing leaks energy every step, and the leak looks like physics until it
destabilises.

COUPLING (coupling.rs). Staggered and Aitken-relaxed subiteration. The
decisive tests reproduce the added-mass effect: at a gain of 2.5 the
fixed-relaxation scheme DIVERGES and is reported as CouplingDiverged
rather than as an exhausted budget, and Aitken recovers the same case. A
partitioned coupling that cannot reproduce its own classic failure mode is
not being tested hard enough. Aitken is exact for a linear fixed point, so
convergence is asserted at <=4 iterations -- pinning that this is the real
delta-squared formula and not an under-relaxation that happens to work.

SCOPE, stated up front in the crate docs: small-displacement transpiration
coupling on a fixed mesh. Deliberately not ALE and not embedded-boundary,
so the Discrete Geometric Conservation Law does not yet apply -- the mesh
does not move. Large motion needs an embedded boundary treatment; that is
the next phase, not an oversight.

External comparator named at entry: Turek-Hron FSI2/FSI3, not yet reached.

26 tests written red-first; cargo test/fmt/clippy -D warnings clean.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>

Cargo.toml

@@ -100,6 +100,7 @@ members = [
     "crates/specialized/rtx-nmf",
     "crates/specialized/rtx-fea",
     "crates/specialized/rtx-cfd",
+    "crates/specialized/rtx-fsi",
 
     # Medical imaging crates (MRI2FE parity)
     "crates/specialized/rtx-medical-core",  # Super-crate consolidating medical imaging I/O
@@ -392,6 +393,7 @@ rtx-science = { path = "crates/specialized/rtx-science", version = "1.0.0" }
 rtx-nmf = { path = "crates/specialized/rtx-nmf", version = "1.0.0" }
 rtx-fea = { path = "crates/specialized/rtx-fea", version = "1.0.0" }
 rtx-cfd = { path = "crates/specialized/rtx-cfd", version = "1.0.0" }
+rtx-fsi = { path = "crates/specialized/rtx-fsi", version = "1.0.0" }
 rtx-medical-core = { path = "crates/specialized/rtx-medical-core", version = "1.0.0" }
 rtx-medical-io = { path = "crates/specialized/rtx-medical-io", version = "1.0.0" }
 rtx-materials = { path = "crates/specialized/rtx-materials", version = "1.0.0" }

crates/specialized/rtx-fsi/Cargo.toml

@@ -0,0 +1,18 @@
+[package]
+name = "rtx-fsi"
+version = "1.0.0"
+edition.workspace = true
+rust-version = "1.92"
+authors = ["RustyTorch Team"]
+license = "MIT OR Apache-2.0"
+description = "Partitioned fluid-structure interaction coupling: conservative load and motion transfer across non-matching interfaces"
+keywords = ["fsi", "fluid-structure", "coupling", "aeroelasticity", "simulation"]
+categories = ["science", "simulation"]
+repository = "https://github.com/rustytorch/rustytorch"
+
+[dependencies]
+nalgebra = { workspace = true }
+thiserror = { workspace = true }
+
+[lints]
+workspace = true

crates/specialized/rtx-fsi/src/coupling.rs

@@ -0,0 +1,399 @@
+//! Driving the fluid and structure to agree at the interface.
+//!
+//! # Why a partitioned coupling needs to iterate at all
+//!
+//! Solving fluid and structure separately means each sees a stale version
+//! of the other. Exchanging once per time step — the *staggered* scheme —
+//! is cheap and works when the structure is much heavier than the fluid it
+//! displaces.
+//!
+//! It fails, badly, when it is not. The fluid's reaction to structural
+//! acceleration behaves like an **added mass**: the structure must
+//! accelerate the fluid around it as well as itself. Once that added mass
+//! exceeds the structural mass, the interface fixed point becomes
+//! repulsive and the staggered scheme diverges no matter how small the
+//! time step. This is the added-mass effect, and it is the reason
+//! partitioned coupling has a literature at all.
+//!
+//! Aeroelasticity in air often escapes it. A flexible structure in water,
+//! or a very light structure in air, does not.
+//!
+//! # Aitken dynamic relaxation
+//!
+//! Under-relaxing the interface update with a fixed factor can stabilise
+//! the iteration, but choosing that factor requires knowing the added-mass
+//! ratio in advance, which is exactly what nobody knows.
+//!
+//! Aitken's delta-squared method infers it from the last two residuals:
+//!
+//! ```text
+//! omega_k = -omega_{k-1} * (r_{k-1} . (r_k - r_{k-1})) / ||r_k - r_{k-1}||^2
+//! ```
+//!
+//! For a linear fixed point this recovers the exact relaxation factor, so
+//! convergence is immediate — which is both why it works and how the
+//! implementation can be tested sharply rather than by "it got there
+//! eventually".
+
+use crate::error::FsiError;
+
+/// Largest interface residual growth tolerated before declaring
+/// divergence rather than iterating to the budget.
+///
+/// A residual an order of magnitude worse than where it started is not a
+/// slow start, and reporting it as divergence names the actual failure —
+/// almost always added mass — instead of an exhausted budget.
+const DIVERGENCE_FACTOR: f64 = 10.0;
+
+/// The interface state after a successful coupling step.
+#[derive(Debug, Clone, PartialEq)]
+pub struct Converged {
+    /// The agreed interface state.
+    pub state: Vec<f64>,
+    /// Final residual norm.
+    pub residual: f64,
+    /// Iterations performed.
+    pub iterations: usize,
+}
+
+/// How the interface update is relaxed between subiterations.
+#[derive(Debug, Clone, Copy, PartialEq)]
+enum Relaxation {
+    /// A constant factor, chosen by the caller.
+    Fixed(f64),
+    /// Aitken delta-squared, inferred from successive residuals.
+    Aitken,
+}
+
+/// Fixed-point driver for the fluid–structure interface.
+#[derive(Debug, Clone, PartialEq)]
+pub struct Subiterated {
+    relaxation: Relaxation,
+    max_iterations: usize,
+    tolerance: f64,
+}
+
+impl Subiterated {
+    /// Iterate with a constant relaxation factor.
+    ///
+    /// A factor of 1.0 is the plain staggered scheme, which is the
+    /// configuration that exhibits the added-mass instability.
+    ///
+    /// # Errors
+    /// [`FsiError::InvalidParameter`] for a non-positive or non-finite
+    /// factor, a zero iteration budget, or a non-positive tolerance.
+    pub fn relaxed(factor: f64, max_iterations: usize, tolerance: f64) -> Result<Self, FsiError> {
+        if !factor.is_finite() || factor <= 0.0 {
+            return Err(FsiError::InvalidParameter {
+                parameter: "relaxation factor",
+                value: factor,
+            });
+        }
+        Self::new(Relaxation::Fixed(factor), max_iterations, tolerance)
+    }
+
+    /// Iterate with Aitken dynamic relaxation.
+    ///
+    /// # Errors
+    /// [`FsiError::InvalidParameter`] for a zero iteration budget or a
+    /// non-positive tolerance.
+    pub fn aitken(max_iterations: usize, tolerance: f64) -> Result<Self, FsiError> {
+        Self::new(Relaxation::Aitken, max_iterations, tolerance)
+    }
+
+    fn new(
+        relaxation: Relaxation,
+        max_iterations: usize,
+        tolerance: f64,
+    ) -> Result<Self, FsiError> {
+        if max_iterations == 0 {
+            return Err(FsiError::InvalidParameter {
+                parameter: "max_iterations",
+                value: 0.0,
+            });
+        }
+        if !tolerance.is_finite() || tolerance <= 0.0 {
+            return Err(FsiError::InvalidParameter {
+                parameter: "tolerance",
+                value: tolerance,
+            });
+        }
+        Ok(Self {
+            relaxation,
+            max_iterations,
+            tolerance,
+        })
+    }
+
+    /// Drive `pass` — one fluid-then-structure exchange — to a fixed
+    /// point on the interface.
+    ///
+    /// # Errors
+    /// - [`FsiError::EmptyInterface`] for an empty initial state.
+    /// - [`FsiError::CountMismatch`] if `pass` returns a different length.
+    /// - [`FsiError::NonFinite`] if `pass` returns a NaN or infinity.
+    /// - [`FsiError::CouplingDiverged`] if the residual runs away, which
+    ///   names added mass rather than blaming the budget.
+    /// - [`FsiError::CouplingNotConverged`] if the budget is exhausted.
+    pub fn solve<F>(&mut self, initial: &[f64], pass: F) -> Result<Converged, FsiError>
+    where
+        F: Fn(&[f64]) -> Vec<f64>,
+    {
+        if initial.is_empty() {
+            return Err(FsiError::EmptyInterface { side: "interface" });
+        }
+
+        let mut state = initial.to_vec();
+        let mut omega = match self.relaxation {
+            Relaxation::Fixed(factor) => factor,
+            // Aitken needs a first step to have residuals to work from.
+            Relaxation::Aitken => 1.0,
+        };
+        let mut previous_residual: Option<Vec<f64>> = None;
+        let mut first_norm = None;
+
+        for iteration in 1..=self.max_iterations {
+            let candidate = pass(&state);
+            if candidate.len() != state.len() {
+                return Err(FsiError::CountMismatch {
+                    field: "coupling pass",
+                    got: candidate.len(),
+                    expected: state.len(),
+                });
+            }
+            if let Some(index) = candidate.iter().position(|v| !v.is_finite()) {
+                return Err(FsiError::NonFinite {
+                    field: "coupling pass",
+                    index,
+                });
+            }
+
+            let residual: Vec<f64> = candidate
+                .iter()
+                .zip(&state)
+                .map(|(new, old)| new - old)
+                .collect();
+            let norm = norm_of(&residual);
+            let first = *first_norm.get_or_insert(norm);
+
+            if norm <= self.tolerance {
+                return Ok(Converged {
+                    state,
+                    residual: norm,
+                    iterations: iteration,
+                });
+            }
+            if norm > first * DIVERGENCE_FACTOR && first > 0.0 {
+                return Err(FsiError::CouplingDiverged {
+                    iterations: iteration,
+                    residual: norm,
+                });
+            }
+
+            if let Relaxation::Aitken = self.relaxation {
+                if let Some(previous) = &previous_residual {
+                    omega = aitken_factor(omega, previous, &residual);
+                }
+            }
+
+            for (value, delta) in state.iter_mut().zip(&residual) {
+                *value += omega * delta;
+            }
+            previous_residual = Some(residual);
+        }
+
+        Err(FsiError::CouplingNotConverged {
+            iterations: self.max_iterations,
+            residual: previous_residual.as_deref().map_or(f64::NAN, norm_of),
+            tolerance: self.tolerance,
+        })
+    }
+}
+
+/// Aitken delta-squared relaxation factor from successive residuals.
+///
+/// Falls back to the previous factor when the residual barely moved,
+/// since the update divides by that difference.
+fn aitken_factor(previous_omega: f64, previous: &[f64], current: &[f64]) -> f64 {
+    let difference: Vec<f64> = current
+        .iter()
+        .zip(previous)
+        .map(|(now, before)| now - before)
+        .collect();
+    let denominator: f64 = difference.iter().map(|d| d * d).sum();
+    if denominator <= f64::EPSILON {
+        return previous_omega;
+    }
+    let numerator: f64 = previous
+        .iter()
+        .zip(&difference)
+        .map(|(before, delta)| before * delta)
+        .sum();
+    -previous_omega * numerator / denominator
+}
+
+fn norm_of(values: &[f64]) -> f64 {
+    values.iter().map(|v| v * v).sum::<f64>().sqrt()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// One fluid-then-structure pass, as a map on the interface state.
+    ///
+    /// A linear map with gain `-gain` stands in for the added-mass
+    /// coupling: the fluid's reaction to structural acceleration is
+    /// proportional to the fluid density, and it opposes the motion. The
+    /// gain is the added-mass ratio, and everything about partitioned
+    /// stability follows from whether it exceeds one.
+    fn added_mass(gain: f64) -> impl Fn(&[f64]) -> Vec<f64> {
+        move |state: &[f64]| state.iter().map(|x| -gain * x).collect()
+    }
+
+    // ---- the added-mass instability ----
+
+    #[test]
+    fn staggered_coupling_survives_a_light_fluid() {
+        // Added mass well below structural mass: the classic staggered
+        // scheme is fine, which is why it is used at all.
+        let mut scheme = Subiterated::relaxed(0.5, 200, 1e-10).expect("valid");
+        let converged = scheme.solve(&[1.0], added_mass(0.3)).expect("converges");
+        assert!(converged.residual < 1e-10);
+    }
+
+    #[test]
+    fn staggered_coupling_diverges_when_the_fluid_is_heavy() {
+        // THE classic partitioned-FSI failure. Once the added mass exceeds
+        // the structural mass the fixed point is repulsive, and no amount
+        // of iterating at unit relaxation recovers it. A coupling that
+        // does not reproduce this is not being tested hard enough.
+        let mut scheme = Subiterated::relaxed(1.0, 200, 1e-10).expect("valid");
+        assert!(matches!(
+            scheme.solve(&[1.0], added_mass(2.5)),
+            Err(FsiError::CouplingDiverged { .. })
+        ));
+    }
+
+    #[test]
+    fn aitken_relaxation_recovers_the_heavy_fluid_case() {
+        // The whole point of dynamic relaxation. Same gain that destroyed
+        // the fixed-relaxation scheme, now converging.
+        let mut scheme = Subiterated::aitken(200, 1e-10).expect("valid");
+        let converged = scheme.solve(&[1.0], added_mass(2.5)).expect("converges");
+        assert!(
+            converged.residual < 1e-10,
+            "residual {}",
+            converged.residual
+        );
+    }
+
+    #[test]
+    fn aitken_is_exact_on_a_linear_map() {
+        // For a linear fixed point Aitken's delta-squared finds the exact
+        // relaxation, so it should land in very few iterations. Asserting
+        // the count pins that the update is the real Aitken formula and
+        // not an under-relaxation that happens to converge.
+        let mut scheme = Subiterated::aitken(200, 1e-12).expect("valid");
+        let converged = scheme.solve(&[1.0], added_mass(2.5)).expect("converges");
+        assert!(
+            converged.iterations <= 4,
+            "expected near-immediate convergence, took {}",
+            converged.iterations
+        );
+    }
+
+    #[test]
+    fn aitken_beats_fixed_relaxation_where_both_converge() {
+        let gain = 0.8;
+        let mut fixed = Subiterated::relaxed(0.5, 500, 1e-10).expect("valid");
+        let mut aitken = Subiterated::aitken(500, 1e-10).expect("valid");
+
+        let slow = fixed.solve(&[1.0], added_mass(gain)).expect("converges");
+        let fast = aitken.solve(&[1.0], added_mass(gain)).expect("converges");
+        assert!(
+            fast.iterations < slow.iterations,
+            "aitken {} vs fixed {}",
+            fast.iterations,
+            slow.iterations
+        );
+    }
+
+    // ---- ordinary behaviour ----
+
+    #[test]
+    fn an_already_converged_interface_does_no_work() {
+        // Zero is the fixed point of the added-mass map. Starting there
+        // must terminate immediately rather than iterating pointlessly.
+        let mut scheme = Subiterated::aitken(100, 1e-10).expect("valid");
+        let converged = scheme
+            .solve(&[0.0, 0.0], added_mass(2.0))
+            .expect("converges");
+        assert_eq!(converged.iterations, 1);
+    }
+
+    #[test]
+    fn a_multi_component_interface_converges_together() {
+        let mut scheme = Subiterated::aitken(200, 1e-10).expect("valid");
+        let converged = scheme
+            .solve(&[1.0, -2.0, 0.5, 3.0], added_mass(1.8))
+            .expect("converges");
+        assert!(converged.state.iter().all(|x| x.abs() < 1e-8));
+    }
+
+    #[test]
+    fn exhausting_the_iteration_budget_is_reported_not_hidden() {
+        // A slowly converging problem cut short must say so. Returning the
+        // unconverged state as if it were converged is how a coupling
+        // silently produces plausible nonsense.
+        let mut scheme = Subiterated::relaxed(0.01, 3, 1e-12).expect("valid");
+        assert!(matches!(
+            scheme.solve(&[1.0], added_mass(0.9)),
+            Err(FsiError::CouplingNotConverged { .. })
+        ));
+    }
+
+    // ---- refusals ----
+
+    #[test]
+    fn an_invalid_relaxation_factor_is_refused() {
+        assert!(Subiterated::relaxed(0.0, 10, 1e-8).is_err());
+        assert!(Subiterated::relaxed(-0.5, 10, 1e-8).is_err());
+        assert!(Subiterated::relaxed(f64::NAN, 10, 1e-8).is_err());
+    }
+
+    #[test]
+    fn a_zero_iteration_budget_is_refused() {
+        assert!(Subiterated::aitken(0, 1e-8).is_err());
+    }
+
+    #[test]
+    fn an_invalid_tolerance_is_refused() {
+        assert!(Subiterated::aitken(10, 0.0).is_err());
+        assert!(Subiterated::aitken(10, -1e-8).is_err());
+    }
+
+    #[test]
+    fn an_empty_interface_state_is_refused() {
+        let mut scheme = Subiterated::aitken(10, 1e-8).expect("valid");
+        assert!(scheme.solve(&[], added_mass(0.5)).is_err());
+    }
+
+    #[test]
+    fn a_map_returning_the_wrong_length_is_refused() {
+        let mut scheme = Subiterated::aitken(10, 1e-8).expect("valid");
+        assert!(matches!(
+            scheme.solve(&[1.0, 2.0], |_: &[f64]| vec![0.0]),
+            Err(FsiError::CountMismatch { .. })
+        ));
+    }
+
+    #[test]
+    fn a_map_returning_non_finite_values_is_refused() {
+        let mut scheme = Subiterated::aitken(10, 1e-8).expect("valid");
+        assert!(matches!(
+            scheme.solve(&[1.0], |_: &[f64]| vec![f64::NAN]),
+            Err(FsiError::NonFinite { .. })
+        ));
+    }
+}

crates/specialized/rtx-fsi/src/error.rs

@@ -0,0 +1,93 @@
+//! Errors raised when an interface or a transfer cannot be formed.
+
+use thiserror::Error;
+
+/// Why a coupling operation was refused.
+#[derive(Debug, Clone, PartialEq, Error)]
+pub enum FsiError {
+    /// One side of the interface was empty.
+    #[error("empty interface: {side} has no entries")]
+    EmptyInterface {
+        /// Which side was empty.
+        side: &'static str,
+    },
+
+    /// Too few structure nodes to reproduce a face centroid exactly.
+    ///
+    /// Linear reproduction in three dimensions imposes four constraints,
+    /// so at least four independent nodes are needed. Falling back to a
+    /// weaker interpolation would conserve force and silently corrupt
+    /// moment, which is worse than refusing.
+    #[error("interface needs at least {needed} structure nodes for linear reproduction, got {got}")]
+    InsufficientNodes {
+        /// Nodes supplied.
+        got: usize,
+        /// Nodes required.
+        needed: usize,
+    },
+
+    /// A field had a different length from the mesh it belongs to.
+    #[error("{field}: expected {expected} values, got {got}")]
+    CountMismatch {
+        /// Which field.
+        field: &'static str,
+        /// Length supplied.
+        got: usize,
+        /// Length required.
+        expected: usize,
+    },
+
+    /// A field value was NaN or infinite.
+    #[error("non-finite value in {field} at index {index}")]
+    NonFinite {
+        /// Which field.
+        field: &'static str,
+        /// Where.
+        index: usize,
+    },
+
+    /// A parameter given to a coupling scheme was outside its valid
+    /// range.
+    #[error("invalid coupling parameter {parameter}: {value}")]
+    InvalidParameter {
+        /// Which parameter.
+        parameter: &'static str,
+        /// The value supplied.
+        value: f64,
+    },
+
+    /// The interface iteration ran away.
+    ///
+    /// In partitioned coupling this is almost always the added-mass
+    /// effect: once the fluid's added mass exceeds the structural mass the
+    /// fixed point becomes repulsive, and a fixed relaxation factor cannot
+    /// recover it. Dynamic relaxation can.
+    #[error("coupling diverged after {iterations} iterations (residual {residual})")]
+    CouplingDiverged {
+        /// Iterations completed before divergence was detected.
+        iterations: usize,
+        /// Residual at that point.
+        residual: f64,
+    },
+
+    /// The iteration budget ran out before the interface agreed.
+    #[error(
+        "coupling did not converge in {iterations} iterations (residual {residual}, tolerance {tolerance})"
+    )]
+    CouplingNotConverged {
+        /// Iterations performed.
+        iterations: usize,
+        /// Residual reached.
+        residual: f64,
+        /// Tolerance required.
+        tolerance: f64,
+    },
+
+    /// The interpolation system was singular — the selected nodes are
+    /// degenerate, for example all collinear.
+    #[error("degenerate node neighbourhood for face {face}: cannot reproduce its centroid")]
+    DegenerateNeighbourhood {
+        /// Which face.
+        face: usize,
+    },
+}

crates/specialized/rtx-fsi/src/lib.rs

@@ -0,0 +1,41 @@
+//! # `rtx-fsi` — partitioned fluid–structure interaction coupling
+//!
+//! Fluid and structure solvers exist separately in this workspace
+//! (`rtx-cfd`, `rtx-fea`) and nothing connects them. This crate is the
+//! coupling layer: it moves loads from the fluid onto the structure and
+//! motion from the structure back onto the fluid, across an interface
+//! where the two meshes do **not** match — which is the normal case and
+//! the source of most fluid–structure interaction errors.
+//!
+//! # Why the coupling is worth building separately
+//!
+//! The properties that make a partitioned coupling correct are testable
+//! **independently of whether the solvers it couples are accurate**.
+//! Conservation of force, moment and interface work are statements about
+//! the transfer operators alone. So this crate can be validated now, on
+//! solvers whose canonical-benchmark validation is still outstanding.
+//!
+//! It therefore depends on neither `rtx-cfd` nor `rtx-fea`. It takes
+//! geometry and field values, and returns field values. Adapters to the
+//! concrete solvers belong above it.
+//!
+//! # Scope, stated up front
+//!
+//! This is a **small-displacement** coupling. Interface velocity is meant
+//! to be applied to the fluid through a boundary condition on a fixed
+//! mesh — a transpiration formulation, valid while displacements are
+//! small relative to a cell. It is deliberately not arbitrary
+//! Lagrangian–Eulerian and not an embedded boundary method, so the
+//! **Discrete Geometric Conservation Law does not yet apply**: the mesh
+//! does not move. Large motion needs an embedded boundary treatment, and
+//! that is the next phase rather than an oversight.
+
+#![forbid(unsafe_code)]
+
+mod coupling;
+mod error;
+mod transfer;
+
+pub use coupling::{Converged, Subiterated};
+pub use error::FsiError;
+pub use transfer::{FluidFace, WettedSurface};

crates/specialized/rtx-fsi/src/transfer.rs

@@ -0,0 +1,544 @@
+//! Conservative load and motion transfer across a non-matching interface.
+//!
+//! # The two constraints that make it conservative
+//!
+//! Each fluid face's load is distributed onto nearby structure nodes with
+//! weights `w_i`. Two conditions decide whether the transfer conserves
+//! anything:
+//!
+//! - `sum(w_i) = 1` — **partition of unity**. Total force is preserved.
+//! - `sum(w_i * x_i) = x_face` — **linear reproduction**. The load arrives
+//!   where it left, in the weighted-average sense, so total **moment** is
+//!   preserved too, about any point.
+//!
+//! The second is the one that gets skipped. Inverse-distance weighting
+//! satisfies partition of unity and generally violates linear
+//! reproduction, which conserves force while quietly corrupting moment —
+//! a defect that shows up as a slow spurious rotation rather than as an
+//! obvious error.
+//!
+//! Both constraints leave the weights underdetermined for more than four
+//! nodes, so this takes the **minimum-norm** solution: `w = Aᵀ(AAᵀ)⁻¹ b`,
+//! where `A` stacks the ones row and the coordinate rows. Minimum norm
+//! keeps the load spread rather than concentrated on whichever node the
+//! solver happened to favour.
+//!
+//! # Why motion transfer is the transpose
+//!
+//! Given the load operator `H` mapping face loads to nodal forces, using
+//! `Hᵀ` to map nodal velocities to face velocities makes interface work
+//! conserved *identically*:
+//!
+//! ```text
+//! (H f)·v  =  f·(Hᵀ v)
+//! ```
+//!
+//! which is just the definition of the transpose. Any other pairing leaks
+//! energy across the interface every step, and the leak looks like physics
+//! until it destabilises.
+
+use nalgebra::{Matrix4, Vector3, Vector4};
+
+use crate::error::FsiError;
+
+/// Structure nodes recruited per fluid face.
+///
+/// Four is the minimum for linear reproduction in three dimensions. More
+/// spreads the load and conditions the system better; too many turn a
+/// local transfer into a global smear.
+const NEIGHBOURS: usize = 8;
+
+/// Singular values below this share of the largest are treated as zero
+/// when forming the pseudo-inverse.
+const SINGULAR_TOLERANCE: f64 = 1e-12;
+
+/// How exactly the partition-of-unity and linear-reproduction constraints
+/// must hold before a face is accepted. Tight, because these are the
+/// properties the conservation guarantees rest on.
+const CONSTRAINT_TOLERANCE: f64 = 1e-9;
+
+/// One fluid-side boundary face on the wetted surface.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct FluidFace {
+    /// Face centroid.
+    pub centroid: Vector3<f64>,
+    /// Outward unit normal, pointing into the fluid.
+    pub normal: Vector3<f64>,
+    /// Face area, used to turn a traction into a force.
+    pub area: f64,
+}
+
+/// The interface between the two meshes, with the transfer operator
+/// already formed.
+#[derive(Debug, Clone, PartialEq)]
+pub struct WettedSurface {
+    /// Per face: the recruited nodes and their weights.
+    rows: Vec<Vec<(usize, f64)>>,
+    node_count: usize,
+}
+
+impl WettedSurface {
+    /// Form the transfer operator between a fluid boundary and a
+    /// structural interface mesh.
+    ///
+    /// # Errors
+    /// - [`FsiError::EmptyInterface`] if either side is empty.
+    /// - [`FsiError::InsufficientNodes`] with fewer than four nodes.
+    /// - [`FsiError::DegenerateNeighbourhood`] if a face's nearest nodes
+    ///   are collinear or coplanar in a way that makes the constraint
+    ///   system singular.
+    pub fn build(
+        fluid_faces: &[FluidFace],
+        structure_nodes: &[Vector3<f64>],
+    ) -> Result<Self, FsiError> {
+        if fluid_faces.is_empty() {
+            return Err(FsiError::EmptyInterface { side: "fluid" });
+        }
+        if structure_nodes.is_empty() {
+            return Err(FsiError::EmptyInterface { side: "structure" });
+        }
+        if structure_nodes.len() < 4 {
+            return Err(FsiError::InsufficientNodes {
+                got: structure_nodes.len(),
+                needed: 4,
+            });
+        }
+
+        let mut rows = Vec::with_capacity(fluid_faces.len());
+        for (face_index, face) in fluid_faces.iter().enumerate() {
+            let recruited = nearest(structure_nodes, face.centroid, NEIGHBOURS);
+            let weights = solve_weights(structure_nodes, &recruited, face.centroid)
+                .ok_or(FsiError::DegenerateNeighbourhood { face: face_index })?;
+            rows.push(recruited.into_iter().zip(weights).collect());
+        }
+
+        Ok(Self {
+            rows,
+            node_count: structure_nodes.len(),
+        })
+    }
+
+    /// Number of fluid faces on the interface.
+    #[must_use]
+    pub fn face_count(&self) -> usize {
+        self.rows.len()
+    }
+
+    /// Number of structure nodes on the interface.
+    #[must_use]
+    pub const fn node_count(&self) -> usize {
+        self.node_count
+    }
+
+    /// The `(node, weight)` pairs a face distributes onto.
+    #[must_use]
+    pub fn weights_for(&self, face: usize) -> &[(usize, f64)] {
+        self.rows.get(face).map_or(&[], Vec::as_slice)
+    }
+
+    /// Fluid tractions to structural nodal forces.
+    ///
+    /// Conserves total force and total moment exactly.
+    ///
+    /// # Errors
+    /// [`FsiError::CountMismatch`] if the field lengths disagree with the
+    /// interface, or [`FsiError::NonFinite`] for a NaN or infinite value.
+    pub fn transfer_load(
+        &self,
+        fluid_faces: &[FluidFace],
+        tractions: &[Vector3<f64>],
+    ) -> Result<Vec<Vector3<f64>>, FsiError> {
+        if fluid_faces.len() != self.rows.len() {
+            return Err(FsiError::CountMismatch {
+                field: "fluid_faces",
+                got: fluid_faces.len(),
+                expected: self.rows.len(),
+            });
+        }
+        if tractions.len() != self.rows.len() {
+            return Err(FsiError::CountMismatch {
+                field: "tractions",
+                got: tractions.len(),
+                expected: self.rows.len(),
+            });
+        }
+        check_finite("tractions", tractions)?;
+
+        let mut nodal = vec![Vector3::zeros(); self.node_count];
+        for ((face, traction), row) in fluid_faces.iter().zip(tractions).zip(&self.rows) {
+            let force = traction * face.area;
+            for (node, weight) in row {
+                nodal[*node] += force * *weight;
+            }
+        }
+        Ok(nodal)
+    }
+
+    /// Structural nodal velocities to fluid face velocities.
+    ///
+    /// Uses the transpose of the load operator, which is what makes
+    /// interface work conserved identically.
+    ///
+    /// # Errors
+    /// [`FsiError::CountMismatch`] or [`FsiError::NonFinite`].
+    pub fn transfer_motion(
+        &self,
+        node_velocities: &[Vector3<f64>],
+    ) -> Result<Vec<Vector3<f64>>, FsiError> {
+        if node_velocities.len() != self.node_count {
+            return Err(FsiError::CountMismatch {
+                field: "node_velocities",
+                got: node_velocities.len(),
+                expected: self.node_count,
+            });
+        }
+        check_finite("node_velocities", node_velocities)?;
+
+        Ok(self
+            .rows
+            .iter()
+            .map(|row| {
+                row.iter()
+                    .map(|(node, weight)| node_velocities[*node] * *weight)
+                    .sum()
+            })
+            .collect())
+    }
+}
+
+/// Indices of the `count` nodes nearest a point.
+fn nearest(nodes: &[Vector3<f64>], point: Vector3<f64>, count: usize) -> Vec<usize> {
+    let mut ordered: Vec<(f64, usize)> = nodes
+        .iter()
+        .enumerate()
+        .map(|(index, node)| ((node - point).norm_squared(), index))
+        .collect();
+    ordered.sort_by(|a, b| a.0.total_cmp(&b.0));
+    ordered
+        .into_iter()
+        .take(count.min(nodes.len()))
+        .map(|(_, index)| index)
+        .collect()
+}
+
+/// Minimum-norm weights satisfying partition of unity and linear
+/// reproduction of `centroid`.
+///
+/// Solves `w = Aᵀ(AAᵀ)⁺ b` with `A` the 4×k constraint matrix, using the
+/// **pseudo-inverse** rather than an inverse.
+///
+/// That is not defensive programming, it is the common case. A wetted
+/// surface is a surface, so its nodes are usually planar — and for a
+/// planar patch the `z` constraint row is an affine multiple of the ones
+/// row, leaving `AAᵀ` genuinely rank-deficient. The constraint is not
+/// unsatisfiable there, it is *redundant*: every partition of unity
+/// reproduces a coordinate that is the same at every node. An ordinary
+/// inverse would reject the most ordinary interface there is.
+///
+/// The constraints are then checked against the weights actually
+/// obtained, because a pseudo-inverse returns a least-squares answer
+/// whether or not the system was consistent. If the target genuinely
+/// cannot be reproduced — a face outside the span of its recruited
+/// nodes — the residual reveals it and the face is refused.
+fn solve_weights(
+    nodes: &[Vector3<f64>],
+    recruited: &[usize],
+    centroid: Vector3<f64>,
+) -> Option<Vec<f64>> {
+    let mut gram = Matrix4::zeros();
+    for index in recruited {
+        let node = nodes[*index];
+        let row = Vector4::new(1.0, node.x, node.y, node.z);
+        gram += row * row.transpose();
+    }
+
+    let target = Vector4::new(1.0, centroid.x, centroid.y, centroid.z);
+    let lambda = gram.pseudo_inverse(SINGULAR_TOLERANCE).ok()? * target;
+
+    let weights: Vec<f64> = recruited
+        .iter()
+        .map(|index| {
+            let node = nodes[*index];
+            Vector4::new(1.0, node.x, node.y, node.z).dot(&lambda)
+        })
+        .collect();
+
+    // Verify what was asked for, rather than trusting the solve.
+    let unity: f64 = weights.iter().sum();
+    let reproduced: Vector3<f64> = recruited
+        .iter()
+        .zip(&weights)
+        .map(|(index, weight)| nodes[*index] * *weight)
+        .sum();
+    if (unity - 1.0).abs() > CONSTRAINT_TOLERANCE
+        || (reproduced - centroid).norm() > CONSTRAINT_TOLERANCE
+    {
+        return None;
+    }
+
+    Some(weights)
+}
+
+fn check_finite(field: &'static str, values: &[Vector3<f64>]) -> Result<(), FsiError> {
+    for (index, value) in values.iter().enumerate() {
+        if !value.iter().all(|component| component.is_finite()) {
+            return Err(FsiError::NonFinite { field, index });
+        }
+    }
+    Ok(())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use nalgebra::Vector3;
+
+    /// A structure-side patch that deliberately does not match the fluid
+    /// side: 4x4 nodes at 0.25 spacing on the z = 0 plane. Matching meshes
+    /// hide every interesting failure, and a planar patch is the ordinary
+    /// case a wetted surface presents.
+    fn structure_nodes() -> Vec<Vector3<f64>> {
+        let mut nodes = Vec::new();
+        for i in 0..4 {
+            for j in 0..4 {
+                nodes.push(Vector3::new(f64::from(i) * 0.25, f64::from(j) * 0.25, 0.0));
+            }
+        }
+        nodes
+    }
+
+    /// Fluid faces at 0.35 spacing, offset so no face sits on a node.
+    fn fluid_faces() -> Vec<FluidFace> {
+        let mut faces = Vec::new();
+        for i in 0..3 {
+            for j in 0..3 {
+                faces.push(FluidFace {
+                    centroid: Vector3::new(0.1 + f64::from(i) * 0.3, 0.1 + f64::from(j) * 0.3, 0.0),
+                    normal: Vector3::new(0.0, 0.0, 1.0),
+                    area: 0.09,
+                });
+            }
+        }
+        faces
+    }
+
+    // ---- the conservation properties ----
+
+    #[test]
+    fn the_weights_form_a_partition_of_unity() {
+        // Sum to one is what conserves total force. Without it the
+        // coupling quietly gains or loses load every step.
+        let surface = WettedSurface::build(&fluid_faces(), &structure_nodes()).expect("buildable");
+        for face in 0..surface.face_count() {
+            let total: f64 = surface.weights_for(face).iter().map(|(_, w)| w).sum();
+            assert!(
+                (total - 1.0).abs() < 1e-12,
+                "face {face} weights sum to {total}"
+            );
+        }
+    }
+
+    #[test]
+    fn the_weights_reproduce_the_face_centroid() {
+        // Linear reproduction. This is what conserves *moment*: the load
+        // must arrive at the same place it left, in the weighted-average
+        // sense. Partition of unity alone conserves force and silently
+        // corrupts the moment.
+        let nodes = structure_nodes();
+        let surface = WettedSurface::build(&fluid_faces(), &nodes).expect("buildable");
+        for (face, expected) in fluid_faces().iter().enumerate() {
+            let reproduced: Vector3<f64> = surface
+                .weights_for(face)
+                .iter()
+                .map(|(node, w)| nodes[*node] * *w)
+                .sum();
+            let error = (reproduced - expected.centroid).norm();
+            assert!(error < 1e-10, "face {face} centroid off by {error}");
+        }
+    }
+
+    #[test]
+    fn total_force_is_preserved_across_the_interface() {
+        let nodes = structure_nodes();
+        let faces = fluid_faces();
+        let surface = WettedSurface::build(&faces, &nodes).expect("buildable");
+
+        // A non-uniform traction, so a bug cannot hide behind symmetry.
+        let tractions: Vec<Vector3<f64>> = (0..faces.len())
+            .map(|i| Vector3::new(1.0 + f64::from(i as i32), -2.0, 0.5))
+            .collect();
+
+        let fluid_total: Vector3<f64> = faces
+            .iter()
+            .zip(&tractions)
+            .map(|(face, traction)| traction * face.area)
+            .sum();
+        let nodal = surface
+            .transfer_load(&faces, &tractions)
+            .expect("transferable");
+        let structure_total: Vector3<f64> = nodal.iter().sum();
+
+        let error = (structure_total - fluid_total).norm();
+        assert!(error < 1e-10, "force lost across interface: {error}");
+    }
+
+    #[test]
+    fn total_moment_is_preserved_across_the_interface() {
+        let nodes = structure_nodes();
+        let faces = fluid_faces();
+        let surface = WettedSurface::build(&faces, &nodes).expect("buildable");
+        let tractions: Vec<Vector3<f64>> = (0..faces.len())
+            .map(|i| Vector3::new(1.0 + f64::from(i as i32), -2.0, 0.5))
+            .collect();
+
+        // About an arbitrary point, deliberately not the origin: a scheme
+        // that only conserves moment about one special point is not
+        // conserving moment.
+        let about = Vector3::new(-0.7, 1.3, 2.1);
+        let fluid_moment: Vector3<f64> = faces
+            .iter()
+            .zip(&tractions)
+            .map(|(face, traction)| (face.centroid - about).cross(&(traction * face.area)))
+            .sum();
+
+        let nodal = surface
+            .transfer_load(&faces, &tractions)
+            .expect("transferable");
+        let structure_moment: Vector3<f64> = nodes
+            .iter()
+            .zip(&nodal)
+            .map(|(position, force)| (position - about).cross(force))
+            .sum();
+
+        let error = (structure_moment - fluid_moment).norm();
+        assert!(error < 1e-10, "moment lost across interface: {error}");
+    }
+
+    #[test]
+    fn interface_work_is_conserved() {
+        // The decisive one. Motion transfer uses the transpose of the load
+        // transfer, so the work the fluid does equals the work the
+        // structure receives, exactly. A coupling that fails this injects
+        // or drains energy every step and will eventually destabilise for
+        // reasons that look like physics.
+        let nodes = structure_nodes();
+        let faces = fluid_faces();
+        let surface = WettedSurface::build(&faces, &nodes).expect("buildable");
+
+        let tractions: Vec<Vector3<f64>> = (0..faces.len())
+            .map(|i| Vector3::new(0.3 * f64::from(i as i32), 1.7, -0.4))
+            .collect();
+        let node_velocities: Vec<Vector3<f64>> = (0..nodes.len())
+            .map(|i| Vector3::new(0.1, 0.2 * f64::from(i as i32), -0.05))
+            .collect();
+
+        let nodal_forces = surface
+            .transfer_load(&faces, &tractions)
+            .expect("transferable");
+        let face_velocities = surface
+            .transfer_motion(&node_velocities)
+            .expect("transferable");
+
+        let fluid_work: f64 = faces
+            .iter()
+            .zip(&tractions)
+            .zip(&face_velocities)
+            .map(|((face, traction), velocity)| (traction * face.area).dot(velocity))
+            .sum();
+        let structure_work: f64 = nodal_forces
+            .iter()
+            .zip(&node_velocities)
+            .map(|(force, velocity)| force.dot(velocity))
+            .sum();
+
+        assert!(
+            (fluid_work - structure_work).abs() < 1e-10,
+            "interface work not conserved: fluid {fluid_work}, structure {structure_work}"
+        );
+    }
+
+    // ---- limits ----
+
+    #[test]
+    fn zero_traction_moves_nothing() {
+        let nodes = structure_nodes();
+        let faces = fluid_faces();
+        let surface = WettedSurface::build(&faces, &nodes).expect("buildable");
+        let tractions = vec![Vector3::zeros(); faces.len()];
+        let nodal = surface
+            .transfer_load(&faces, &tractions)
+            .expect("transferable");
+        assert!(nodal.iter().all(|force| force.norm() < 1e-15));
+    }
+
+    #[test]
+    fn a_motionless_structure_leaves_the_fluid_boundary_still() {
+        let nodes = structure_nodes();
+        let surface = WettedSurface::build(&fluid_faces(), &nodes).expect("buildable");
+        let velocities = vec![Vector3::zeros(); nodes.len()];
+        let face_velocities = surface.transfer_motion(&velocities).expect("transferable");
+        assert!(face_velocities.iter().all(|v| v.norm() < 1e-15));
+    }
+
+    #[test]
+    fn a_rigid_translation_transfers_unchanged() {
+        // Every structure node moving together must give every fluid face
+        // the same velocity. This follows from partition of unity, and it
+        // is the check that catches a normalisation slip.
+        let nodes = structure_nodes();
+        let surface = WettedSurface::build(&fluid_faces(), &nodes).expect("buildable");
+        let rigid = Vector3::new(0.4, -1.1, 0.9);
+        let velocities = vec![rigid; nodes.len()];
+
+        for velocity in surface.transfer_motion(&velocities).expect("transferable") {
+            assert!(
+                (velocity - rigid).norm() < 1e-12,
+                "rigid translation distorted: {velocity:?}"
+            );
+        }
+    }
+
+    // ---- refusals ----
+
+    #[test]
+    fn an_interface_with_too_few_nodes_is_refused() {
+        // Linear reproduction in three dimensions needs four independent
+        // nodes. Fewer cannot satisfy the constraint, and pretending
+        // otherwise would silently break moment conservation.
+        let nodes: Vec<Vector3<f64>> = (0..3)
+            .map(|i| Vector3::new(f64::from(i), f64::from(i) * 0.5, 0.0))
+            .collect();
+        assert!(matches!(
+            WettedSurface::build(&fluid_faces(), &nodes),
+            Err(FsiError::InsufficientNodes { .. })
+        ));
+    }
+
+    #[test]
+    fn an_empty_interface_is_refused() {
+        assert!(WettedSurface::build(&[], &structure_nodes()).is_err());
+        assert!(WettedSurface::build(&fluid_faces(), &[]).is_err());
+    }
+
+    #[test]
+    fn a_traction_count_mismatch_is_refused() {
+        let surface = WettedSurface::build(&fluid_faces(), &structure_nodes()).expect("buildable");
+        let wrong = vec![Vector3::new(1.0, 0.0, 0.0); 2];
+        assert!(matches!(
+            surface.transfer_load(&fluid_faces(), &wrong),
+            Err(FsiError::CountMismatch { .. })
+        ));
+    }
+
+    #[test]
+    fn a_non_finite_traction_is_refused() {
+        let faces = fluid_faces();
+        let surface = WettedSurface::build(&faces, &structure_nodes()).expect("buildable");
+        let mut tractions = vec![Vector3::new(1.0, 0.0, 0.0); faces.len()];
+        tractions[1].y = f64::NAN;
+        assert!(matches!(
+            surface.transfer_load(&faces, &tractions),
+            Err(FsiError::NonFinite { .. })
+        ));
+    }
+}