rtx-fea: ECSW model-order reduction — POD-Galerkin plus hyper-reduction, verified end to end

The third Farhat gap. New rtx_fea::mor module:

- pod::pod_basis — orthonormal SVD basis with an energy-criterion
  truncation. Verified: rank-2 data yields exactly 2 orthonormal modes that
  reconstruct every snapshot to machine precision; a loose tolerance
  truncates a dominant-mode-plus-noise set to one mode.
- nnls — Lawson-Hanson non-negative least squares with the early stop that
  makes ECSW work: iteration ends at the requested residual, and the
  active-set structure caps the support at one column per outer iteration,
  so sparsity falls out of the stopping tolerance. Verified against KKT
  conditions, exact positive solutions, negative-clipping, and a
  sparsity-vs-tolerance case. Its thresholds are RELATIVE to the problem's
  own scales — the first version used absolute cutoffs (1e-14) that
  silently ended the iteration on ECSW's small-magnitude training systems
  at 1.2e-3 instead of the requested 1e-4.
- ecsw::train_ecsw — element weights such that a small subset reproduces
  the reduced internal force (the virtual work against the basis) over the
  training snapshots. w = 1 solves the system exactly by construction, so
  it is always consistent; nonnegativity is what keeps a sampled element
  from producing energy.
- reduced::ReducedNonlinearModel — Newton in POD coordinates, assembling
  either every element (POD-Galerkin) or the ECSW sample, on the same
  per-element force/tangent machinery the nonlinear analysis uses.

End-to-end verification (tests/ecsw_mor.rs): a clamped nonlinear block,
snapshots from a 4-point load sweep, evaluated at an UNSEEN load factor:

    POD modes: 2         ECSW sample: 5 of 24 elements
    training residual 2.2e-7 (requested 1e-4)
    error vs full solve: POD-Galerkin 3.09e-7, ECSW 3.08e-7
    hyper-reduction cost (ECSW vs full ROM): 1.6e-8

And the assertion with the most teeth: the same 5 elements with their
weights forced to 1 read a relative error of 1.22 — a completely wrong
field — so the accuracy is carried by the WEIGHTS, not by the subset
happening to be representative.

Scope, stated plainly: geometrically linear, materially nonlinear,
homogeneous Dirichlet only (no lifting); the basis lives on the free DOFs.

559 rtx-fea tests, 0 failing.

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

crates/specialized/rtx-fea/src/mor/ecsw.rs

@@ -0,0 +1,94 @@
+//! Energy-conserving sampling and weighting (Farhat, Chapman & Avery 2014).
+//!
+//! A POD-Galerkin reduced model still assembles the internal force over
+//! *every* element, so its cost scales with the full mesh no matter how few
+//! modes it keeps. ECSW replaces the sum over all elements by a weighted sum
+//! over a small subset, with the weights chosen so that the subset
+//! reproduces the *reduced* internal force — the virtual work of the
+//! elements against the basis — on the training snapshots:
+//!
+//! ```text
+//!     min ||C w - b||  s.t.  w >= 0
+//!     C[(s,·), e] = V_e' f_e(u_s)      (element e's reduced force)
+//!     b[(s,·)]    = sum_e V_e' f_e(u_s)
+//! ```
+//!
+//! `w = 1` solves this exactly by construction, so the system is always
+//! consistent; the sparsity comes from stopping the non-negative
+//! least-squares solve as soon as the residual is below `tolerance * ||b||`.
+//! Nonnegativity is what makes the sampled model inherit stability: a
+//! negative weight would let a sampled element *produce* energy.
+
+use super::nnls::nnls;
+use super::reduced::ElementOperator;
+use crate::assembly::dof_mapping::AdvancedDofNumbering;
+use crate::error::{FeaError, FeaResult};
+use crate::materials::MaterialDatabase;
+use crate::mesh::{ElementId, Mesh};
+use nalgebra::{DMatrix, DVector};
+
+/// The trained element sample: which elements carry the reduced internal
+/// force, and with what weights.
+#[derive(Debug, Clone)]
+pub struct EcswModel {
+    /// Elements with strictly positive weight, in mesh iteration order.
+    pub weights: Vec<(ElementId, f64)>,
+    /// `||C w - b|| / ||b||` at the returned weights — the relative error of
+    /// the sampled reduced force over the training set.
+    pub training_residual: f64,
+}
+
+/// Train ECSW weights for `basis` over the given snapshots (full free-DOF
+/// vectors). `tolerance` is relative: the sampled reduced force matches the
+/// exact one to `tolerance * ||b||` over the training set.
+pub fn train_ecsw(
+    mesh: &Mesh,
+    materials: &MaterialDatabase,
+    dof_numbering: &AdvancedDofNumbering,
+    basis: &DMatrix<f64>,
+    snapshots: &[DVector<f64>],
+    tolerance: f64,
+) -> FeaResult<EcswModel> {
+    if snapshots.is_empty() {
+        return Err(FeaError::InvalidInput(
+            "ECSW training needs at least one snapshot".to_string(),
+        ));
+    }
+    let modes = basis.ncols();
+    let operators = ElementOperator::build_all(mesh, materials, dof_numbering, basis)?;
+    let num_elements = operators.len();
+
+    let mut c = DMatrix::zeros(snapshots.len() * modes, num_elements);
+    let mut b = DVector::zeros(snapshots.len() * modes);
+
+    for (s, snapshot) in snapshots.iter().enumerate() {
+        for (e, operator) in operators.iter().enumerate() {
+            let reduced_force =
+                operator.reduced_internal_force(materials, mesh.spatial_dimension, snapshot)?;
+            for m in 0..modes {
+                c[(s * modes + m, e)] = reduced_force[m];
+                b[s * modes + m] += reduced_force[m];
+            }
+        }
+    }
+
+    let scale = b.norm();
+    let weights = nnls(&c, &b, tolerance * scale)?;
+    let training_residual = if scale > 0.0 {
+        (&c * &weights - &b).norm() / scale
+    } else {
+        0.0
+    };
+
+    let selected = operators
+        .iter()
+        .enumerate()
+        .filter(|(e, _)| weights[*e] > 0.0)
+        .map(|(e, operator)| (operator.element_id, weights[e]))
+        .collect();
+
+    Ok(EcswModel {
+        weights: selected,
+        training_residual,
+    })
+}

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

@@ -0,0 +1,29 @@
+//! Model-order reduction: POD-Galerkin projection with ECSW hyper-reduction.
+//!
+//! The pipeline, all offline steps verified by their own invariants:
+//!
+//! 1. Collect full-order solution snapshots (the caller's job — typically a
+//!    parameter or load sweep of [`crate::analysis::NonlinearStaticAnalysis`]).
+//! 2. [`pod::pod_basis`] — orthonormal basis by SVD, truncated at an energy
+//!    criterion.
+//! 3. [`ecsw::train_ecsw`] — nonnegative element weights so that a small
+//!    element subset reproduces the reduced internal force over the
+//!    training set ([`nnls`] with an early stop; sparsity comes from the
+//!    stopping tolerance).
+//! 4. [`reduced::ReducedNonlinearModel`] — Newton in reduced coordinates,
+//!    assembling only the sampled elements.
+//!
+//! Scope: geometrically linear, materially nonlinear (what the nonlinear
+//! analysis supports), homogeneous Dirichlet data. Lifting for
+//! inhomogeneous boundary values is not implemented and the basis is over
+//! the free DOFs only.
+
+pub mod ecsw;
+pub mod nnls;
+pub mod pod;
+pub mod reduced;
+
+pub use ecsw::{EcswModel, train_ecsw};
+pub use nnls::nnls;
+pub use pod::pod_basis;
+pub use reduced::ReducedNonlinearModel;

crates/specialized/rtx-fea/src/mor/nnls.rs

@@ -0,0 +1,180 @@
+//! Non-negative least squares by Lawson–Hanson active sets, with the early
+//! stop that makes ECSW work: iteration ends as soon as the residual is
+//! below the requested tolerance, and the number of positive weights at that
+//! point — not the full-accuracy solution — is what delivers the sparsity.
+
+use crate::error::{FeaError, FeaResult};
+use nalgebra::{DMatrix, DVector};
+
+/// Solve `min ||A w - b||` subject to `w >= 0`, stopping as soon as
+/// `||A w - b|| <= tolerance`. Returns the weight vector (dense storage,
+/// mostly zeros).
+///
+/// The active-set structure guarantees at most one column enters per outer
+/// iteration, so an early stop after `s` iterations has at most `s` nonzero
+/// weights — sparsity falls out of the stopping tolerance.
+pub fn nnls(a: &DMatrix<f64>, b: &DVector<f64>, tolerance: f64) -> FeaResult<DVector<f64>> {
+    let n = a.ncols();
+    if a.nrows() != b.len() {
+        return Err(FeaError::InvalidInput(format!(
+            "nnls: A is {}x{}, b has {} rows",
+            a.nrows(),
+            n,
+            b.len()
+        )));
+    }
+
+    let mut weights = DVector::zeros(n);
+    let mut passive: Vec<usize> = Vec::new();
+    let mut residual = b.clone();
+
+    // Scale for the KKT test: an absolute gradient threshold silently ends
+    // the iteration early on small-magnitude problems (ECSW's training
+    // gradients are ~1e-8 at the start), so the threshold is relative to the
+    // largest column times the current residual.
+    let column_scale = (0..n)
+        .map(|i| a.column(i).norm())
+        .fold(0.0_f64, f64::max)
+        .max(1e-300);
+
+    // Each outer iteration adds one column, so n iterations would already
+    // reproduce plain least squares; 2n allows for removals.
+    for _outer in 0..2 * n {
+        if residual.norm() <= tolerance {
+            break;
+        }
+
+        // Most-descending inactive column.
+        let gradient = a.transpose() * &residual;
+        let candidate = (0..n)
+            .filter(|i| !passive.contains(i))
+            .max_by(|&i, &j| gradient[i].total_cmp(&gradient[j]));
+        let Some(candidate) = candidate else {
+            break;
+        };
+        if gradient[candidate] <= 1e-12 * column_scale * residual.norm() {
+            break; // KKT: no inactive column can meaningfully reduce the residual.
+        }
+        passive.push(candidate);
+
+        // Inner loop: least squares on the passive set, stepping back and
+        // dropping columns whenever a coefficient would go negative.
+        loop {
+            let sub = a.select_columns(passive.iter());
+            let svd = sub.svd(true, true);
+            let cutoff = 1e-12 * svd.singular_values.max();
+            let ls = svd
+                .solve(b, cutoff)
+                .map_err(|e| FeaError::InvalidInput(format!("nnls least squares: {e}")))?;
+
+            if ls.iter().all(|&z| z > 0.0) {
+                for (k, &column) in passive.iter().enumerate() {
+                    weights[column] = ls[k];
+                }
+                break;
+            }
+
+            // Step from w toward z until the first coefficient hits zero.
+            let mut alpha = f64::INFINITY;
+            for (k, &column) in passive.iter().enumerate() {
+                if ls[k] <= 0.0 {
+                    let denominator = weights[column] - ls[k];
+                    if denominator > 0.0 {
+                        alpha = alpha.min(weights[column] / denominator);
+                    }
+                }
+            }
+            if !alpha.is_finite() {
+                // Degenerate: drop the newest column and give up on it.
+                let dropped = passive.pop().expect("passive set is nonempty");
+                weights[dropped] = 0.0;
+                break;
+            }
+            for (k, &column) in passive.iter().enumerate() {
+                weights[column] += alpha * (ls[k] - weights[column]);
+            }
+            passive.retain(|&column| {
+                if weights[column] <= 1e-14 {
+                    weights[column] = 0.0;
+                    false
+                } else {
+                    true
+                }
+            });
+        }
+
+        residual = b - a * &weights;
+    }
+
+    Ok(weights)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Identity system: the constrained answer clips negatives to zero.
+    #[test]
+    fn identity_clips_negatives() {
+        let a = DMatrix::identity(3, 3);
+        let b = DVector::from_vec(vec![2.0, -1.0, 3.0]);
+        let w = nnls(&a, &b, 1e-12).unwrap();
+        assert!((w[0] - 2.0).abs() < 1e-10);
+        assert!(w[1].abs() < 1e-12);
+        assert!((w[2] - 3.0).abs() < 1e-10);
+    }
+
+    /// A target inside the positive cone is met exactly.
+    #[test]
+    fn exact_positive_solution_is_found() {
+        let a = DMatrix::from_row_slice(3, 2, &[1.0, 0.0, 1.0, 1.0, 0.0, 2.0]);
+        let exact = DVector::from_vec(vec![1.5, 0.5]);
+        let b = &a * &exact;
+        let w = nnls(&a, &b, 1e-12).unwrap();
+        assert!((&a * &w - &b).norm() < 1e-10);
+        assert!((w - exact).norm() < 1e-9);
+    }
+
+    /// The unconstrained optimum has a negative component; the constrained
+    /// solution must satisfy the KKT conditions instead: nonnegative weights
+    /// and a gradient that is nonpositive wherever the weight is zero.
+    #[test]
+    fn kkt_holds_when_the_constraint_binds() {
+        let a = DMatrix::from_row_slice(2, 2, &[1.0, 1.0, 0.0, 1.0]);
+        let b = DVector::from_vec(vec![1.0, -1.0]);
+        let w = nnls(&a, &b, 1e-12).unwrap();
+        assert!(w.iter().all(|&x| x >= 0.0));
+        let gradient = a.transpose() * (&b - &a * &w);
+        for i in 0..2 {
+            if w[i] == 0.0 {
+                assert!(
+                    gradient[i] <= 1e-10,
+                    "KKT violated: g[{i}] = {}",
+                    gradient[i]
+                );
+            } else {
+                assert!(gradient[i].abs() < 1e-10);
+            }
+        }
+    }
+
+    /// The early stop trades accuracy for sparsity: with a loose tolerance
+    /// on a system whose columns are redundant, fewer columns enter.
+    #[test]
+    fn early_stop_is_sparse() {
+        // 4 columns, the first alone explains 99% of b.
+        let a = DMatrix::from_row_slice(
+            3,
+            4,
+            &[1.0, 0.9, 0.5, 0.1, 1.0, 1.1, 0.4, 0.2, 1.0, 1.0, 0.6, 0.15],
+        );
+        let b = DVector::from_vec(vec![10.0, 10.05, 10.02]);
+        let loose = nnls(&a, &b, 0.05 * b.norm()).unwrap();
+        let tight = nnls(&a, &b, 1e-10).unwrap();
+        let loose_support = loose.iter().filter(|&&x| x > 0.0).count();
+        let tight_support = tight.iter().filter(|&&x| x > 0.0).count();
+        assert!(loose_support <= tight_support);
+        assert!(loose_support <= 2, "loose support {loose_support}");
+        assert!((&a * &loose - &b).norm() <= 0.05 * b.norm() + 1e-12);
+    }
+}

crates/specialized/rtx-fea/src/mor/pod.rs

@@ -0,0 +1,95 @@
+//! Proper orthogonal decomposition of solution snapshots.
+
+use crate::error::{FeaError, FeaResult};
+use nalgebra::{DMatrix, DVector};
+
+/// Orthonormal POD basis of the snapshot set, truncated at the smallest rank
+/// whose retained singular-value energy `sum(s_i^2)` is at least
+/// `1 - energy_tolerance` of the total.
+pub fn pod_basis(snapshots: &[DVector<f64>], energy_tolerance: f64) -> FeaResult<DMatrix<f64>> {
+    let Some(first) = snapshots.first() else {
+        return Err(FeaError::InvalidInput(
+            "POD needs at least one snapshot".to_string(),
+        ));
+    };
+    let rows = first.len();
+    if snapshots.iter().any(|s| s.len() != rows) {
+        return Err(FeaError::InvalidInput(
+            "POD snapshots have inconsistent lengths".to_string(),
+        ));
+    }
+
+    let mut matrix = DMatrix::zeros(rows, snapshots.len());
+    for (column, snapshot) in snapshots.iter().enumerate() {
+        matrix.set_column(column, snapshot);
+    }
+
+    let svd = matrix.svd(true, false);
+    let u = svd
+        .u
+        .ok_or_else(|| FeaError::InvalidInput("SVD did not produce U".to_string()))?;
+    let singular = &svd.singular_values;
+
+    let total_energy: f64 = singular.iter().map(|s| s * s).sum();
+    if total_energy <= 0.0 {
+        return Err(FeaError::InvalidInput(
+            "POD snapshots are all zero".to_string(),
+        ));
+    }
+    let mut retained = 0.0;
+    let mut rank = 0;
+    for s in singular.iter() {
+        retained += s * s;
+        rank += 1;
+        if retained >= (1.0 - energy_tolerance) * total_energy {
+            break;
+        }
+    }
+
+    Ok(u.columns(0, rank).into_owned())
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Rank-2 data yields exactly 2 modes, orthonormal, that reconstruct the
+    /// snapshots to machine precision.
+    #[test]
+    fn low_rank_data_is_reconstructed_exactly() {
+        let m1 = DVector::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
+        let m2 = DVector::from_vec(vec![5.0, 3.0, 1.0, -1.0, -3.0]);
+        let snapshots: Vec<DVector<f64>> = (0..6)
+            .map(|k| &m1 * (k as f64 + 1.0) + &m2 * (0.5 - 0.3 * k as f64))
+            .collect();
+
+        let basis = pod_basis(&snapshots, 1e-12).unwrap();
+        assert_eq!(basis.ncols(), 2, "rank-2 data must give 2 modes");
+
+        // Orthonormality.
+        let gram = basis.transpose() * &basis;
+        assert!((gram - DMatrix::identity(2, 2)).norm() < 1e-12);
+
+        // Projection reproduces every snapshot.
+        for snapshot in &snapshots {
+            let reconstructed = &basis * (basis.transpose() * snapshot);
+            assert!((reconstructed - snapshot).norm() < 1e-9 * snapshot.norm());
+        }
+    }
+
+    /// A loose energy tolerance truncates: dominant mode plus noise keeps
+    /// one mode.
+    #[test]
+    fn energy_tolerance_truncates() {
+        let dominant = DVector::from_vec(vec![1.0, 1.0, 1.0, 1.0]);
+        let snapshots: Vec<DVector<f64>> = (0..5)
+            .map(|k| {
+                let mut s = &dominant * (10.0 + k as f64);
+                s[k % 4] += 1e-3;
+                s
+            })
+            .collect();
+        let basis = pod_basis(&snapshots, 1e-4).unwrap();
+        assert_eq!(basis.ncols(), 1);
+    }
+}

crates/specialized/rtx-fea/src/mor/reduced.rs

@@ -0,0 +1,237 @@
+//! The reduced nonlinear model: Newton in POD coordinates, with the internal
+//! force assembled either over every element (POD-Galerkin) or over the
+//! ECSW-sampled subset.
+
+use super::ecsw::EcswModel;
+use crate::assembly::dof_mapping::AdvancedDofNumbering;
+use crate::elements::{ElementMatrixComputer, StandardFiniteElement};
+use crate::error::{AnalysisError, FeaError, FeaResult};
+use crate::materials::MaterialDatabase;
+use crate::mesh::{ElementId, Mesh};
+use nalgebra::{DMatrix, DVector, Vector3};
+
+/// One element's contribution to the reduced model: its finite element, its
+/// gather map into the free-DOF vector, and its slice of the basis.
+pub(crate) struct ElementOperator {
+    pub element_id: ElementId,
+    finite_element: StandardFiniteElement,
+    node_coords: Vec<Vector3<f64>>,
+    /// For each element-local DOF, its index in the free-DOF vector, or
+    /// `None` for a constrained DOF (held at zero — the reduced model
+    /// supports homogeneous Dirichlet data only).
+    free_positions: Vec<Option<usize>>,
+    /// Element-local rows of the basis: `(element dofs) x modes`, with zero
+    /// rows for constrained DOFs.
+    local_basis: DMatrix<f64>,
+    material_index: crate::mesh::MaterialId,
+}
+
+impl ElementOperator {
+    pub(crate) fn build_all(
+        mesh: &Mesh,
+        materials: &MaterialDatabase,
+        dof_numbering: &AdvancedDofNumbering,
+        basis: &DMatrix<f64>,
+    ) -> FeaResult<Vec<Self>> {
+        let modes = basis.ncols();
+        let total = dof_numbering.total_dofs;
+        let mut free_index: Vec<Option<usize>> = vec![None; total];
+        for (i, &dof) in dof_numbering.free_dofs.iter().enumerate() {
+            free_index[dof] = Some(i);
+        }
+        if basis.nrows() != dof_numbering.free_dofs.len() {
+            return Err(FeaError::InvalidInput(format!(
+                "basis has {} rows, free-DOF space has {}",
+                basis.nrows(),
+                dof_numbering.free_dofs.len()
+            )));
+        }
+
+        let mut operators = Vec::with_capacity(mesh.num_elements());
+        for (&element_id, element) in &mesh.elements {
+            if materials.get_material(element.material_id).is_none() {
+                return Err(AnalysisError::InvalidConfiguration(format!(
+                    "Material {} not found",
+                    element.material_id.0
+                ))
+                .into());
+            }
+
+            let node_coords: Vec<Vector3<f64>> = element
+                .nodes
+                .iter()
+                .map(|id| mesh.get_node(*id).unwrap().position())
+                .collect();
+            let finite_element =
+                StandardFiniteElement::new(element.element_type, node_coords.clone());
+
+            let dofs: Vec<usize> = element
+                .nodes
+                .iter()
+                .flat_map(|node| dof_numbering.get_node_dofs(*node))
+                .collect();
+            let free_positions: Vec<Option<usize>> =
+                dofs.iter().map(|&dof| free_index[dof]).collect();
+
+            let mut local_basis = DMatrix::zeros(dofs.len(), modes);
+            for (row, position) in free_positions.iter().enumerate() {
+                if let Some(free) = position {
+                    for m in 0..modes {
+                        local_basis[(row, m)] = basis[(*free, m)];
+                    }
+                }
+            }
+
+            operators.push(Self {
+                element_id,
+                finite_element,
+                node_coords,
+                free_positions,
+                local_basis,
+                material_index: element.material_id,
+            });
+        }
+        Ok(operators)
+    }
+
+    fn gather(&self, free_vector: &DVector<f64>) -> DVector<f64> {
+        let mut local = DVector::zeros(self.free_positions.len());
+        for (row, position) in self.free_positions.iter().enumerate() {
+            if let Some(free) = position {
+                local[row] = free_vector[*free];
+            }
+        }
+        local
+    }
+
+    fn force_and_tangent(
+        &self,
+        materials: &MaterialDatabase,
+        local_displacement: &DVector<f64>,
+        spatial_dim: usize,
+    ) -> FeaResult<(DVector<f64>, DMatrix<f64>)> {
+        let material = materials
+            .get_material(self.material_index)
+            .expect("checked at build time");
+        let constitutive = crate::materials::reduced_constitutive(material, spatial_dim)?;
+        ElementMatrixComputer::compute_internal_force_and_tangent(
+            &self.finite_element,
+            &self.node_coords,
+            local_displacement,
+            constitutive.as_ref(),
+            None,
+        )
+    }
+
+    /// `V_e' f_e(u)` for a full free-DOF displacement — used by training.
+    pub(crate) fn reduced_internal_force(
+        &self,
+        materials: &MaterialDatabase,
+        spatial_dim: usize,
+        free_displacement: &DVector<f64>,
+    ) -> FeaResult<DVector<f64>> {
+        let local = self.gather(free_displacement);
+        let (force, _) = self.force_and_tangent(materials, &local, spatial_dim)?;
+        Ok(self.local_basis.transpose() * force)
+    }
+}
+
+/// Reduced nonlinear model over a POD basis, optionally ECSW-sampled.
+pub struct ReducedNonlinearModel<'a> {
+    mesh: &'a Mesh,
+    materials: &'a MaterialDatabase,
+    basis: DMatrix<f64>,
+    operators: Vec<ElementOperator>,
+    /// `(operator index, weight)` — every element at weight 1 for a plain
+    /// POD-Galerkin model, or the ECSW sample.
+    active: Vec<(usize, f64)>,
+}
+
+impl<'a> ReducedNonlinearModel<'a> {
+    /// Build a POD-Galerkin model: every element, weight 1.
+    pub fn new(
+        mesh: &'a Mesh,
+        materials: &'a MaterialDatabase,
+        dof_numbering: &AdvancedDofNumbering,
+        basis: DMatrix<f64>,
+    ) -> FeaResult<Self> {
+        let operators = ElementOperator::build_all(mesh, materials, dof_numbering, &basis)?;
+        let active = (0..operators.len()).map(|e| (e, 1.0)).collect();
+        Ok(Self {
+            mesh,
+            materials,
+            basis,
+            operators,
+            active,
+        })
+    }
+
+    /// Restrict assembly to the ECSW sample.
+    #[must_use]
+    pub fn with_ecsw(mut self, model: &EcswModel) -> Self {
+        let mut active = Vec::with_capacity(model.weights.len());
+        for &(element_id, weight) in &model.weights {
+            if let Some(index) = self
+                .operators
+                .iter()
+                .position(|op| op.element_id == element_id)
+            {
+                active.push((index, weight));
+            }
+        }
+        self.active = active;
+        self
+    }
+
+    /// Number of elements the model assembles over.
+    pub fn active_elements(&self) -> usize {
+        self.active.len()
+    }
+
+    /// Newton in reduced coordinates: solve
+    /// `V' f_int(V q) = V' f_ext` and return the expanded free-DOF
+    /// displacement `V q`.
+    pub fn solve(
+        &self,
+        external_force_free: &DVector<f64>,
+        force_tolerance: f64,
+        max_iterations: usize,
+    ) -> FeaResult<DVector<f64>> {
+        let modes = self.basis.ncols();
+        let spatial_dim = self.mesh.spatial_dimension;
+        let reduced_external = self.basis.transpose() * external_force_free;
+        let scale = reduced_external.norm().max(1.0);
+
+        let mut q = DVector::zeros(modes);
+        for _iteration in 0..max_iterations {
+            let free_displacement = &self.basis * &q;
+
+            let mut residual = reduced_external.clone();
+            let mut jacobian = DMatrix::zeros(modes, modes);
+            for &(index, weight) in &self.active {
+                let operator = &self.operators[index];
+                let local = operator.gather(&free_displacement);
+                let (force, tangent) =
+                    operator.force_and_tangent(self.materials, &local, spatial_dim)?;
+                residual -= operator.local_basis.transpose() * force * weight;
+                jacobian +=
+                    operator.local_basis.transpose() * tangent * &operator.local_basis * weight;
+            }
+
+            if residual.norm() < force_tolerance * scale {
+                return Ok(&self.basis * &q);
+            }
+
+            let delta = jacobian
+                .lu()
+                .solve(&residual)
+                .ok_or_else(|| FeaError::InvalidInput("singular reduced tangent".to_string()))?;
+            q += delta;
+        }
+
+        Err(AnalysisError::ConvergenceFailed {
+            iterations: max_iterations,
+        }
+        .into())
+    }
+}