R8-b: 3-D Hex20 flag structure (rtx-fea analysis::flag3d) + plane-strain self-consistency gates

New code only: Flag3d builds the Turek-Hron flag as a structured Hex20
plate (x-major serendipity lattice), the root clamp with free or
plane-strain (u_z = 0) lateral faces, the TL-SVK Newmark analysis with the
2-D harness's settings (beta = (gamma + 1/2)^2/4, 60 Newton), the wetted
Quad8 faces with outward orientation, and consistent nodal forces of a
traction field on the current faces. The dynamic stepper and the TL path
were already dimension-generic; nothing existing changes.

Gates (tests/flag3d_structure.rs): mass = rho V, face-force totals;
plane-strain 3-D reproduces the 2-D 35x2 Quad8 CSM1 to rounding and the
first 60 CSM3 steps to 9e-16 m. Instruments (#[ignore]): CSM1 table,
CSM3 march, modal K/M dump.

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

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

@@ -0,0 +1,475 @@
+//! The Turek–Hron flag as a 3-D solid: a Hex20 plate `[x0, x1] × [y0, y1]
+//! × [z0, z1]` (length × thickness × span), root face `x = x0` clamped,
+//! for R8's coupled 3-D FSI (omni-cortex roadmap, item R8-b).
+//!
+//! Nothing here is a new solver: [`NonlinearDynamicAnalysis`] and the
+//! total-Lagrangian St. Venant–Kirchhoff path are dimension-generic, so
+//! the 3-D flag steps through exactly the machinery the 2-D harness uses
+//! (`set_nodal_forces`, `step(state) → state`, Newmark γ/β, the
+//! line-search/subdivision rescue). This module adds what the 2-D harness
+//! builds by hand:
+//!
+//! * the structured Hex20 mesh on the `(2nx+1) × (2ny+1) × (2nz+1)`
+//!   serendipity lattice (x-major node order, so the sequential DOF
+//!   numbering keeps the banded LU's bandwidth at one x-slab);
+//! * the root clamp, with the lateral faces either **free** (the real
+//!   plate, the physical 3-D problem) or **plane strain** (`u_z = 0` at
+//!   every node: a z-independent field is exactly representable by Hex20
+//!   and the 3-D energy then equals span × the 2-D plane-strain energy, so
+//!   this reproduces the 2-D Quad8 model to rounding — the self-consistency
+//!   gate, "run the 3-D problem as the 2-D problem first");
+//! * the wetted surface (bottom, top, tip, and the two lateral faces) as
+//!   Quad8 faces with outward orientation, and the consistent nodal forces
+//!   of a traction field integrated over the *current* (deformed) faces —
+//!   the interface load a partitioned coupling feeds to
+//!   [`NonlinearDynamicStepper::set_nodal_forces`](super::NonlinearDynamicStepper::set_nodal_forces).
+//!
+//! New code only: no existing solver path changes.
+
+use super::{AnalysisConfig, ConvergenceCriteria, NonlinearDynamicAnalysis};
+use crate::assembly::dof_mapping::DofComponent;
+use crate::boundary::dirichlet::{DirichletBC, DirichletType};
+use crate::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
+use crate::error::FeaResult;
+use crate::materials::{LinearElastic, MaterialDatabase};
+use crate::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+use nalgebra::{DVector, Vector3};
+use std::collections::BTreeMap;
+
+/// Geometry and resolution of the plate. `nx`, `ny`, `nz` are Hex20
+/// element counts along length, thickness and span.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Flag3dSpec {
+    pub x0: f64,
+    pub x1: f64,
+    pub y0: f64,
+    pub y1: f64,
+    pub z0: f64,
+    pub z1: f64,
+    pub nx: usize,
+    pub ny: usize,
+    pub nz: usize,
+}
+
+impl Flag3dSpec {
+    /// The Turek–Hron flag (`[0.25, 0.6] × [0.19, 0.21]`) extruded over
+    /// `z ∈ [z0, z0 + span]`.
+    pub fn turek_hron(span: f64, z0: f64, nx: usize, ny: usize, nz: usize) -> Self {
+        Self {
+            x0: 0.25,
+            x1: 0.6,
+            y0: 0.19,
+            y1: 0.21,
+            z0,
+            z1: z0 + span,
+            nx,
+            ny,
+            nz,
+        }
+    }
+
+    pub fn span(&self) -> f64 {
+        self.z1 - self.z0
+    }
+}
+
+/// How the lateral faces `z = z0, z1` are held.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum LateralFaces {
+    /// Traction-free: the physical 3-D plate.
+    Free,
+    /// `u_z = 0` at every node: exact 2-D plane strain (the 2-D model's
+    /// definition), for the self-consistency gate.
+    PlaneStrain,
+}
+
+/// Which wetted face of the plate a [`SurfaceFace`] lies on.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+pub enum FlagSide {
+    /// `y = y0`
+    Bottom,
+    /// `y = y1`
+    Top,
+    /// `x = x1`
+    Tip,
+    /// `z = z0`
+    SideLow,
+    /// `z = z1`
+    SideHigh,
+}
+
+/// A Quad8 face of the Hex20 mesh: corners counter-clockwise seen from
+/// outside, then mid-edge nodes `(0-1, 1-2, 2-3, 3-0)` — the library's
+/// Quad8 order, oriented so `∂x/∂ξ × ∂x/∂η` points out of the solid.
+#[derive(Debug, Clone, PartialEq)]
+pub struct SurfaceFace {
+    pub side: FlagSide,
+    pub nodes: [NodeId; 8],
+}
+
+/// The structured Hex20 plate.
+#[derive(Debug, Clone)]
+pub struct Flag3d {
+    pub spec: Flag3dSpec,
+    pub mesh: Mesh,
+    lattice: Vec<Option<NodeId>>,
+    dims: [usize; 3],
+}
+
+const GAUSS3: [(f64, f64); 3] = [
+    (-0.774_596_669_241_483_4, 5.0 / 9.0),
+    (0.0, 8.0 / 9.0),
+    (0.774_596_669_241_483_4, 5.0 / 9.0),
+];
+
+/// Quad8 serendipity shape functions and their `(ξ, η)` derivatives, the
+/// library's node order.
+fn quad8(xi: f64, eta: f64) -> ([f64; 8], [[f64; 2]; 8]) {
+    let corners = [(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)];
+    let mut n = [0.0; 8];
+    let mut d = [[0.0; 2]; 8];
+    for (a, &(xa, ya)) in corners.iter().enumerate() {
+        let (p, q) = (1.0 + xa * xi, 1.0 + ya * eta);
+        let r = xa * xi + ya * eta - 1.0;
+        n[a] = 0.25 * p * q * r;
+        d[a][0] = 0.25 * xa * (q * r + p * q);
+        d[a][1] = 0.25 * ya * (p * r + p * q);
+    }
+    // Mid-edge nodes: (0, -1), (1, 0), (0, 1), (-1, 0).
+    n[4] = 0.5 * (1.0 - xi * xi) * (1.0 - eta);
+    d[4] = [-xi * (1.0 - eta), -0.5 * (1.0 - xi * xi)];
+    n[5] = 0.5 * (1.0 + xi) * (1.0 - eta * eta);
+    d[5] = [0.5 * (1.0 - eta * eta), -(1.0 + xi) * eta];
+    n[6] = 0.5 * (1.0 - xi * xi) * (1.0 + eta);
+    d[6] = [-xi * (1.0 + eta), 0.5 * (1.0 - xi * xi)];
+    n[7] = 0.5 * (1.0 - xi) * (1.0 - eta * eta);
+    d[7] = [-0.5 * (1.0 - eta * eta), -(1.0 - xi) * eta];
+    (n, d)
+}
+
+impl Flag3d {
+    /// Build the mesh (material 0 on every element).
+    pub fn build(spec: Flag3dSpec) -> FeaResult<Self> {
+        assert!(
+            spec.nx > 0 && spec.ny > 0 && spec.nz > 0,
+            "element counts must be positive"
+        );
+        let dims = [2 * spec.nx + 1, 2 * spec.ny + 1, 2 * spec.nz + 1];
+        let mut mesh = Mesh::new(3)?;
+        let mut lattice = vec![None; dims[0] * dims[1] * dims[2]];
+        for i in 0..dims[0] {
+            for j in 0..dims[1] {
+                for k in 0..dims[2] {
+                    if (i % 2) + (j % 2) + (k % 2) > 1 {
+                        continue; // not a serendipity node
+                    }
+                    let x = spec.x0 + (spec.x1 - spec.x0) * i as f64 / (dims[0] - 1) as f64;
+                    let y = spec.y0 + (spec.y1 - spec.y0) * j as f64 / (dims[1] - 1) as f64;
+                    let z = spec.z0 + (spec.z1 - spec.z0) * k as f64 / (dims[2] - 1) as f64;
+                    lattice[(i * dims[1] + j) * dims[2] + k] =
+                        Some(mesh.add_node(Node::new_3d(x, y, z)));
+                }
+            }
+        }
+        let mut flag = Self {
+            spec,
+            mesh,
+            lattice,
+            dims,
+        };
+        for ex in 0..spec.nx {
+            for ey in 0..spec.ny {
+                for ez in 0..spec.nz {
+                    let (a, b, c) = (2 * ex, 2 * ey, 2 * ez);
+                    let at = |i, j, k| flag.lattice_node(i, j, k).expect("serendipity node");
+                    let nodes = vec![
+                        at(a, b, c),
+                        at(a + 2, b, c),
+                        at(a + 2, b + 2, c),
+                        at(a, b + 2, c),
+                        at(a, b, c + 2),
+                        at(a + 2, b, c + 2),
+                        at(a + 2, b + 2, c + 2),
+                        at(a, b + 2, c + 2),
+                        at(a + 1, b, c),
+                        at(a + 2, b + 1, c),
+                        at(a + 1, b + 2, c),
+                        at(a, b + 1, c),
+                        at(a + 1, b, c + 2),
+                        at(a + 2, b + 1, c + 2),
+                        at(a + 1, b + 2, c + 2),
+                        at(a, b + 1, c + 2),
+                        at(a, b, c + 1),
+                        at(a + 2, b, c + 1),
+                        at(a + 2, b + 2, c + 1),
+                        at(a, b + 2, c + 1),
+                    ];
+                    flag.mesh.add_element(Element::new(
+                        ElementType::Hex20,
+                        nodes,
+                        MaterialId(0),
+                    )?)?;
+                }
+            }
+        }
+        Ok(flag)
+    }
+
+    /// The node at lattice index `(i, j, k)` (`0..=2n` per direction), if
+    /// that lattice point carries a serendipity node.
+    pub fn lattice_node(&self, i: usize, j: usize, k: usize) -> Option<NodeId> {
+        if i >= self.dims[0] || j >= self.dims[1] || k >= self.dims[2] {
+            return None;
+        }
+        self.lattice[(i * self.dims[1] + j) * self.dims[2] + k]
+    }
+
+    /// Lattice sizes `(2nx+1, 2ny+1, 2nz+1)`.
+    pub fn lattice_dims(&self) -> [usize; 3] {
+        self.dims
+    }
+
+    /// The node nearest to `p` (reference coordinates).
+    pub fn nearest_node(&self, p: Vector3<f64>) -> NodeId {
+        self.mesh
+            .nodes
+            .iter()
+            .min_by(|a, b| {
+                let da = (a.1.position() - p).norm();
+                let db = (b.1.position() - p).norm();
+                da.partial_cmp(&db).unwrap()
+            })
+            .map(|(&id, _)| id)
+            .expect("non-empty mesh")
+    }
+
+    /// The Turek–Hron point A `(x1, (y0+y1)/2)` on the mid-span line
+    /// (the nearest lattice node: exact for even `nz`, else the nearer of
+    /// the two mid-span candidates).
+    pub fn point_a(&self) -> NodeId {
+        let s = &self.spec;
+        self.nearest_node(Vector3::new(s.x1, 0.5 * (s.y0 + s.y1), 0.5 * (s.z0 + s.z1)))
+    }
+
+    /// Nodes on the clamped root face `x = x0`.
+    pub fn root_nodes(&self) -> Vec<NodeId> {
+        let mut out = Vec::new();
+        for j in 0..self.dims[1] {
+            for k in 0..self.dims[2] {
+                if let Some(id) = self.lattice_node(0, j, k) {
+                    out.push(id);
+                }
+            }
+        }
+        out
+    }
+
+    /// The root clamp (`u = 0` on `x = x0`) plus the lateral condition.
+    pub fn clamp_root(&self, lateral: LateralFaces) -> BoundaryConditionSet {
+        let zero = || DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0)));
+        let dirichlet = |nodes: Vec<NodeId>, component| {
+            BoundaryCondition::Dirichlet(DirichletBC {
+                nodes,
+                components: vec![component],
+                condition_type: zero(),
+                time_range: None,
+                ramping_factor: 1.0,
+                gradual_enforcement: false,
+            })
+        };
+        let root = self.root_nodes();
+        let mut set = BoundaryConditionSet::new();
+        for component in [
+            DofComponent::DisplacementX,
+            DofComponent::DisplacementY,
+            DofComponent::DisplacementZ,
+        ] {
+            set.add_condition(dirichlet(root.clone(), component));
+        }
+        if lateral == LateralFaces::PlaneStrain {
+            let all: Vec<NodeId> = self.mesh.nodes.keys().copied().collect();
+            set.add_condition(dirichlet(all, DofComponent::DisplacementZ));
+        }
+        set
+    }
+
+    /// One linear-elastic material (the TL path reads its Lamé pair and
+    /// density).
+    pub fn materials(e: f64, nu: f64, rho: f64) -> MaterialDatabase {
+        let mut db = MaterialDatabase::new();
+        db.add_material(
+            MaterialId(0),
+            LinearElastic::new(e, nu).with_density(rho),
+            None,
+        );
+        db
+    }
+
+    /// The coupled march's structure: total-Lagrangian SVK, Newmark with
+    /// the given `γ` and `β = (γ + ½)²/4`, 60 Newton iterations — the 2-D
+    /// harness's settings. `num_steps` only matters for
+    /// [`NonlinearDynamicAnalysis::run`]; a coupling uses
+    /// [`NonlinearDynamicAnalysis::stepper`].
+    pub fn dynamic_analysis(
+        &self,
+        e: f64,
+        nu: f64,
+        rho: f64,
+        lateral: LateralFaces,
+        dt: f64,
+        num_steps: usize,
+        gamma: f64,
+    ) -> NonlinearDynamicAnalysis {
+        let beta = (gamma + 0.5).powi(2) / 4.0;
+        NonlinearDynamicAnalysis::new(
+            self.mesh.clone(),
+            Self::materials(e, nu, rho),
+            self.clamp_root(lateral),
+            dt,
+            num_steps,
+            AnalysisConfig::default(),
+        )
+        .with_total_lagrangian()
+        .with_convergence_criteria(ConvergenceCriteria {
+            max_iterations: 60,
+            ..ConvergenceCriteria::default()
+        })
+        .with_newmark_parameters(gamma, beta)
+    }
+
+    /// The Quad8 faces of the given sides (the root face is never
+    /// wetted). Order: side, then element index.
+    pub fn surface_faces(&self, sides: &[FlagSide]) -> Vec<SurfaceFace> {
+        let [di, dj, dk] = self.dims;
+        let (ilast, jlast, klast) = (di - 1, dj - 1, dk - 1);
+        let mut faces = Vec::new();
+        // base lattice point + the two in-face axes (ξ, η) as unit steps.
+        let mut push = |side, base: [usize; 3], u: [usize; 3], v: [usize; 3]| {
+            let p = |cu: usize, cv: usize| {
+                let q = [
+                    base[0] + cu * u[0] + cv * v[0],
+                    base[1] + cu * u[1] + cv * v[1],
+                    base[2] + cu * u[2] + cv * v[2],
+                ];
+                self.lattice_node(q[0], q[1], q[2]).expect("face node")
+            };
+            faces.push(SurfaceFace {
+                side,
+                nodes: [
+                    p(0, 0),
+                    p(2, 0),
+                    p(2, 2),
+                    p(0, 2),
+                    p(1, 0),
+                    p(2, 1),
+                    p(1, 2),
+                    p(0, 1),
+                ],
+            });
+        };
+        const X: [usize; 3] = [1, 0, 0];
+        const Y: [usize; 3] = [0, 1, 0];
+        const Z: [usize; 3] = [0, 0, 1];
+        for &side in sides {
+            match side {
+                FlagSide::Bottom => {
+                    for a in (0..ilast).step_by(2) {
+                        for c in (0..klast).step_by(2) {
+                            push(side, [a, 0, c], X, Z); // x × z = −y
+                        }
+                    }
+                }
+                FlagSide::Top => {
+                    for a in (0..ilast).step_by(2) {
+                        for c in (0..klast).step_by(2) {
+                            push(side, [a, jlast, c], Z, X); // z × x = +y
+                        }
+                    }
+                }
+                FlagSide::Tip => {
+                    for b in (0..jlast).step_by(2) {
+                        for c in (0..klast).step_by(2) {
+                            push(side, [ilast, b, c], Y, Z); // y × z = +x
+                        }
+                    }
+                }
+                FlagSide::SideLow => {
+                    for a in (0..ilast).step_by(2) {
+                        for b in (0..jlast).step_by(2) {
+                            push(side, [a, b, 0], Y, X); // y × x = −z
+                        }
+                    }
+                }
+                FlagSide::SideHigh => {
+                    for a in (0..ilast).step_by(2) {
+                        for b in (0..jlast).step_by(2) {
+                            push(side, [a, b, klast], X, Y); // x × y = +z
+                        }
+                    }
+                }
+            }
+        }
+        faces
+    }
+
+    /// Every wetted face: bottom, top, tip, and both lateral faces.
+    pub fn wetted_faces(&self) -> Vec<SurfaceFace> {
+        self.surface_faces(&[
+            FlagSide::Bottom,
+            FlagSide::Top,
+            FlagSide::Tip,
+            FlagSide::SideLow,
+            FlagSide::SideHigh,
+        ])
+    }
+
+    /// Consistent nodal forces `f_a = ∫ N_a t(x, n) da` of a traction
+    /// field over the given faces in the configuration `X + u`
+    /// (`displacement` in the stepper's global DOF numbering, read through
+    /// `node_dofs`; `None` = the reference configuration). `traction`
+    /// receives the current point and the current outward unit normal —
+    /// a pressure `p` is `|x, n| -p(x) * n`. 3 × 3 Gauss per face.
+    /// Returns one entry per touched node, sorted by node id.
+    pub fn face_nodal_forces(
+        &self,
+        faces: &[SurfaceFace],
+        displacement: Option<(&DVector<f64>, &dyn Fn(NodeId) -> Vec<usize>)>,
+        traction: &dyn Fn(Vector3<f64>, Vector3<f64>) -> Vector3<f64>,
+    ) -> Vec<(NodeId, Vector3<f64>)> {
+        let position = |id: NodeId| -> Vector3<f64> {
+            let x = self.mesh.get_node(id).expect("face node").position();
+            match displacement {
+                Some((u, dofs)) => {
+                    let d = dofs(id);
+                    x + Vector3::new(u[d[0]], u[d[1]], u[d[2]])
+                }
+                None => x,
+            }
+        };
+        let mut out: BTreeMap<NodeId, Vector3<f64>> = BTreeMap::new();
+        for face in faces {
+            let xs: Vec<Vector3<f64>> = face.nodes.iter().map(|&id| position(id)).collect();
+            for &(xi, wx) in &GAUSS3 {
+                for &(eta, wy) in &GAUSS3 {
+                    let (n, d) = quad8(xi, eta);
+                    let mut x = Vector3::zeros();
+                    let mut t1 = Vector3::zeros();
+                    let mut t2 = Vector3::zeros();
+                    for a in 0..8 {
+                        x += xs[a] * n[a];
+                        t1 += xs[a] * d[a][0];
+                        t2 += xs[a] * d[a][1];
+                    }
+                    let cross = t1.cross(&t2);
+                    let jac = cross.norm();
+                    let t = traction(x, cross / jac);
+                    let w = wx * wy * jac;
+                    for a in 0..8 {
+                        *out.entry(face.nodes[a]).or_insert_with(Vector3::zeros) += t * (n[a] * w);
+                    }
+                }
+            }
+        }
+        out.into_iter().collect()
+    }
+}

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

@@ -7,6 +7,7 @@
 //! coordinating all lower-level components into complete workflows.
 
 pub mod dynamic_analysis;
+pub mod flag3d;
 pub mod modal_analysis;
 pub mod nonlinear_analysis;
 pub mod nonlinear_dynamic;

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

@@ -0,0 +1,638 @@
+//! R8-b: the Turek–Hron flag as a 3-D Hex20 total-Lagrangian SVK solid
+//! ([`rtx_fea::analysis::flag3d`]).
+//!
+//! Suite (fast, run by default):
+//!
+//! 1. `flag3d_mesh_mass_and_surface_forces` — node/element counts, the
+//!    consistent mass sums to `ρ V`, and the face integrator's totals
+//!    (uniform traction on the top face = `t · L · span`; a uniform
+//!    pressure on bottom + top cancels; a rigid translation of the
+//!    configuration changes nothing).
+//! 2. `plane_strain_3d_reproduces_the_2d_csm1` — CSM1 (static, gravity)
+//!    with `u_z = 0` everywhere reproduces the 2-D 35×2 Quad8 plane-strain
+//!    model to rounding (same Newton, same banded LU).
+//! 3. `plane_strain_3d_reproduces_the_2d_csm3_start` — the first 60
+//!    Newmark steps of CSM3 agree with the 2-D stepper to rounding.
+//!
+//! Instruments (`#[ignore]`, env-driven, write under `FLAG3D_OUT`):
+//!
+//! * `flag3d_csm1_table` — CSM1 tip displacement per configuration.
+//! * `flag3d_csm3_march` — the full CSM3 oscillation, CSV of point A and
+//!   of the tip's lateral corners.
+//! * `flag3d_modes_dump` — the linearised operators (TL tangent at u = 0,
+//!   consistent mass) on the free DOFs, for an outside eigen-solve.
+
+use std::io::Write as _;
+
+use nalgebra::{DVector, Vector3};
+use rtx_fea::analysis::flag3d::{Flag3d, Flag3dSpec, FlagSide, LateralFaces};
+use rtx_fea::analysis::{
+    ConvergenceCriteria, DynamicState, NonlinearDynamicAnalysis, NonlinearDynamicStepper,
+};
+use rtx_fea::assembly::dof_mapping::DofComponent;
+use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
+use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
+use rtx_fea::elements::total_lagrangian::{internal_force_and_tangent, saint_venant_kirchhoff};
+use rtx_fea::elements::{ElementMatrixComputer, StandardFiniteElement};
+use rtx_fea::materials::{LinearElastic, Material as _};
+use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+
+const E_MOD: f64 = 1.4e6;
+const NU: f64 = 0.4;
+/// CSM1/CSM3 density and gravity (FSI2's structure is ρ_s = 1e4).
+const RHO_CSM: f64 = 1000.0;
+const G: f64 = 2.0;
+
+fn env_str(name: &str, default: &str) -> String {
+    std::env::var(name).unwrap_or_else(|_| default.to_string())
+}
+
+fn env_num(name: &str, default: f64) -> f64 {
+    std::env::var(name)
+        .map(|v| v.parse().expect(name))
+        .unwrap_or(default)
+}
+
+/// The 2-D flag, exactly as the FSI2 harness builds it.
+fn quad8_flag(nx: usize, ny: usize) -> Mesh {
+    let (x0, x1, y0, y1) = (0.25, 0.6, 0.19, 0.21);
+    let mut mesh = Mesh::new(2).unwrap();
+    let (lx, ly) = (2 * nx + 1, 2 * ny + 1);
+    let mut grid = vec![vec![None; ly]; lx];
+    for (i, column) in grid.iter_mut().enumerate() {
+        for (j, slot) in column.iter_mut().enumerate() {
+            if i % 2 == 1 && j % 2 == 1 {
+                continue;
+            }
+            let x = x0 + (x1 - x0) * i as f64 / (2 * nx) as f64;
+            let y = y0 + (y1 - y0) * j as f64 / (2 * ny) as f64;
+            *slot = Some(mesh.add_node(Node::new_2d(x, y)));
+        }
+    }
+    for i in 0..nx {
+        for j in 0..ny {
+            let (a, b) = (2 * i, 2 * j);
+            let nodes = vec![
+                grid[a][b].unwrap(),
+                grid[a + 2][b].unwrap(),
+                grid[a + 2][b + 2].unwrap(),
+                grid[a][b + 2].unwrap(),
+                grid[a + 1][b].unwrap(),
+                grid[a + 2][b + 1].unwrap(),
+                grid[a + 1][b + 2].unwrap(),
+                grid[a][b + 1].unwrap(),
+            ];
+            mesh.add_element(Element::new(ElementType::Quad8, nodes, MaterialId(0)).unwrap())
+                .unwrap();
+        }
+    }
+    mesh
+}
+
+fn clamp_2d(mesh: &Mesh) -> BoundaryConditionSet {
+    let clamped: Vec<NodeId> = mesh
+        .nodes
+        .iter()
+        .filter(|(_, node)| (node.position().x - 0.25).abs() < 1e-12)
+        .map(|(&id, _)| id)
+        .collect();
+    let mut set = BoundaryConditionSet::new();
+    for component in [DofComponent::DisplacementX, DofComponent::DisplacementY] {
+        set.add_condition(BoundaryCondition::Dirichlet(DirichletBC {
+            nodes: clamped.clone(),
+            components: vec![component],
+            condition_type: DirichletType::Spatial(SpatialFunction(Box::new(|_| 0.0))),
+            time_range: None,
+            ramping_factor: 1.0,
+            gradual_enforcement: false,
+        }));
+    }
+    set
+}
+
+fn point_2d(mesh: &Mesh, x: f64, y: f64) -> NodeId {
+    mesh.nodes
+        .iter()
+        .find(|(_, n)| (n.position().x - x).abs() < 1e-12 && (n.position().y - y).abs() < 1e-12)
+        .map(|(&id, _)| id)
+        .unwrap()
+}
+
+/// Consistent gravity nodal forces `∫ N_a ρ g dV` (per unit depth in 2-D).
+fn gravity_forces(mesh: &Mesh, rho: f64, g: f64) -> Vec<(NodeId, Vector3<f64>)> {
+    let dim = mesh.spatial_dimension;
+    let mut acc: std::collections::BTreeMap<NodeId, Vector3<f64>> = Default::default();
+    for element in mesh.elements.values() {
+        let coords: Vec<Vector3<f64>> = element
+            .nodes
+            .iter()
+            .map(|id| mesh.get_node(*id).unwrap().position())
+            .collect();
+        let fe = StandardFiniteElement::new(element.element_type, coords.clone());
+        let f = ElementMatrixComputer::compute_body_force_vector(
+            &fe,
+            &coords,
+            &|_| Vector3::new(0.0, -rho * g, 0.0),
+            None,
+        )
+        .unwrap();
+        for (a, id) in element.nodes.iter().enumerate() {
+            let e = acc.entry(*id).or_insert_with(Vector3::zeros);
+            for c in 0..dim {
+                e[c] += f[a * dim + c];
+            }
+        }
+    }
+    acc.into_iter().collect()
+}
+
+/// Newton to rounding: the static comparisons are otherwise limited by
+/// the default 1e-6 stopping rule (whose force scale differs between the
+/// 2-D per-unit-depth and the 3-D per-span loads).
+fn static_criteria() -> ConvergenceCriteria {
+    ConvergenceCriteria {
+        force_tolerance: 1e-12,
+        displacement_tolerance: 1e-14,
+        max_iterations: 60,
+        ..ConvergenceCriteria::default()
+    }
+}
+
+/// Static equilibrium through the dynamic stepper: one "Newmark step" of
+/// `dt = 1e4 s` from `u = v = a = 0` is Newton on
+/// `f_int(u) + M u/(β Δt²) = F` — the static problem up to a mass term
+/// 1e-9 of the stiffness. `load_steps` ramps the load, each step starting
+/// from the previous equilibrium with zero velocity and acceleration.
+fn static_solve<'a>(
+    analysis: &'a NonlinearDynamicAnalysis,
+    forces: &[(NodeId, Vector3<f64>)],
+    load_steps: usize,
+) -> (DVector<f64>, NonlinearDynamicStepper<'a>, usize) {
+    let mut stepper = analysis.stepper().unwrap();
+    let n = stepper.rest_state().unwrap().displacement.len();
+    let mut u = DVector::zeros(n);
+    let mut iterations = 0;
+    for s in 1..=load_steps {
+        let scale = s as f64 / load_steps as f64;
+        let scaled: Vec<_> = forces.iter().map(|(id, f)| (*id, f * scale)).collect();
+        stepper.set_nodal_forces(&scaled);
+        let state = DynamicState {
+            displacement: u.clone(),
+            velocity: DVector::zeros(n),
+            acceleration: DVector::zeros(n),
+        };
+        let (next, it) = stepper.step(&state).unwrap();
+        iterations += it;
+        u = next.displacement;
+    }
+    (u, stepper, iterations)
+}
+
+fn static_2d_csm1(nx: usize, ny: usize) -> (f64, f64, usize) {
+    let mesh = quad8_flag(nx, ny);
+    let a = point_2d(&mesh, 0.6, 0.2);
+    let forces = gravity_forces(&mesh, RHO_CSM, G);
+    let analysis = NonlinearDynamicAnalysis::new(
+        mesh.clone(),
+        Flag3d::materials(E_MOD, NU, RHO_CSM),
+        clamp_2d(&mesh),
+        1e4,
+        1,
+        Default::default(),
+    )
+    .with_total_lagrangian()
+    .with_convergence_criteria(static_criteria());
+    let (u, stepper, it) = static_solve(&analysis, &forces, 5);
+    let d = stepper.node_dofs(a);
+    (u[d[0]], u[d[1]], it)
+}
+
+struct Tip3d {
+    a: Vector3<f64>,
+    /// Point A's line at the two lateral faces (z0, z1).
+    side_low: Vector3<f64>,
+    side_high: Vector3<f64>,
+    iterations: usize,
+    dofs: usize,
+}
+
+fn static_3d_csm1(spec: Flag3dSpec, lateral: LateralFaces, load_steps: usize) -> Tip3d {
+    let flag = Flag3d::build(spec).unwrap();
+    let forces = gravity_forces(&flag.mesh, RHO_CSM, G);
+    let analysis = flag
+        .dynamic_analysis(E_MOD, NU, RHO_CSM, lateral, 1e4, 1, 0.5)
+        .with_convergence_criteria(static_criteria());
+    let (u, stepper, iterations) = static_solve(&analysis, &forces, load_steps);
+    let read = |id: NodeId| {
+        let d = stepper.node_dofs(id);
+        Vector3::new(u[d[0]], u[d[1]], u[d[2]])
+    };
+    let ym = 0.5 * (spec.y0 + spec.y1);
+    Tip3d {
+        a: read(flag.point_a()),
+        side_low: read(flag.nearest_node(Vector3::new(spec.x1, ym, spec.z0))),
+        side_high: read(flag.nearest_node(Vector3::new(spec.x1, ym, spec.z1))),
+        iterations,
+        dofs: u.len(),
+    }
+}
+
+fn rel(a: f64, b: f64) -> f64 {
+    ((a - b) / b).abs()
+}
+
+#[test]
+fn flag3d_mesh_mass_and_surface_forces() {
+    let spec = Flag3dSpec::turek_hron(0.1, -0.05, 7, 2, 3);
+    let flag = Flag3d::build(spec).unwrap();
+    // Serendipity lattice: points with at most one odd index.
+    let [di, dj, dk] = flag.lattice_dims();
+    let mut expected = 0;
+    for i in 0..di {
+        for j in 0..dj {
+            for k in 0..dk {
+                if (i % 2) + (j % 2) + (k % 2) <= 1 {
+                    expected += 1;
+                }
+            }
+        }
+    }
+    assert_eq!(flag.mesh.nodes.len(), expected);
+    assert_eq!(flag.mesh.elements.len(), 7 * 2 * 3);
+    // Consistent mass sums to ρ V; every Jacobian is positive.
+    let mut mass = 0.0;
+    for element in flag.mesh.elements.values() {
+        let coords: Vec<Vector3<f64>> = element
+            .nodes
+            .iter()
+            .map(|id| flag.mesh.get_node(*id).unwrap().position())
+            .collect();
+        let fe = StandardFiniteElement::new(element.element_type, coords.clone());
+        let m =
+            ElementMatrixComputer::compute_consistent_mass_matrix(&fe, &coords, 1e4, None).unwrap();
+        mass += m.matrix.sum();
+    }
+    let volume = 0.35 * 0.02 * 0.1;
+    assert!(
+        rel(mass, 1e4 * volume) < 1e-12,
+        "mass {mass} vs {}",
+        1e4 * volume
+    );
+
+    // Uniform traction on the top face: total = t · L · span.
+    let top = flag.surface_faces(&[FlagSide::Top]);
+    let t0 = Vector3::new(3.0, -2.0, 0.5);
+    let total: Vector3<f64> = flag
+        .face_nodal_forces(&top, None, &|_, _| t0)
+        .iter()
+        .map(|(_, f)| f)
+        .sum();
+    assert!((total - t0 * (0.35 * 0.1)).norm() < 1e-12, "{total:?}");
+    // Normals point out: a pressure p on the top pushes down, on the tip
+    // pushes −x, on the side faces pushes inward; bottom + top cancel.
+    let p = 7.0;
+    let pressure = |_: Vector3<f64>, n: Vector3<f64>| -p * n;
+    let sum = |sides: &[FlagSide]| -> Vector3<f64> {
+        flag.face_nodal_forces(&flag.surface_faces(sides), None, &pressure)
+            .iter()
+            .map(|(_, f)| f)
+            .sum()
+    };
+    assert!((sum(&[FlagSide::Top]) - Vector3::new(0.0, -p * 0.035, 0.0)).norm() < 1e-12);
+    assert!((sum(&[FlagSide::Tip]) - Vector3::new(-p * 0.002, 0.0, 0.0)).norm() < 1e-12);
+    assert!((sum(&[FlagSide::SideHigh]) - Vector3::new(0.0, 0.0, -p * 0.007)).norm() < 1e-12);
+    assert!((sum(&[FlagSide::SideLow]) - Vector3::new(0.0, 0.0, p * 0.007)).norm() < 1e-12);
+    assert!(sum(&[FlagSide::Bottom, FlagSide::Top]).norm() < 1e-12);
+    // All five wetted faces + the root would close; without the root the
+    // pressure resultant is the root's missing +x share.
+    let wetted: Vector3<f64> = flag
+        .face_nodal_forces(&flag.wetted_faces(), None, &pressure)
+        .iter()
+        .map(|(_, f)| f)
+        .sum();
+    assert!((wetted - Vector3::new(-p * 0.002, 0.0, 0.0)).norm() < 1e-12);
+    // A rigid translation of the configuration changes nothing.
+    let analysis = flag.dynamic_analysis(1.4e6, 0.4, 1e4, LateralFaces::Free, 1e-3, 1, 0.5);
+    let stepper = analysis.stepper().unwrap();
+    let mut u = DVector::zeros(3 * flag.mesh.nodes.len());
+    for id in flag.mesh.nodes.keys() {
+        let d = stepper.node_dofs(*id);
+        u[d[0]] = 0.01;
+        u[d[1]] = -0.03;
+        u[d[2]] = 0.02;
+    }
+    let dofs = |id: NodeId| stepper.node_dofs(id);
+    let moved = flag.face_nodal_forces(&top, Some((&u, &dofs)), &|_, n| -p * n);
+    let still = flag.face_nodal_forces(&top, None, &|_, n| -p * n);
+    for ((ia, fa), (ib, fb)) in moved.iter().zip(&still) {
+        assert_eq!(ia, ib);
+        assert!((fa - fb).norm() < 1e-14);
+    }
+}
+
+#[test]
+fn plane_strain_3d_reproduces_the_2d_csm1() {
+    let (ux2, uy2, it2) = static_2d_csm1(35, 2);
+    let tip = static_3d_csm1(
+        Flag3dSpec::turek_hron(0.05, 0.0, 35, 2, 1),
+        LateralFaces::PlaneStrain,
+        5,
+    );
+    println!(
+        "  CSM1 35x2: 2-D Quad8 u(A) = ({ux2:.9e}, {uy2:.9e}) [{it2} Newton]; 3-D Hex20 \
+         35x2x1 plane strain u(A) = ({:.9e}, {:.9e}, {:.2e}) [{} Newton, {} DOFs]; \
+         reference (−7.18777e-3, −66.1029e-3)",
+        tip.a.x, tip.a.y, tip.a.z, tip.iterations, tip.dofs
+    );
+    assert!(rel(tip.a.x, ux2) < 1e-8, "ux {} vs 2-D {ux2}", tip.a.x);
+    assert!(rel(tip.a.y, uy2) < 1e-8, "uy {} vs 2-D {uy2}", tip.a.y);
+    assert!(tip.a.z.abs() < 1e-15);
+    // Span-uniform: both lateral faces carry the mid-span value.
+    assert!((tip.side_low - tip.a).norm() < 1e-9 * tip.a.norm());
+    assert!((tip.side_high - tip.a).norm() < 1e-9 * tip.a.norm());
+    // And the 2-D model is the one pinned against FEATFLOW (1% short in
+    // u_y at 35x2, total_lagrangian_svk.rs).
+    assert!(rel(uy2, -66.1029e-3) < 0.02 && rel(ux2, -7.18777e-3) < 0.04);
+}
+
+fn csm3_2d(dt: f64) -> (NonlinearDynamicAnalysis, NodeId) {
+    let mesh = quad8_flag(35, 2);
+    let a = point_2d(&mesh, 0.6, 0.2);
+    let mut analysis = NonlinearDynamicAnalysis::new(
+        mesh.clone(),
+        Flag3d::materials(E_MOD, NU, RHO_CSM),
+        clamp_2d(&mesh),
+        dt,
+        1,
+        Default::default(),
+    )
+    .with_total_lagrangian();
+    analysis.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0));
+    (analysis, a)
+}
+
+#[test]
+fn plane_strain_3d_reproduces_the_2d_csm3_start() {
+    let dt = 0.005;
+    let steps = 60;
+    let (a2d, node2) = csm3_2d(dt);
+    let mut s2 = a2d.stepper().unwrap();
+    let flag = Flag3d::build(Flag3dSpec::turek_hron(0.05, 0.0, 35, 2, 1)).unwrap();
+    let mut a3d = flag.dynamic_analysis(E_MOD, NU, RHO_CSM, LateralFaces::PlaneStrain, dt, 1, 0.5);
+    a3d.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0));
+    let mut s3 = a3d.stepper().unwrap();
+    let node3 = flag.point_a();
+    let (d2, d3) = (s2.node_dofs(node2), s3.node_dofs(node3));
+    let mut st2 = s2.rest_state().unwrap();
+    let mut st3 = s3.rest_state().unwrap();
+    let mut worst: f64 = 0.0;
+    let mut peak: f64 = 0.0;
+    for _ in 0..steps {
+        st2 = s2.step(&st2).unwrap().0;
+        st3 = s3.step(&st3).unwrap().0;
+        for c in 0..2 {
+            worst = worst.max((st2.displacement[d2[c]] - st3.displacement[d3[c]]).abs());
+            peak = peak.max(st2.displacement[d2[c]].abs());
+        }
+    }
+    println!(
+        "  CSM3 first {steps} steps (t = {:.2} s): max |u_3D − u_2D| at A {worst:.3e} m, \
+         peak |u| {peak:.3e} m",
+        steps as f64 * dt
+    );
+    assert!(peak > 1e-2, "the flag must have moved: {peak}");
+    assert!(
+        worst < 1e-8 * peak,
+        "3-D plane strain departs from 2-D: {worst:.3e}"
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Instruments
+// ---------------------------------------------------------------------------
+
+/// `NXxNYxNZ:span:free|ps` entries, comma-separated.
+fn parse_configs(spec: &str) -> Vec<(usize, usize, usize, f64, LateralFaces)> {
+    spec.split(',')
+        .map(|entry| {
+            let mut parts = entry.trim().split(':');
+            let mesh = parts.next().unwrap();
+            let span: f64 = parts.next().unwrap().parse().unwrap();
+            let lateral = match parts.next().unwrap() {
+                "free" => LateralFaces::Free,
+                "ps" => LateralFaces::PlaneStrain,
+                other => panic!("lateral {other}"),
+            };
+            let n: Vec<usize> = mesh.split('x').map(|t| t.parse().unwrap()).collect();
+            (n[0], n[1], n[2], span, lateral)
+        })
+        .collect()
+}
+
+fn tag(nx: usize, ny: usize, nz: usize, span: f64, lateral: LateralFaces) -> String {
+    let l = if lateral == LateralFaces::Free {
+        "free"
+    } else {
+        "ps"
+    };
+    format!("{nx}x{ny}x{nz}_s{span}_{l}")
+}
+
+#[test]
+#[ignore = "instrument: CSM1 tip displacement per configuration"]
+fn flag3d_csm1_table() {
+    let out = env_str("FLAG3D_OUT", ".");
+    let configs = parse_configs(&env_str(
+        "FLAG3D_CONFIGS",
+        "35x2x1:0.05:ps,35x2x4:0.41:free",
+    ));
+    let steps = env_num("FLAG3D_LOAD_STEPS", 5.0) as usize;
+    let mut table = std::fs::OpenOptions::new()
+        .create(true)
+        .append(true)
+        .open(format!("{out}/csm1_table.txt"))
+        .unwrap();
+    for (nx, ny, nz, span, lateral) in configs {
+        let start = std::time::Instant::now();
+        let tip = static_3d_csm1(
+            Flag3dSpec::turek_hron(span, -0.5 * span, nx, ny, nz),
+            lateral,
+            steps,
+        );
+        let line = format!(
+            "CSM1 {} dofs {} newton {} | A ux {:.6e} uy {:.6e} uz {:.3e} | side_low uy {:.6e} \
+             side_high uy {:.6e} | {:.1} s | ref ux -7.18777e-3 uy -66.1029e-3",
+            tag(nx, ny, nz, span, lateral),
+            tip.dofs,
+            tip.iterations,
+            tip.a.x,
+            tip.a.y,
+            tip.a.z,
+            tip.side_low.y,
+            tip.side_high.y,
+            start.elapsed().as_secs_f64()
+        );
+        println!("{line}");
+        writeln!(table, "{line}").unwrap();
+    }
+    if env_str("FLAG3D_2D", "1") == "1" {
+        for (nx, ny) in [(35, 2), (70, 4)] {
+            let (ux, uy, it) = static_2d_csm1(nx, ny);
+            let line =
+                format!("CSM1 2-D {nx}x{ny} Quad8 | A ux {ux:.6e} uy {uy:.6e} [{it} Newton]");
+            println!("{line}");
+            writeln!(table, "{line}").unwrap();
+        }
+    }
+}
+
+#[test]
+#[ignore = "instrument: the CSM3 oscillation on the 3-D flag"]
+fn flag3d_csm3_march() {
+    let out = env_str("FLAG3D_OUT", ".");
+    let configs = parse_configs(&env_str("FLAG3D_CONFIGS", "35x2x1:0.05:ps"));
+    let dt = env_num("FLAG3D_DT", 0.005);
+    let steps = env_num("FLAG3D_STEPS", 2000.0) as usize;
+    for (nx, ny, nz, span, lateral) in configs {
+        let spec = Flag3dSpec::turek_hron(span, -0.5 * span, nx, ny, nz);
+        let flag = Flag3d::build(spec).unwrap();
+        let mut analysis = flag.dynamic_analysis(E_MOD, NU, RHO_CSM, lateral, dt, steps, 0.5);
+        analysis.set_body_force(|_| Vector3::new(0.0, -RHO_CSM * G, 0.0));
+        let mut stepper = analysis.stepper().unwrap();
+        let ym = 0.5 * (spec.y0 + spec.y1);
+        let probes = [
+            flag.point_a(),
+            flag.nearest_node(Vector3::new(spec.x1, ym, spec.z0)),
+            flag.nearest_node(Vector3::new(spec.x1, ym, spec.z1)),
+        ];
+        let dofs: Vec<Vec<usize>> = probes.iter().map(|p| stepper.node_dofs(*p)).collect();
+        let name = tag(nx, ny, nz, span, lateral);
+        let path = format!("{out}/csm3_{name}_dt{dt}.csv");
+        let mut csv = std::fs::File::create(&path).unwrap();
+        writeln!(csv, "t,ax,ay,az,low_y,high_y,low_z,high_z,newton").unwrap();
+        let mut state = stepper.rest_state().unwrap();
+        let start = std::time::Instant::now();
+        let mut total = 0usize;
+        for k in 1..=steps {
+            let (next, it) = stepper.step(&state).unwrap();
+            state = next;
+            total += it;
+            let u = &state.displacement;
+            writeln!(
+                csv,
+                "{:.6e},{:.12e},{:.12e},{:.12e},{:.12e},{:.12e},{:.12e},{:.12e},{it}",
+                k as f64 * dt,
+                u[dofs[0][0]],
+                u[dofs[0][1]],
+                u[dofs[0][2]],
+                u[dofs[1][1]],
+                u[dofs[2][1]],
+                u[dofs[1][2]],
+                u[dofs[2][2]]
+            )
+            .unwrap();
+        }
+        println!(
+            "CSM3 {name} dt {dt}: {steps} steps, {total} Newton, rescues {:?}, {:.0} s → {path}",
+            stepper.rescue_counts(),
+            start.elapsed().as_secs_f64()
+        );
+    }
+}
+
+#[test]
+#[ignore = "instrument: the 3-D flag's K and M for an outside eigen-solve"]
+fn flag3d_modes_dump() {
+    let out = env_str("FLAG3D_OUT", ".");
+    let configs = parse_configs(&env_str("FLAG3D_CONFIGS", "35x2x1:0.05:ps"));
+    let rho = env_num("FLAG3D_RHO", 1e4);
+    let (lambda, mu) = LinearElastic::new(E_MOD, NU)
+        .with_density(rho)
+        .properties()
+        .lame_parameters();
+    let constitutive = saint_venant_kirchhoff(lambda, mu, 3);
+    for (nx, ny, nz, span, lateral) in configs {
+        let spec = Flag3dSpec::turek_hron(span, -0.5 * span, nx, ny, nz);
+        let flag = Flag3d::build(spec).unwrap();
+        let mut ids: Vec<NodeId> = flag.mesh.nodes.keys().copied().collect();
+        ids.sort();
+        let mut free = std::collections::HashMap::new();
+        let mut dof_lines = Vec::new();
+        for id in &ids {
+            let p = flag.mesh.get_node(*id).unwrap().position();
+            if (p.x - spec.x0).abs() < 1e-12 {
+                continue;
+            }
+            let comps = if lateral == LateralFaces::PlaneStrain {
+                2
+            } else {
+                3
+            };
+            for c in 0..comps {
+                free.insert((*id, c), dof_lines.len());
+                dof_lines.push(format!(
+                    "{} {} {:.12e} {:.12e} {:.12e} {c}",
+                    dof_lines.len(),
+                    id.0,
+                    p.x,
+                    p.y,
+                    p.z
+                ));
+            }
+        }
+        let mut k_trip: std::collections::BTreeMap<(usize, usize), f64> = Default::default();
+        let mut m_trip: std::collections::BTreeMap<(usize, usize), f64> = Default::default();
+        for element in flag.mesh.elements.values() {
+            let coords: Vec<Vector3<f64>> = element
+                .nodes
+                .iter()
+                .map(|id| flag.mesh.get_node(*id).unwrap().position())
+                .collect();
+            let fe = StandardFiniteElement::new(element.element_type, coords.clone());
+            let zero = DVector::zeros(3 * element.nodes.len());
+            let (_, k_e) =
+                internal_force_and_tangent(&fe, &coords, &zero, constitutive.as_ref(), None)
+                    .unwrap();
+            let m_s =
+                ElementMatrixComputer::compute_consistent_mass_matrix(&fe, &coords, rho, None)
+                    .unwrap();
+            let local: Vec<Option<usize>> = element
+                .nodes
+                .iter()
+                .flat_map(|n| (0..3).map(move |c| (*n, c)))
+                .map(|key| free.get(&key).copied())
+                .collect();
+            for (a, ga) in local.iter().enumerate() {
+                let Some(ga) = ga else { continue };
+                for (b, gb) in local.iter().enumerate() {
+                    let Some(gb) = gb else { continue };
+                    *k_trip.entry((*ga, *gb)).or_default() += k_e[(a, b)];
+                    if a % 3 == b % 3 {
+                        *m_trip.entry((*ga, *gb)).or_default() += m_s.matrix[(a / 3, b / 3)];
+                    }
+                }
+            }
+        }
+        let name = tag(nx, ny, nz, span, lateral);
+        let write = |kind: &str, trip: &std::collections::BTreeMap<(usize, usize), f64>| {
+            let mut f = std::fs::File::create(format!("{out}/{kind}_{name}.coo")).unwrap();
+            for ((i, j), v) in trip {
+                if *v != 0.0 {
+                    writeln!(f, "{i} {j} {v:.17e}").unwrap();
+                }
+            }
+        };
+        write("k", &k_trip);
+        write("m", &m_trip);
+        std::fs::write(
+            format!("{out}/dofs_{name}.txt"),
+            dof_lines.join("\n") + "\n",
+        )
+        .unwrap();
+        println!(
+            "modes dump {name}: {} free DOFs, K nnz {}",
+            dof_lines.len(),
+            k_trip.len()
+        );
+    }
+}