rtx-fea + rtx-cfd: the single-step seams FSI2 stands on

rtx-fea: NonlinearDynamicAnalysis refactored onto a NonlinearDynamicStepper
- set_nodal_forces on both (the interface load of a coupling subiteration,
  replaceable between steps and between subiterations of one step);
- step(&DynamicState) is a pure function of the start-of-step state and
  the current forces - commits nothing, so a partitioned coupling re-runs
  one Newmark step to the interface fixed point (the piston semantics);
- run() marches through the same stepper: one code path, pinned from both
  ends (linear limit, CSM3, and a new manual-drive == run() assertion);
- new test: a nodal step load oscillates about the *static* nonlinear
  analysis's deflection (cross-code-path, mean within 3%, amplitude 6%),
  with re-run determinism and force-swap sensitivity asserted mid-march
  (a one-step response to a force change is ~ beta dt^2 - the first
  assertion draft demanded 10% and was corrected against the physics).

rtx-cfd: the subiteration seam and the moving no-slip closure
- EmbeddedPisoSolver::snapshot()/restore() (mask + time + init flag; the
  mask is now Clone): re-running a fluid step within a subiteration is
  bit-identical to never having diverted - proven on a moving body with
  cells flipping in the re-run window;
- polygon_interface_velocity: nearest-edge linear interpolation of
  per-vertex velocities, exact for the linear-along-edge boundary data a
  finite-element interface hands over - the no-slip closure that replaces
  FSI1's zero-velocity polygon.

Suites: rtx-fea 567, rtx-cfd 325, rtx-fsi piston+transfer - all green.

Co-Authored-By: Claude Fable 5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01Lnyrw33Lu6rUhW42E9KHwq

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

@@ -87,6 +87,14 @@ impl Default for EmbeddedParameters {
     }
 }
 
+/// A snapshot of [`EmbeddedPisoSolver`]'s per-step state, for re-running a
+/// step within a coupling subiteration. See [`EmbeddedPisoSolver::snapshot`].
+pub struct EmbeddedSolverState {
+    mask: Option<EmbeddedMask>,
+    time: f64,
+    initialized: bool,
+}
+
 /// Result of one embedded PISO step.
 #[derive(Debug, Clone)]
 pub struct EmbeddedResult {
@@ -191,6 +199,29 @@ impl EmbeddedPisoSolver {
         self.time
     }
 
+    /// Snapshot of the solver's own per-step state — the mask, the
+    /// accumulated time and the initialization flag. A coupling
+    /// subiteration re-runs one step from the same start: clone the
+    /// [`FlowField`], take this snapshot, and [`Self::restore`] both
+    /// before every re-run — otherwise the moving-body path's fresh-cell
+    /// detection compares against the *previous subiteration's* mask
+    /// instead of the committed step-start mask.
+    pub fn snapshot(&self) -> EmbeddedSolverState {
+        EmbeddedSolverState {
+            mask: self.mask.clone(),
+            time: self.time,
+            initialized: self.initialized,
+        }
+    }
+
+    /// Restore a [`Self::snapshot`]. The snapshot is cloned, so one
+    /// snapshot serves any number of re-runs.
+    pub fn restore(&mut self, state: &EmbeddedSolverState) {
+        self.mask = state.mask.clone();
+        self.time = state.time;
+        self.initialized = state.initialized;
+    }
+
     /// Reset the accumulated time.
     pub fn set_time(&mut self, t: f64) {
         self.time = t;

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

@@ -174,8 +174,7 @@ impl EmbeddedBody {
             * 0.5;
         let ccw = signed_area > 0.0;
         let sdf_vertices = vertices.clone();
-        let mut body =
+        let mut body = Self::from_sdf(move |x, y, _| polygon_signed_distance(&sdf_vertices, x, y));
-            Self::from_sdf(move |x, y, _| polygon_signed_distance(&sdf_vertices, x, y));
         let sampler_vertices = vertices;
         body.sampler = Some(Box::new(move |ds| {
             let n = sampler_vertices.len();
@@ -291,7 +290,6 @@ pub struct SurfaceForce {
     pub skipped: usize,
 }
 
-
 /// Signed distance to a closed polygon (negative inside, either winding):
 /// minimum distance over the edges, sign by the even-odd ray-crossing rule.
 /// Public so a coupling loop can build a time-dependent body from a shared,
@@ -324,6 +322,50 @@ pub fn polygon_signed_distance(vertices: &[(f64, f64)], x: f64, y: f64) -> f64 {
     if inside { -dist } else { dist }
 }
 
+/// Velocity of the point on a closed polygon nearest to `(x, y)`, where
+/// the vertices carry velocities: the nearest edge point is found exactly
+/// as in [`polygon_signed_distance`], and that edge's endpoint velocities
+/// are interpolated linearly along it. This is the no-slip closure of a
+/// deforming body whose boundary nodes move with known velocities — exact
+/// wherever the boundary velocity is linear along an edge, which is what a
+/// finite-element interface hands over. `velocities` must have one entry
+/// per vertex.
+#[must_use]
+pub fn polygon_interface_velocity(
+    vertices: &[(f64, f64)],
+    velocities: &[(f64, f64)],
+    x: f64,
+    y: f64,
+) -> (f64, f64) {
+    assert_eq!(
+        vertices.len(),
+        velocities.len(),
+        "one velocity per polygon vertex"
+    );
+    let n = vertices.len();
+    let mut best = (f64::MAX, 0usize, 0.0f64);
+    for k in 0..n {
+        let (ax, ay) = vertices[k];
+        let (bx, by) = vertices[(k + 1) % n];
+        let (ex, ey) = (bx - ax, by - ay);
+        let len2 = ex * ex + ey * ey;
+        let s = if len2 > 0.0 {
+            (((x - ax) * ex + (y - ay) * ey) / len2).clamp(0.0, 1.0)
+        } else {
+            0.0
+        };
+        let (qx, qy) = (ax + s * ex - x, ay + s * ey - y);
+        let d2 = qx * qx + qy * qy;
+        if d2 < best.0 {
+            best = (d2, k, s);
+        }
+    }
+    let (_, k, s) = best;
+    let (vax, vay) = velocities[k];
+    let (vbx, vby) = velocities[(k + 1) % n];
+    (vax + s * (vbx - vax), vay + s * (vby - vay))
+}
+
 /// What a velocity face is.
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
 pub enum FaceKind {
@@ -369,7 +411,10 @@ struct Ghost {
     flux_sign: f64,
 }
 
-/// Classification of a grid against a body at one instant.
+/// Classification of a grid against a body at one instant. `Clone` so a
+/// coupling loop can snapshot the solver's step state and re-run a step
+/// within a subiteration ([`super::EmbeddedPisoSolver::snapshot`]).
+#[derive(Clone)]
 pub struct EmbeddedMask {
     nx: usize,
     ny: usize,
@@ -1125,6 +1170,28 @@ mod tests {
         }
     }
 
+    #[test]
+    fn interface_velocity_interpolates_along_the_nearest_edge() {
+        // Unit square, CCW; each vertex carries a distinct velocity.
+        let vertices = vec![(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)];
+        let velocities = vec![(0.0, 0.0), (1.0, -1.0), (2.0, 4.0), (3.0, 9.0)];
+
+        // Near the bottom edge at s = 0.25: linear interpolation of the
+        // edge's endpoint velocities, regardless of the offset off the edge.
+        let (u, v) = polygon_interface_velocity(&vertices, &velocities, 0.25, -0.3);
+        assert!((u - 0.25).abs() < 1e-14 && (v + 0.25).abs() < 1e-14);
+        let (u, v) = polygon_interface_velocity(&vertices, &velocities, 0.25, 0.1);
+        assert!((u - 0.25).abs() < 1e-14 && (v + 0.25).abs() < 1e-14);
+
+        // Near a vertex (outside the corner): the vertex velocity.
+        let (u, v) = polygon_interface_velocity(&vertices, &velocities, 1.2, 1.3);
+        assert!((u - 2.0).abs() < 1e-14 && (v - 4.0).abs() < 1e-14);
+
+        // Midpoint of the right edge.
+        let (u, v) = polygon_interface_velocity(&vertices, &velocities, 1.4, 0.5);
+        assert!((u - 1.5).abs() < 1e-14 && (v - 1.5).abs() < 1e-14);
+    }
+
     #[test]
     fn body_touching_the_boundary_is_refused() {
         let body = EmbeddedBody::circle(0.0, 0.5, 0.2);

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

@@ -39,9 +39,10 @@ pub use ale::{
 pub use boundary_conditions::{
     BoundaryCondition, BoundaryConditions, BoundaryLocation, BoundaryType,
 };
-pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult};
+pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult, EmbeddedSolverState};
 pub use embedded_body::{
-    EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample, polygon_signed_distance,
+    EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample, polygon_interface_velocity,
+    polygon_signed_distance,
 };
 pub use flow_field::FlowField;
 pub use piso::{PisoParameters, PisoResult, PisoSolver};

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

@@ -168,6 +168,73 @@ async fn a_stationary_body_through_the_moving_path_is_bit_identical() -> CfdResu
     Ok(())
 }
 
+/// The subiteration seam: snapshot the solver + clone the field mid-run of
+/// a MOVING body, advance further (a discarded coupling candidate), then
+/// restore and advance the same steps again — the re-run must be
+/// bit-identical to a run that never diverted. This is what lets an FSI
+/// coupling re-run one fluid step under updated interface geometry.
+#[tokio::test]
+async fn snapshot_restore_rerun_is_bit_identical() -> CfdResult<()> {
+    let n = 24;
+    let dt = time_step(n);
+    let mover = || phantom_circle(|t| 0.42 + 0.30 * t, |t| 0.48 + 0.15 * t, 0.2);
+
+    // Reference: an uninterrupted run of 30 steps.
+    let mut reference = solver(n)?;
+    reference.set_moving_body(mover());
+    let mut ref_field = exact_field(n)?;
+    for _ in 0..30 {
+        reference.advance(&mut ref_field, dt).await?;
+    }
+
+    // Diverted run: 18 steps, snapshot, 12 steps of a discarded candidate,
+    // restore, the real 12 steps.
+    let mut solver_d = solver(n)?;
+    solver_d.set_moving_body(mover());
+    let mut field = exact_field(n)?;
+    for _ in 0..18 {
+        solver_d.advance(&mut field, dt).await?;
+    }
+    let saved_state = solver_d.snapshot();
+    let saved_field = field.clone();
+    for _ in 0..12 {
+        solver_d.advance(&mut field, dt).await?; // discarded candidate
+    }
+    solver_d.restore(&saved_state);
+    field = saved_field;
+    let mut rerun_fresh = 0usize;
+    for _ in 0..12 {
+        rerun_fresh += solver_d.advance(&mut field, dt).await?.fresh_cells;
+    }
+    // Cells must actually flip in the re-run window, or the restore of the
+    // mask was never exercised against a mask that changes.
+    assert!(
+        rerun_fresh > 0,
+        "no cells flipped after the restore — the test is vacuous"
+    );
+
+    assert_eq!(
+        solver_d.time().to_bits(),
+        reference.time().to_bits(),
+        "restored time diverges"
+    );
+    let mut max_diff: f64 = 0.0;
+    for (x, y) in field.u.iter().zip(ref_field.u.iter()) {
+        max_diff = max_diff.max((x - y).abs());
+    }
+    for (x, y) in field.v.iter().zip(ref_field.v.iter()) {
+        max_diff = max_diff.max((x - y).abs());
+    }
+    for (x, y) in field.p.iter().zip(ref_field.p.iter()) {
+        max_diff = max_diff.max((x - y).abs());
+    }
+    assert!(
+        max_diff == 0.0,
+        "restored re-run differs from the uninterrupted run by {max_diff:.3e}"
+    );
+    Ok(())
+}
+
 /// Claim 2: the translating phantom circle. Static steady-state baselines
 /// at n = 32 (upwind, from `tests/embedded_mms.rs`): L2 u 8.489e-3,
 /// L2 p 2.22e-2.

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

@@ -22,12 +22,28 @@
 //! a total-Lagrangian setting); no damping (Rayleigh damping can be added
 //! when something needs it — the Turek–Hron CSM3 benchmark is undamped).
 //!
+//! # Two ways to drive it
+//!
+//! [`NonlinearDynamicAnalysis::run`] marches `num_steps` steps from rest —
+//! the benchmark shape (CSM3: gravity switched on at rest).
+//!
+//! [`NonlinearDynamicAnalysis::stepper`] hands out the same machinery one
+//! step at a time, for a partitioned coupling loop: the caller owns the
+//! state ([`DynamicState`]), sets the interface load with
+//! [`NonlinearDynamicStepper::set_nodal_forces`], and calls
+//! [`NonlinearDynamicStepper::step`] — which reads the start-of-step state
+//! and *does not commit anything*, so a subiteration can re-run the same
+//! step from the same state under an updated load as many times as the
+//! interface fixed point takes (the semantics the coupled piston benchmark
+//! established). `run` is implemented on the stepper, so the benchmark
+//! tests pin both.
+//!
 //! Limits, stated up front: Dirichlet conditions must be homogeneous
-//! (`u = 0` — a clamped edge; prescribed motion belongs to the FSI rung
+//! (`u = 0` — a clamped edge); the body force is constant in time, applied
-//! and enters through `set_prescribed_history` when that lands); the body
+//! fully from `t = 0` (CSM3's definition: gravity switched on at rest, the
-//! force is constant in time, applied fully from `t = 0` (CSM3's
+//! structure oscillates about its static deflection). Nodal forces may
-//! definition: gravity switched on at rest, the structure oscillates about
+//! change between steps (and between subiterations of one step) through
-//! its static deflection).
+//! the stepper.
 
 use super::{AnalysisConfig, ConvergenceCriteria};
 use crate::assembly::SparseMatrix;
@@ -61,6 +77,30 @@ pub struct NonlinearDynamicResults {
     pub max_iterations_per_step: usize,
 }
 
+/// The full kinematic state at one instant: displacement, velocity and
+/// acceleration as full-length vectors under the analysis's DOF numbering
+/// (constrained entries zero). The caller owns it; a coupling loop clones
+/// the committed state and re-steps from it freely.
+#[derive(Debug, Clone)]
+pub struct DynamicState {
+    /// Displacement.
+    pub displacement: DVector<f64>,
+    /// Velocity.
+    pub velocity: DVector<f64>,
+    /// Acceleration.
+    pub acceleration: DVector<f64>,
+}
+
+/// Per-element setup computed once: coordinates, DOFs, the DOF-expanded
+/// consistent mass (configuration-independent), and the material.
+struct ElementCache {
+    coords: Vec<Vector3<f64>>,
+    dofs: Vec<usize>,
+    element_type: crate::mesh::ElementType,
+    mass: DMatrix<f64>,
+    material_id: crate::mesh::MaterialId,
+}
+
 /// Nonlinear Newmark transient analysis. See the module docs.
 pub struct NonlinearDynamicAnalysis {
     mesh: Mesh,
@@ -76,6 +116,7 @@ pub struct NonlinearDynamicAnalysis {
     total_lagrangian: bool,
     #[allow(clippy::type_complexity)]
     body_force: Option<Box<dyn Fn(Vector3<f64>) -> Vector3<f64> + Send + Sync>>,
+    nodal_forces: Vec<(NodeId, Vector3<f64>)>,
     tracked_nodes: Vec<NodeId>,
 }
 
@@ -102,6 +143,7 @@ impl NonlinearDynamicAnalysis {
             beta: 0.25,
             total_lagrangian: false,
             body_force: None,
+            nodal_forces: Vec::new(),
             tracked_nodes: Vec::new(),
         }
     }
@@ -137,23 +179,113 @@ impl NonlinearDynamicAnalysis {
         self.body_force = Some(Box::new(f));
     }
 
-    /// Record this node's displacement at every step.
+    /// Concentrated nodal forces, added to the external force (as on
+    /// [`super::NonlinearStaticAnalysis`]). For a load that changes in
+    /// time, use [`Self::stepper`] and set the forces before each step.
+    pub fn set_nodal_forces(&mut self, forces: Vec<(NodeId, Vector3<f64>)>) {
+        self.nodal_forces = forces;
+    }
+
+    /// Record this node's displacement at every step of [`Self::run`].
     pub fn track_node(&mut self, node: NodeId) {
         self.tracked_nodes.push(node);
     }
 
-    /// March `num_steps` steps of `dt` from rest.
+    /// Build the single-step driver: DOF numbering, element caches, the
-    #[allow(clippy::too_many_lines)]
+    /// consistent mass, and the constant external force, assembled once.
+    pub fn stepper(&self) -> FeaResult<NonlinearDynamicStepper<'_>> {
+        NonlinearDynamicStepper::build(self)
+    }
+
+    /// March `num_steps` steps of `dt` from rest, via the same stepper a
+    /// coupling loop would drive.
     pub fn run(&mut self) -> FeaResult<NonlinearDynamicResults> {
-        let dim = self.mesh.spatial_dimension;
+        let mut stepper = self.stepper()?;
-        let mut dof_numbering =
+        let mut state = stepper.rest_state()?;
+
+        let mut times = Vec::with_capacity(self.num_steps);
+        let mut tracked: Vec<Vec<DVector<f64>>> =
+            vec![Vec::with_capacity(self.num_steps); self.tracked_nodes.len()];
+        let mut total_iterations = 0usize;
+        let mut max_iterations_per_step = 0usize;
+
+        for step in 1..=self.num_steps {
+            let (new_state, iterations) = stepper.step(&state)?;
+            state = new_state;
+            total_iterations += iterations;
+            max_iterations_per_step = max_iterations_per_step.max(iterations);
+
+            times.push(step as f64 * self.dt);
+            for (slot, node) in self.tracked_nodes.iter().enumerate() {
+                let dofs = stepper.node_dofs(*node);
+                let mut value = DVector::zeros(dofs.len());
+                for (c, &dof) in dofs.iter().enumerate() {
+                    value[c] = state.displacement[dof];
+                }
+                tracked[slot].push(value);
+            }
+        }
+
+        Ok(NonlinearDynamicResults {
+            times,
+            tracked,
+            displacement: state.displacement,
+            velocity: state.velocity,
+            acceleration: state.acceleration,
+            total_iterations,
+            max_iterations_per_step,
+        })
+    }
+
+    /// The DOF indices of a node under the analysis's own numbering, for
+    /// reading the returned full-length vectors.
+    pub fn node_dofs(&self, node: NodeId) -> FeaResult<Vec<usize>> {
+        let numbering =
             AdvancedDofNumbering::displacement_only(&self.mesh, DofMappingStrategy::Sequential)?;
+        Ok(numbering.get_node_dofs(node))
+    }
+}
+
+/// The single-step Newmark–Newton driver behind
+/// [`NonlinearDynamicAnalysis`]. Holds everything assembled once (DOF
+/// numbering, element caches, consistent mass, the constant body-force
+/// vector); the mutable pieces are the nodal forces and the linear solver.
+///
+/// [`Self::step`] is a pure function of the start-of-step [`DynamicState`]
+/// and the current forces: nothing is committed, so a partitioned coupling
+/// can re-run one step under updated interface loads until the interface
+/// converges, then keep the accepted state.
+pub struct NonlinearDynamicStepper<'a> {
+    analysis: &'a NonlinearDynamicAnalysis,
+    dof_numbering: AdvancedDofNumbering,
+    free_dofs: Vec<usize>,
+    free_index: Vec<Option<usize>>,
+    total_dofs: usize,
+    caches: Vec<ElementCache>,
+    /// Free-free consistent mass, for consistent initial accelerations.
+    mass_free: SparseMatrix,
+    /// The body-force part of the external force (constant).
+    external_body: DVector<f64>,
+    /// Body force plus the current nodal forces.
+    external: DVector<f64>,
+    solver: LuDirect,
+    solver_options: SolverOptions,
+}
+
+impl<'a> NonlinearDynamicStepper<'a> {
+    #[allow(clippy::too_many_lines)]
+    fn build(analysis: &'a NonlinearDynamicAnalysis) -> FeaResult<Self> {
+        let dim = analysis.mesh.spatial_dimension;
+        let mut dof_numbering = AdvancedDofNumbering::displacement_only(
+            &analysis.mesh,
+            DofMappingStrategy::Sequential,
+        )?;
 
         // Homogeneous Dirichlet only (see the module docs).
-        for condition in self.boundary_conditions.conditions() {
+        for condition in analysis.boundary_conditions.conditions() {
             if let BoundaryCondition::Dirichlet(dirichlet) = condition {
                 for &node in &dirichlet.nodes {
-                    let position = self
+                    let position = analysis
                         .mesh
                         .get_node(node)
                         .ok_or_else(|| {
@@ -191,28 +323,19 @@ impl NonlinearDynamicAnalysis {
             free_index[dof] = Some(k);
         }
 
-        // Per-element caches: coordinates, DOFs, the DOF-expanded consistent
+        let mut caches = Vec::with_capacity(analysis.mesh.elements.len());
-        // mass (configuration-independent), and the constitutive choice.
+        for element in analysis.mesh.elements.values() {
-        struct ElementCache {
-            coords: Vec<Vector3<f64>>,
-            dofs: Vec<usize>,
-            element_type: crate::mesh::ElementType,
-            mass: DMatrix<f64>,
-            material_id: crate::mesh::MaterialId,
-        }
-        let mut caches = Vec::with_capacity(self.mesh.elements.len());
-        for element in self.mesh.elements.values() {
             let coords: Vec<Vector3<f64>> = element
                 .nodes
                 .iter()
-                .map(|id| self.mesh.get_node(*id).unwrap().position())
+                .map(|id| analysis.mesh.get_node(*id).unwrap().position())
                 .collect();
             let dofs: Vec<usize> = element
                 .nodes
                 .iter()
                 .flat_map(|node| dof_numbering.get_node_dofs(*node))
                 .collect();
-            let material = self
+            let material = analysis
                 .materials
                 .get_material(element.material_id)
                 .ok_or_else(|| {
@@ -244,86 +367,28 @@ impl NonlinearDynamicAnalysis {
             });
         }
 
-        // Internal force and (optionally) tangent at a full displacement
+        // Free-free consistent mass (for initial accelerations).
-        // vector, reduced to the free DOFs.
+        let mut mass_free = SparseMatrix::new(num_free, num_free);
-        let assemble = |solution: &DVector<f64>,
+        for cache in &caches {
-                        with_tangent: bool|
+            for (local_row, &dof_row) in cache.dofs.iter().enumerate() {
-         -> FeaResult<(DVector<f64>, SparseMatrix)> {
+                let Some(free_row) = free_index[dof_row] else {
-            let mut internal = DVector::zeros(num_free);
+                    continue;
-            let mut tangent = SparseMatrix::new(num_free, num_free);
-            let inv_beta_dt2 = 1.0 / (self.beta * self.dt * self.dt);
-            for cache in &caches {
-                let material = self.materials.get_material(cache.material_id).unwrap();
-                let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone());
-                let mut element_displacement = DVector::zeros(cache.dofs.len());
-                for (local, &dof) in cache.dofs.iter().enumerate() {
-                    element_displacement[local] = solution[dof];
-                }
-                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,
-                        &cache.coords,
-                        &element_displacement,
-                        constitutive.as_ref(),
-                        None,
-                    )?
-                } else {
-                    let constitutive = reduced_constitutive(material, dim)?;
-                    ElementMatrixComputer::compute_internal_force_and_tangent(
-                        &fe,
-                        &cache.coords,
-                        &element_displacement,
-                        constitutive.as_ref(),
-                        None,
-                    )?
                 };
-                for (local_row, &dof_row) in cache.dofs.iter().enumerate() {
+                for (local_col, &dof_col) in cache.dofs.iter().enumerate() {
-                    let Some(free_row) = free_index[dof_row] else {
+                    if let Some(free_col) = free_index[dof_col] {
-                        continue;
+                        let value = cache.mass[(local_row, local_col)];
-                    };
+                        if value != 0.0 {
-                    internal[free_row] += f_int[local_row];
+                            mass_free.add_entry(free_row, free_col, value)?;
-                    if with_tangent {
-                        for (local_col, &dof_col) in cache.dofs.iter().enumerate() {
-                            if let Some(free_col) = free_index[dof_col] {
-                                let value = k_t[(local_row, local_col)]
-                                    + inv_beta_dt2 * cache.mass[(local_row, local_col)];
-                                if value != 0.0 {
-                                    tangent.add_entry(free_row, free_col, value)?;
-                                }
-                            }
                         }
                     }
                 }
             }
-            if with_tangent {
+        }
-                tangent.finalize()?;
+        mass_free.finalize()?;
-            }
-            Ok((internal, tangent))
-        };
 
-        // M times a full-length vector, reduced to the free DOFs.
+        // Constant consistent external force from the body-force field.
-        let mass_times = |a_full: &DVector<f64>| -> DVector<f64> {
+        let mut external_body: DVector<f64> = DVector::zeros(num_free);
-            let mut out = DVector::zeros(num_free);
+        if let Some(force) = &analysis.body_force {
-            for cache in &caches {
-                let mut a_e = DVector::zeros(cache.dofs.len());
-                for (local, &dof) in cache.dofs.iter().enumerate() {
-                    a_e[local] = a_full[dof];
-                }
-                let m_a = &cache.mass * a_e;
-                for (local, &dof) in cache.dofs.iter().enumerate() {
-                    if let Some(free) = free_index[dof] {
-                        out[free] += m_a[local];
-                    }
-                }
-            }
-            out
-        };
-
-        // Constant consistent external force.
-        let mut external: DVector<f64> = DVector::zeros(num_free);
-        if let Some(force) = &self.body_force {
             for cache in &caches {
                 let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone());
                 let f_e = ElementMatrixComputer::compute_body_force_vector(
@@ -334,141 +399,220 @@ impl NonlinearDynamicAnalysis {
                 )?;
                 for (local, &dof) in cache.dofs.iter().enumerate() {
                     if let Some(free) = free_index[dof] {
-                        external[free] += f_e[local];
+                        external_body[free] += f_e[local];
                     }
                 }
             }
         }
-        let force_scale = external.norm().max(1.0);
 
-        // State (full-length vectors; constrained entries stay zero).
+        let mut stepper = Self {
-        let mut u = DVector::zeros(total_dofs);
+            analysis,
-        let mut v = DVector::zeros(total_dofs);
+            dof_numbering,
-        let mut a = DVector::zeros(total_dofs);
+            free_dofs,
+            free_index,
+            total_dofs,
+            caches,
+            mass_free,
+            external_body: external_body.clone(),
+            external: external_body,
+            solver: LuDirect::new(),
+            solver_options: SolverOptions::default(),
+        };
+        stepper.set_nodal_forces(&analysis.nodal_forces);
+        Ok(stepper)
+    }
 
-        let mut solver = LuDirect::new();
+    /// Replace the concentrated nodal forces (the interface load of a
-        let solver_options = SolverOptions::default();
+    /// coupling subiteration). The body-force part is unaffected.
+    pub fn set_nodal_forces(&mut self, forces: &[(NodeId, Vector3<f64>)]) {
+        self.external.copy_from(&self.external_body);
+        for (node, force) in forces {
+            let dofs = self.dof_numbering.get_node_dofs(*node);
+            for (component, &dof) in dofs.iter().enumerate() {
+                if let Some(free) = self.free_index[dof] {
+                    self.external[free] += force[component];
+                }
+            }
+        }
+    }
 
-        // Initial acceleration from rest: M a0 = F_ext - f_int(0). The mass
+    /// The state at rest under the *current* external force: `u = v = 0`,
-        // block is SPD; reuse the Newton machinery by solving with the
+    /// the acceleration consistent with `M a0 = F_ext - f_int(0)`.
-        // effective matrix at beta dt^2 = 1 scaling of the mass alone —
+    pub fn rest_state(&mut self) -> FeaResult<DynamicState> {
-        // assemble M (free-free) once through the tangent path with a zero
+        let u = DVector::zeros(self.total_dofs);
-        // stiffness contribution is not available, so build it directly.
+        let (f_int0, _) = self.assemble(&u, false)?;
-        {
+        let residual0 = &self.external - &f_int0;
-            let mut mass_free = SparseMatrix::new(num_free, num_free);
+        let (a0_free, _) = self
-            for cache in &caches {
+            .solver
-                for (local_row, &dof_row) in cache.dofs.iter().enumerate() {
+            .solve(&self.mass_free, &residual0, &self.solver_options)?;
-                    let Some(free_row) = free_index[dof_row] else {
+        let mut a = DVector::zeros(self.total_dofs);
-                        continue;
+        for (k, &dof) in self.free_dofs.iter().enumerate() {
-                    };
+            a[dof] = a0_free[k];
+        }
+        Ok(DynamicState {
+            displacement: DVector::zeros(self.total_dofs),
+            velocity: DVector::zeros(self.total_dofs),
+            acceleration: a,
+        })
+    }
+
+    /// One Newmark step of the analysis's `dt` from `state` under the
+    /// current forces. Returns the end-of-step state and the Newton
+    /// iteration count; commits nothing — calling again with the same
+    /// state and forces returns the identical result.
+    pub fn step(&mut self, state: &DynamicState) -> FeaResult<(DynamicState, usize)> {
+        let dt = self.analysis.dt;
+        let gamma = self.analysis.gamma;
+        let beta = self.analysis.beta;
+        let criteria = &self.analysis.criteria;
+        let force_scale = self.external.norm().max(1.0);
+
+        let mut u_pred = DVector::zeros(self.total_dofs);
+        for &dof in &self.free_dofs {
+            u_pred[dof] = state.displacement[dof]
+                + dt * state.velocity[dof]
+                + dt * dt * (0.5 - beta) * state.acceleration[dof];
+        }
+        let inv_beta_dt2 = 1.0 / (beta * dt * dt);
+
+        // Newton on the end-of-step displacement, starting from the
+        // predictor (a_new = 0 there).
+        let mut u_iter = u_pred.clone();
+        let mut step_converged = false;
+        let mut iterations = 0usize;
+        for _ in 0..criteria.max_iterations {
+            let mut a_new = DVector::zeros(self.total_dofs);
+            for &dof in &self.free_dofs {
+                a_new[dof] = inv_beta_dt2 * (u_iter[dof] - u_pred[dof]);
+            }
+            let (f_int, tangent) = self.assemble(&u_iter, true)?;
+            let residual = &self.external - &f_int - self.mass_times(&a_new);
+            if residual.norm() < criteria.force_tolerance * force_scale {
+                step_converged = true;
+                break;
+            }
+            iterations += 1;
+            let (delta, _) = self
+                .solver
+                .solve(&tangent, &residual, &self.solver_options)?;
+            for (k, &dof) in self.free_dofs.iter().enumerate() {
+                u_iter[dof] += delta[k];
+            }
+            if delta.norm() < criteria.displacement_tolerance * u_iter.norm().max(1.0) {
+                step_converged = true;
+                break;
+            }
+        }
+        if !step_converged {
+            return Err(AnalysisError::ConvergenceFailed { iterations }.into());
+        }
+
+        let mut a_new = DVector::zeros(self.total_dofs);
+        let mut v_new = DVector::zeros(self.total_dofs);
+        for &dof in &self.free_dofs {
+            a_new[dof] = inv_beta_dt2 * (u_iter[dof] - u_pred[dof]);
+            v_new[dof] = state.velocity[dof]
+                + dt * ((1.0 - gamma) * state.acceleration[dof] + gamma * a_new[dof]);
+        }
+        Ok((
+            DynamicState {
+                displacement: u_iter,
+                velocity: v_new,
+                acceleration: a_new,
+            },
+            iterations,
+        ))
+    }
+
+    /// The DOF indices of a node, for reading [`DynamicState`] vectors.
+    pub fn node_dofs(&self, node: NodeId) -> Vec<usize> {
+        self.dof_numbering.get_node_dofs(node)
+    }
+
+    /// Internal force and (optionally) tangent at a full displacement
+    /// vector, reduced to the free DOFs. The tangent includes the Newmark
+    /// mass term `M / (β Δt²)`.
+    fn assemble(
+        &self,
+        solution: &DVector<f64>,
+        with_tangent: bool,
+    ) -> FeaResult<(DVector<f64>, SparseMatrix)> {
+        let dim = self.analysis.mesh.spatial_dimension;
+        let num_free = self.free_dofs.len();
+        let mut internal = DVector::zeros(num_free);
+        let mut tangent = SparseMatrix::new(num_free, num_free);
+        let inv_beta_dt2 = 1.0 / (self.analysis.beta * self.analysis.dt * self.analysis.dt);
+        for cache in &self.caches {
+            let material = self
+                .analysis
+                .materials
+                .get_material(cache.material_id)
+                .unwrap();
+            let fe = StandardFiniteElement::new(cache.element_type, cache.coords.clone());
+            let mut element_displacement = DVector::zeros(cache.dofs.len());
+            for (local, &dof) in cache.dofs.iter().enumerate() {
+                element_displacement[local] = solution[dof];
+            }
+            let (f_int, k_t) = if self.analysis.total_lagrangian {
+                let (lambda, mu) = material.properties().lame_parameters();
+                let constitutive = saint_venant_kirchhoff(lambda, mu, dim);
+                total_lagrangian::internal_force_and_tangent(
+                    &fe,
+                    &cache.coords,
+                    &element_displacement,
+                    constitutive.as_ref(),
+                    None,
+                )?
+            } else {
+                let constitutive = reduced_constitutive(material, dim)?;
+                ElementMatrixComputer::compute_internal_force_and_tangent(
+                    &fe,
+                    &cache.coords,
+                    &element_displacement,
+                    constitutive.as_ref(),
+                    None,
+                )?
+            };
+            for (local_row, &dof_row) in cache.dofs.iter().enumerate() {
+                let Some(free_row) = self.free_index[dof_row] else {
+                    continue;
+                };
+                internal[free_row] += f_int[local_row];
+                if with_tangent {
                     for (local_col, &dof_col) in cache.dofs.iter().enumerate() {
-                        if let Some(free_col) = free_index[dof_col] {
+                        if let Some(free_col) = self.free_index[dof_col] {
-                            let value = cache.mass[(local_row, local_col)];
+                            let value = k_t[(local_row, local_col)]
+                                + inv_beta_dt2 * cache.mass[(local_row, local_col)];
                             if value != 0.0 {
-                                mass_free.add_entry(free_row, free_col, value)?;
+                                tangent.add_entry(free_row, free_col, value)?;
                             }
                         }
                     }
                 }
             }
-            mass_free.finalize()?;
-            let (f_int0, _) = assemble(&u, false)?;
-            let residual0 = &external - &f_int0;
-            let (a0_free, _) = solver.solve(&mass_free, &residual0, &solver_options)?;
-            for (k, &dof) in free_dofs.iter().enumerate() {
-                a[dof] = a0_free[k];
-            }
         }
-
+        if with_tangent {
-        let mut times = Vec::with_capacity(self.num_steps);
+            tangent.finalize()?;
-        let mut tracked: Vec<Vec<DVector<f64>>> =
-            vec![Vec::with_capacity(self.num_steps); self.tracked_nodes.len()];
-        let mut total_iterations = 0usize;
-        let mut max_iterations_per_step = 0usize;
-
-        for step in 1..=self.num_steps {
-            // Predictor and the Newmark kinematics.
-            let mut u_pred = DVector::zeros(total_dofs);
-            for &dof in &free_dofs {
-                u_pred[dof] =
-                    u[dof] + self.dt * v[dof] + self.dt * self.dt * (0.5 - self.beta) * a[dof];
-            }
-            let inv_beta_dt2 = 1.0 / (self.beta * self.dt * self.dt);
-
-            // Newton on the end-of-step displacement, starting from the
-            // predictor (a_new = 0 there).
-            let mut u_iter = u_pred.clone();
-            let mut step_converged = false;
-            let mut iterations = 0usize;
-            for _ in 0..self.criteria.max_iterations {
-                let mut a_new = DVector::zeros(total_dofs);
-                for &dof in &free_dofs {
-                    a_new[dof] = inv_beta_dt2 * (u_iter[dof] - u_pred[dof]);
-                }
-                let (f_int, tangent) = assemble(&u_iter, true)?;
-                let residual = &external - &f_int - mass_times(&a_new);
-                if residual.norm() < self.criteria.force_tolerance * force_scale {
-                    step_converged = true;
-                    break;
-                }
-                iterations += 1;
-                let (delta, _) = solver.solve(&tangent, &residual, &solver_options)?;
-                for (k, &dof) in free_dofs.iter().enumerate() {
-                    u_iter[dof] += delta[k];
-                }
-                if delta.norm() < self.criteria.displacement_tolerance * u_iter.norm().max(1.0) {
-                    step_converged = true;
-                    break;
-                }
-            }
-            if !step_converged {
-                return Err(AnalysisError::ConvergenceFailed {
-                    iterations: total_iterations + iterations,
-                }
-                .into());
-            }
-            total_iterations += iterations;
-            max_iterations_per_step = max_iterations_per_step.max(iterations);
-
-            // Accept the step.
-            let mut a_new = DVector::zeros(total_dofs);
-            let mut v_new = DVector::zeros(total_dofs);
-            for &dof in &free_dofs {
-                a_new[dof] = inv_beta_dt2 * (u_iter[dof] - u_pred[dof]);
-                v_new[dof] =
-                    v[dof] + self.dt * ((1.0 - self.gamma) * a[dof] + self.gamma * a_new[dof]);
-            }
-            u = u_iter;
-            v = v_new;
-            a = a_new;
-
-            times.push(step as f64 * self.dt);
-            for (slot, node) in self.tracked_nodes.iter().enumerate() {
-                let dofs = dof_numbering.get_node_dofs(*node);
-                let mut value = DVector::zeros(dofs.len());
-                for (c, &dof) in dofs.iter().enumerate() {
-                    value[c] = u[dof];
-                }
-                tracked[slot].push(value);
-            }
         }
-
+        Ok((internal, tangent))
-        Ok(NonlinearDynamicResults {
-            times,
-            tracked,
-            displacement: u,
-            velocity: v,
-            acceleration: a,
-            total_iterations,
-            max_iterations_per_step,
-        })
     }
 
-    /// The DOF indices of a node under the analysis's own numbering, for
+    /// M times a full-length vector, reduced to the free DOFs.
-    /// reading the returned full-length vectors.
+    fn mass_times(&self, a_full: &DVector<f64>) -> DVector<f64> {
-    pub fn node_dofs(&self, node: NodeId) -> FeaResult<Vec<usize>> {
+        let num_free = self.free_dofs.len();
-        let numbering =
+        let mut out = DVector::zeros(num_free);
-            AdvancedDofNumbering::displacement_only(&self.mesh, DofMappingStrategy::Sequential)?;
+        for cache in &self.caches {
-        Ok(numbering.get_node_dofs(node))
+            let mut a_e = DVector::zeros(cache.dofs.len());
+            for (local, &dof) in cache.dofs.iter().enumerate() {
+                a_e[local] = a_full[dof];
+            }
+            let m_a = &cache.mass * a_e;
+            for (local, &dof) in cache.dofs.iter().enumerate() {
+                if let Some(free) = self.free_index[dof] {
+                    out[free] += m_a[local];
+                }
+            }
+        }
+        out
     }
 }

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

@@ -15,7 +15,10 @@
 //!    35×2 Quad8 the static CSM tests bounded at ~1.5%.
 
 use nalgebra::{DMatrix, DVector, Vector3};
-use rtx_fea::analysis::{AnalysisConfig, NewmarkStepper, NonlinearDynamicAnalysis};
+use rtx_fea::analysis::{
+    Analysis, AnalysisConfig, NewmarkStepper, NonlinearConfig, NonlinearDynamicAnalysis,
+    NonlinearStaticAnalysis,
+};
 use rtx_fea::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
 use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
 use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
@@ -330,3 +333,117 @@ fn turek_hron_csm3_oscillation() {
         results.max_iterations_per_step
     );
 }
+
+/// 3. The stepper under a nodal step load — the FSI2 seam. Three claims:
+///
+/// a. Driving the stepper by hand (set the nodal force, step, commit) is
+///    bit-identical to `run()` with the same force set on the analysis —
+///    one code path, verified from both ends.
+/// b. `step` commits nothing: repeating a step from the same state under
+///    the same force is bit-identical; changing the force between the
+///    repeats changes the answer (the subiteration a coupling loop needs).
+/// c. The undamped step response oscillates about the static deflection:
+///    `mid_amp` of the tip trajectory must match the *static* nonlinear
+///    analysis under the identical nodal force — a different code path —
+///    in both mean and amplitude (`u(t) ≈ u_s (1 − cos ωt)` while the
+///    first mode dominates a tip-loaded cantilever).
+#[test]
+fn stepper_nodal_step_load_oscillates_about_the_static_deflection() {
+    let mesh = quad8_rect_mesh(0.25, 0.6, 0.19, 0.21, 10, 2);
+    let dt = 0.005;
+    let steps = 400; // 2 s: two periods of the ~1 Hz first mode
+    let a_node = point_a(&mesh, 0.6, 0.2);
+    let tip_force = Vector3::new(0.0, -0.1, 0.0);
+
+    // run() with the force set on the analysis.
+    let mut analysis = NonlinearDynamicAnalysis::new(
+        mesh.clone(),
+        materials(),
+        clamp_left(&mesh, 0.25),
+        dt,
+        steps,
+        AnalysisConfig::default(),
+    )
+    .with_total_lagrangian();
+    analysis.set_nodal_forces(vec![(a_node, tip_force)]);
+    analysis.track_node(a_node);
+    let results = analysis.run().unwrap();
+    let uy_run: Vec<f64> = results.tracked[0].iter().map(|u| u[1]).collect();
+
+    // The same march, driven by hand through the stepper.
+    let mut stepper = analysis.stepper().unwrap();
+    stepper.set_nodal_forces(&[(a_node, tip_force)]);
+    let mut state = stepper.rest_state().unwrap();
+    let a_dofs = stepper.node_dofs(a_node);
+    let mut uy_manual = Vec::with_capacity(steps);
+    for step in 0..steps {
+        if step == 7 {
+            // b. Re-running the same step is bit-identical; a different
+            // force from the same state gives a different answer and
+            // leaves no trace once the force is restored.
+            let (first, _) = stepper.step(&state).unwrap();
+            let (again, _) = stepper.step(&state).unwrap();
+            assert_eq!(
+                first.displacement, again.displacement,
+                "re-running a step from the same state changed the answer"
+            );
+            stepper.set_nodal_forces(&[(a_node, 2.0 * tip_force)]);
+            let (other, _) = stepper.step(&state).unwrap();
+            // One step's response to an extra force is ~ ΔF β Δt² / m_modal
+            // (≈ 1% of the accumulated displacement here), downward.
+            let moved = other.displacement[a_dofs[1]] - first.displacement[a_dofs[1]];
+            assert!(
+                moved < -1e-3 * first.displacement[a_dofs[1]].abs(),
+                "doubling the interface force did not move the step down: \
+                 delta {moved:.3e} vs u {:.3e}",
+                first.displacement[a_dofs[1]]
+            );
+            stepper.set_nodal_forces(&[(a_node, tip_force)]);
+        }
+        let (new_state, _) = stepper.step(&state).unwrap();
+        state = new_state;
+        uy_manual.push(state.displacement[a_dofs[1]]);
+    }
+    // a. One code path, verified from both ends.
+    assert_eq!(
+        uy_run, uy_manual,
+        "manual stepper drive deviates from run()"
+    );
+
+    // c. Static deflection under the identical nodal force, from the
+    // nonlinear *static* analysis.
+    let mut static_analysis = NonlinearStaticAnalysis::new(
+        mesh.clone(),
+        materials(),
+        clamp_left(&mesh, 0.25),
+        NonlinearConfig::default(),
+        AnalysisConfig::default(),
+    )
+    .with_total_lagrangian();
+    static_analysis.set_nodal_forces(vec![(a_node, tip_force)]);
+    let static_results = static_analysis.run().unwrap();
+    assert!(static_results.convergence.converged);
+    let numbering =
+        AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
+    let uy_static = static_results.displacements[numbering.get_node_dofs(a_node)[1]];
+
+    let (uy_mean, uy_amp) = mid_amp(&uy_manual);
+    println!(
+        "  step load: static uy = {uy_static:.4e}, dynamic mid ± amp = \
+         {uy_mean:.4e} ± {uy_amp:.4e}"
+    );
+    assert!(
+        uy_static < -1e-4,
+        "static deflection suspiciously small: {uy_static:.3e}"
+    );
+    let rel = |a: f64, b: f64| ((a - b) / b).abs();
+    assert!(
+        rel(uy_mean, uy_static) < 0.03,
+        "oscillation midpoint {uy_mean:.4e} vs static deflection {uy_static:.4e}"
+    );
+    assert!(
+        rel(uy_amp, -uy_static) < 0.06,
+        "oscillation amplitude {uy_amp:.4e} vs |static| {:.4e}",
+        -uy_static
+    );
+}