rtx-fea: wire NonlinearStaticAnalysis — Newton on the consistent tangent, MMS-verified at second order

NonlinearStaticAnalysis::run returned DVector::zeros unconditionally, like
ModalAnalysis and DynamicAnalysis before their repair. It is now full
Newton-Raphson on R(u) = f_ext - f_int(u):

- ElementMatrixComputer::compute_internal_force_and_tangent integrates
  f_int = int(B' sigma dV) and K_T = int(B' D_T B dV) in ONE quadrature
  sweep from a constitutive closure in the element's reduced Voigt space —
  computing both together is what keeps the tangent consistent with the
  stress, which is what quadratic convergence rides on.
- materials::reduced_constitutive bridges the Material trait (Voigt-6) to
  that closure: 3-D passes the total strain straight through; 2-D supports
  the linear plane-stress closed form and refuses nonlinear materials
  explicitly, since plane-stress condensation of a general law needs a
  per-point iteration that is not implemented yet.
- Dirichlet DOFs are held at their (load-scaled) values and Newton runs on
  the free DOFs, so the prescribed motion enters through f_int itself — no
  K_fc bookkeeping to get wrong. Body force enters via set_body_force, the
  same hook pattern the CFD solvers use for manufactured solutions. Uniform
  load stepping; other strategies and quasi-Newton refuse explicitly.
- StandardFiniteElement::compute_internal_forces, previously a zeros stub,
  now delegates to the same machinery.
- Mesh::validate is now called in run() (the old TODO), and NonlinearConfig
  gained a Default.

Verified two ways (tests/nonlinear_static.rs):

- Equivalence: with LinearElastic the loop lands on the directly assembled
  linear solution to 1e-10 in exactly one Newton step — same B, quadrature
  and solver, so any disagreement is the nonlinear assembly.
- Manufactured solution with a genuinely nonlinear material (energy
  W = 1/2 e'De + alpha/3 I1^3, so stress and tangent are exact derivatives;
  body force by central differences of the closed-form stress): L2 errors
  6.032e-2, 1.780e-2, 4.595e-3 on 2/4/8 Hex8 — observed orders 1.76 and
  1.95, climbing to the theoretical 2. The forcing contains the nonlinear
  term, so the order is reachable only if it is solved; an inconsistent
  tangent is caught separately by the iteration-count bound.

This unblocks ECSW model-order reduction, which needs a working nonlinear
solve underneath it. 551 rtx-fea tests, 0 failing.

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

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

@@ -429,6 +429,17 @@ pub struct NonlinearConfig {
     pub convergence_criteria: ConvergenceCriteria,
 }
 
+impl Default for NonlinearConfig {
+    fn default() -> Self {
+        Self {
+            max_load_steps: 1,
+            load_stepping: LoadSteppingStrategy::Uniform,
+            solver_type: NonlinearSolverType::NewtonRaphson,
+            convergence_criteria: ConvergenceCriteria::default(),
+        }
+    }
+}
+
 /// Load stepping strategies.
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
 pub enum LoadSteppingStrategy {

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

@@ -1,17 +1,41 @@
 // Copyright (c) 2024 RustyTorch++ Team
 // Licensed under the Apache License, Version 2.0
 
-//! Nonlinear finite element analysis.
+//! Nonlinear static finite element analysis.
+//!
+//! Full Newton–Raphson on the residual `R(u) = f_ext - f_int(u)`, with the
+//! consistent tangent `K_T(u) = ∂f_int/∂u` assembled per element by
+//! [`ElementMatrixComputer::compute_internal_force_and_tangent`] from the
+//! material's own stress and tangent. Geometrically linear (small strain);
+//! the nonlinearity is the constitutive law.
+//!
+//! Dirichlet data comes from the [`BoundaryConditionSet`]; prescribed DOFs
+//! are held at their (load-scaled) values and the Newton system is reduced to
+//! the free DOFs, so no `K_fc` bookkeeping is needed — the prescribed motion
+//! enters the residual through `f_int(u)` itself. External load enters as a
+//! consistent body-force vector via [`NonlinearStaticAnalysis::set_body_force`],
+//! the same hook pattern the CFD solvers use for manufactured solutions.
+//!
+//! For a *linear* material the loop converges in one Newton step to exactly
+//! the linear solution — `tests/nonlinear_static.rs` pins that equivalence,
+//! and verifies the genuinely nonlinear path by manufactured solution.
+//!
+//! This previously returned `DVector::zeros` unconditionally, like
+//! `ModalAnalysis` and `DynamicAnalysis` before they were wired.
 
-use super::{Analysis, AnalysisConfig, AnalysisResults, NonlinearConfig};
+use super::{Analysis, AnalysisConfig, AnalysisResults, ConvergenceData, NonlinearConfig};
-use crate::boundary::BoundaryConditionSet;
+use crate::analysis::{LoadSteppingStrategy, NonlinearSolverType};
-use crate::error::FeaResult;
+use crate::assembly::SparseMatrix;
-use crate::materials::MaterialDatabase;
+use crate::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
+use crate::boundary::{BoundaryCondition, BoundaryConditionSet};
+use crate::elements::{ElementMatrixComputer, StandardFiniteElement};
+use crate::error::{AnalysisError, FeaResult};
+use crate::materials::{MaterialDatabase, reduced_constitutive};
 use crate::mesh::Mesh;
-use nalgebra::DVector;
+use crate::solvers::{LinearSolver, LuDirect, SolverOptions};
+use nalgebra::{DVector, Vector3};
 
 /// Nonlinear static finite element analysis.
-#[derive(Debug)]
 pub struct NonlinearStaticAnalysis {
     mesh: Mesh,
     materials: MaterialDatabase,
@@ -20,6 +44,20 @@ pub struct NonlinearStaticAnalysis {
     config: AnalysisConfig,
     progress: f64,
     complete: bool,
+    /// Optional body force per unit volume, integrated consistently
+    /// (`∫ N_i f dV`) into the external force vector.
+    #[allow(clippy::type_complexity)]
+    body_force: Option<Box<dyn Fn(Vector3<f64>) -> Vector3<f64> + Send + Sync>>,
+}
+
+impl std::fmt::Debug for NonlinearStaticAnalysis {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.debug_struct("NonlinearStaticAnalysis")
+            .field("mesh", &self.mesh.num_elements())
+            .field("progress", &self.progress)
+            .field("complete", &self.complete)
+            .finish()
+    }
 }
 
 impl NonlinearStaticAnalysis {
@@ -38,23 +76,260 @@ impl NonlinearStaticAnalysis {
             config,
             progress: 0.0,
             complete: false,
+            body_force: None,
         }
     }
+
+    /// Set a body force per unit volume. See the module docs.
+    pub fn set_body_force<F>(&mut self, f: F)
+    where
+        F: Fn(Vector3<f64>) -> Vector3<f64> + Send + Sync + 'static,
+    {
+        self.body_force = Some(Box::new(f));
+    }
+
+    /// Assemble the global internal force and tangent at the current
+    /// displacement, reduced to the free DOFs.
+    fn assemble_reduced(
+        &self,
+        solution: &DVector<f64>,
+        dof_numbering: &AdvancedDofNumbering,
+        free_index: &[Option<usize>],
+        num_free: usize,
+    ) -> FeaResult<(DVector<f64>, SparseMatrix)> {
+        let dim = self.mesh.spatial_dimension;
+        let mut internal_force = DVector::zeros(num_free);
+        let mut tangent = SparseMatrix::new(num_free, num_free);
+
+        for element in self.mesh.elements.values() {
+            let material = self
+                .materials
+                .get_material(element.material_id)
+                .ok_or_else(|| {
+                    AnalysisError::InvalidConfiguration(format!(
+                        "Material {} not found",
+                        element.material_id.0
+                    ))
+                })?;
+            let constitutive = reduced_constitutive(material, dim)?;
+
+            let node_coords: Vec<Vector3<f64>> = element
+                .nodes
+                .iter()
+                .map(|id| self.mesh.get_node(*id).unwrap().position())
+                .collect();
+            let fe = 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 mut element_displacement = DVector::zeros(dofs.len());
+            for (local, &dof) in dofs.iter().enumerate() {
+                element_displacement[local] = solution[dof];
+            }
+
+            let (f_int, k_t) = ElementMatrixComputer::compute_internal_force_and_tangent(
+                &fe,
+                &node_coords,
+                &element_displacement,
+                constitutive.as_ref(),
+                None,
+            )?;
+
+            for (local_row, &dof_row) in dofs.iter().enumerate() {
+                let Some(free_row) = free_index[dof_row] else {
+                    continue;
+                };
+                internal_force[free_row] += f_int[local_row];
+                for (local_col, &dof_col) in dofs.iter().enumerate() {
+                    if let Some(free_col) = free_index[dof_col] {
+                        let value = k_t[(local_row, local_col)];
+                        if value != 0.0 {
+                            tangent.add_entry(free_row, free_col, value)?;
+                        }
+                    }
+                }
+            }
+        }
+
+        tangent.finalize()?;
+        Ok((internal_force, tangent))
+    }
+
+    /// Consistent external force from the body-force field, reduced to the
+    /// free DOFs.
+    fn assemble_external_force(
+        &self,
+        dof_numbering: &AdvancedDofNumbering,
+        free_index: &[Option<usize>],
+        num_free: usize,
+    ) -> FeaResult<DVector<f64>> {
+        let mut external = DVector::zeros(num_free);
+        let Some(force) = self.body_force.as_ref() else {
+            return Ok(external);
+        };
+
+        for element in self.mesh.elements.values() {
+            let node_coords: Vec<Vector3<f64>> = element
+                .nodes
+                .iter()
+                .map(|id| self.mesh.get_node(*id).unwrap().position())
+                .collect();
+            let fe = StandardFiniteElement::new(element.element_type, node_coords.clone());
+            let local = ElementMatrixComputer::compute_body_force_vector(
+                &fe,
+                &node_coords,
+                force.as_ref(),
+                None,
+            )?;
+
+            let dofs: Vec<usize> = element
+                .nodes
+                .iter()
+                .flat_map(|node| dof_numbering.get_node_dofs(*node))
+                .collect();
+            for (local_index, &dof) in dofs.iter().enumerate() {
+                if let Some(free) = free_index[dof] {
+                    external[free] += local[local_index];
+                }
+            }
+        }
+        Ok(external)
+    }
 }
 
 impl Analysis for NonlinearStaticAnalysis {
     fn run(&mut self) -> FeaResult<AnalysisResults> {
-        // Simplified nonlinear analysis
+        self.progress = 0.05;
-        let num_dofs = self.mesh.num_nodes() * 3;
+        self.mesh.validate()?;
-        let solution = DVector::zeros(num_dofs);
+        if !matches!(
+            self.nonlinear_config.solver_type,
+            NonlinearSolverType::NewtonRaphson
+        ) {
+            return Err(AnalysisError::InvalidConfiguration(format!(
+                "nonlinear solver {:?} is not implemented; use NewtonRaphson",
+                self.nonlinear_config.solver_type
+            ))
+            .into());
+        }
+        if !matches!(
+            self.nonlinear_config.load_stepping,
+            LoadSteppingStrategy::Uniform
+        ) {
+            return Err(AnalysisError::InvalidConfiguration(format!(
+                "load stepping {:?} is not implemented; use Uniform",
+                self.nonlinear_config.load_stepping
+            ))
+            .into());
+        }
+
+        // DOF numbering, with Dirichlet DOFs constrained and their full-load
+        // values recorded.
+        let dim = self.mesh.spatial_dimension;
+        let mut dof_numbering =
+            AdvancedDofNumbering::displacement_only(&self.mesh, DofMappingStrategy::Sequential)?;
+        let mut prescribed: Vec<(usize, f64)> = Vec::new();
+        for bc in self.boundary_conditions.conditions() {
+            let BoundaryCondition::Dirichlet(dirichlet) = bc else {
+                continue;
+            };
+            for &node_id in &dirichlet.nodes {
+                let position = self
+                    .mesh
+                    .get_node(node_id)
+                    .ok_or(crate::error::MeshError::NodeNotFound { node_id: node_id.0 })?
+                    .position();
+                let value = dirichlet.get_value(0.0, &position);
+                for component in &dirichlet.components {
+                    if component.canonical_index() >= dim {
+                        continue;
+                    }
+                    if let Some(dof) = dof_numbering.get_dof(node_id, *component) {
+                        dof_numbering.constrain_dof(dof)?;
+                        prescribed.push((dof, value));
+                    }
+                }
+            }
+        }
+
+        let total_dofs = dof_numbering.total_dofs;
+        let mut free_index: Vec<Option<usize>> = vec![None; total_dofs];
+        for (i, &dof) in dof_numbering.free_dofs.iter().enumerate() {
+            free_index[dof] = Some(i);
+        }
+        let num_free = dof_numbering.free_dofs.len();
+
+        let full_external = self.assemble_external_force(&dof_numbering, &free_index, num_free)?;
+        let force_scale = full_external.norm().max(1.0);
+
+        let criteria = &self.nonlinear_config.convergence_criteria;
+        let load_steps = self.nonlinear_config.max_load_steps.max(1);
+        let mut solution = DVector::zeros(total_dofs);
+        let mut solver = LuDirect::new();
+        let solver_options = SolverOptions::default();
+        let mut residual_history = Vec::new();
+        let mut converged = true;
+        let mut total_iterations = 0;
+
+        for step in 1..=load_steps {
+            let load_factor = step as f64 / load_steps as f64;
+            for &(dof, value) in &prescribed {
+                solution[dof] = load_factor * value;
+            }
+            let external = &full_external * load_factor;
+
+            let mut step_converged = false;
+            for _iteration in 0..criteria.max_iterations {
+                let (internal, tangent) =
+                    self.assemble_reduced(&solution, &dof_numbering, &free_index, num_free)?;
+                let residual = &external - &internal;
+                let residual_norm = residual.norm();
+                residual_history.push(residual_norm);
+
+                if residual_norm < criteria.force_tolerance * force_scale {
+                    step_converged = true;
+                    break;
+                }
+
+                total_iterations += 1;
+                let (delta, _) = solver.solve(&tangent, &residual, &solver_options)?;
+                for (i, &dof) in dof_numbering.free_dofs.iter().enumerate() {
+                    solution[dof] += delta[i];
+                }
+
+                if delta.norm() < criteria.displacement_tolerance * solution.norm().max(1.0) {
+                    step_converged = true;
+                    break;
+                }
+            }
+            if !step_converged {
+                converged = false;
+                break;
+            }
+            self.progress = 0.1 + 0.9 * step as f64 / load_steps as f64;
+        }
+
+        let final_residual = residual_history.last().copied().unwrap_or(0.0);
+        let mut results = AnalysisResults::new("Nonlinear Static".to_string(), solution);
+        results.convergence = ConvergenceData {
+            converged,
+            iterations: total_iterations,
+            final_residual,
+            residual_history,
+        };
+
+        if !converged {
+            return Err(AnalysisError::ConvergenceFailed {
+                iterations: total_iterations,
+            }
+            .into());
+        }
 
         self.progress = 1.0;
         self.complete = true;
-
+        Ok(results)
-        Ok(AnalysisResults::new(
-            "Nonlinear Static".to_string(),
-            solution,
-        ))
     }
 
     fn analysis_type(&self) -> &'static str {

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

@@ -179,6 +179,93 @@ impl ElementMatrixComputer {
         ))
     }
 
+    /// Internal force vector and consistent tangent stiffness for a
+    /// (materially) nonlinear element, in one quadrature sweep:
+    ///
+    /// ```text
+    ///     f_int = ∫ Bᵀ σ(ε) dV,      K_T = ∫ Bᵀ D_T(ε) B dV,      ε = B u_e
+    /// ```
+    ///
+    /// `constitutive` maps the strain at a quadrature point to the stress and
+    /// consistent tangent, all in the element's reduced Voigt space — 3
+    /// components in 2-D, 6 in 3-D, engineering shear, matching
+    /// [`Self::compute_stiffness_matrix`]'s `B`. Computing both in the same
+    /// sweep is what keeps them consistent: a tangent evaluated at different
+    /// points than the stress destroys Newton's quadratic convergence
+    /// silently.
+    ///
+    /// For a linear material (`σ = D ε`, `D_T = D`) this reduces to
+    /// `f_int = K u_e` with `K` bit-identical to
+    /// [`Self::compute_stiffness_matrix`] under the same quadrature, which is
+    /// what the nonlinear analysis's equivalence test pins.
+    #[allow(clippy::type_complexity)]
+    pub fn compute_internal_force_and_tangent(
+        element: &dyn FiniteElement,
+        node_coords: &[Vector3<f64>],
+        element_displacement: &DVector<f64>,
+        constitutive: &dyn Fn(&DVector<f64>) -> FeaResult<(DVector<f64>, DMatrix<f64>)>,
+        quadrature_order: Option<usize>,
+    ) -> FeaResult<(DVector<f64>, DMatrix<f64>)> {
+        let quad_rule = element.quadrature_rule(quadrature_order)?;
+        let spatial_dim = element.spatial_dimension();
+        let num_nodes = element.num_nodes();
+        let total_dofs = num_nodes * spatial_dim;
+
+        if element_displacement.len() != total_dofs {
+            return Err(ElementError::MatrixComputationFailed {
+                reason: format!(
+                    "element displacement has {} entries, element has {} DOFs",
+                    element_displacement.len(),
+                    total_dofs
+                ),
+            }
+            .into());
+        }
+
+        let strain_components = if spatial_dim == 2 { 3 } else { 6 };
+        let mut internal_force = DVector::zeros(total_dofs);
+        let mut tangent = DMatrix::zeros(total_dofs, total_dofs);
+
+        for point in &quad_rule.points {
+            let shape_eval = element.shape_functions(&point.coords)?;
+            let jacobian_eval = element.jacobian(&point.coords, node_coords)?;
+            if !jacobian_eval.is_valid() {
+                return Err(ElementError::JacobianSingular {
+                    det: jacobian_eval.determinant,
+                }
+                .into());
+            }
+
+            let physical_derivatives =
+                jacobian_eval.transform_derivatives(&shape_eval.derivatives)?;
+            let b_matrix = Self::strain_displacement_matrix(&physical_derivatives, spatial_dim)?;
+
+            let strain = &b_matrix * element_displacement;
+            let (stress, material_tangent) = constitutive(&strain)?;
+            if stress.len() != strain_components
+                || material_tangent.nrows() != strain_components
+                || material_tangent.ncols() != strain_components
+            {
+                return Err(ElementError::MatrixComputationFailed {
+                    reason: format!(
+                        "constitutive closure returned stress of {} and tangent {}x{}, \
+                         expected {strain_components} components",
+                        stress.len(),
+                        material_tangent.nrows(),
+                        material_tangent.ncols()
+                    ),
+                }
+                .into());
+            }
+
+            let scale = jacobian_eval.determinant().abs() * point.weight;
+            internal_force += b_matrix.transpose() * stress * scale;
+            tangent += b_matrix.transpose() * material_tangent * b_matrix * scale;
+        }
+
+        Ok((internal_force, tangent))
+    }
+
     /// Stiffness of the four-node quadrilateral with incompatible modes —
     /// Wilson's Q6 in Taylor's QM6 form.
     ///

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

@@ -395,18 +395,26 @@ impl StandardFiniteElement {
     }
 
     /// Compute internal forces for an element
+    /// Internal force `∫ Bᵀ σ(ε(u)) dV` for this element at the given element
+    /// displacement, using the material's own stress. This previously
+    /// returned a zero vector unconditionally — the same never-connected
+    /// pattern `compute_element_matrices` had before its repair.
     pub fn compute_internal_forces(
         &self,
-        _material: &dyn crate::materials::Material,
+        material: &dyn crate::materials::Material,
-        _displacement: &nalgebra::DVector<f64>,
+        displacement: &nalgebra::DVector<f64>,
         _time: f64,
     ) -> FeaResult<nalgebra::DVector<f64>> {
-        let num_nodes = self.num_nodes();
+        let constitutive =
-        let dofs_per_node = self.spatial_dimension();
+            crate::materials::reduced_constitutive(material, self.spatial_dimension())?;
-        let total_dofs = num_nodes * dofs_per_node;
+        let (internal_force, _tangent) = ElementMatrixComputer::compute_internal_force_and_tangent(
-
+            self,
-        // Return placeholder force vector
+            &self.node_coordinates,
-        Ok(nalgebra::DVector::zeros(total_dofs))
+            displacement,
+            constitutive.as_ref(),
+            None,
+        )?;
+        Ok(internal_force)
     }
 
     /// Get volume quadrature points

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

@@ -272,6 +272,75 @@ impl MaterialResponse {
 }
 
 /// Base trait for all material models.
+/// Constitutive closure in an element's reduced Voigt space, for the
+/// nonlinear assembly path (`ElementMatrixComputer::
+/// compute_internal_force_and_tangent`).
+///
+/// In 3-D the reduced space *is* the material's full Voigt-6 space, so the
+/// closure passes the total strain straight to
+/// [`Material::compute_response`] (from a virgin state — the analyses using
+/// this are path-independent for now) and hands back its stress and
+/// consistent tangent.
+///
+/// In 2-D the element works in plane stress with 3 strain components, and
+/// condensing a *general* nonlinear material to plane stress requires a
+/// per-point iteration on the out-of-plane strain that is not implemented
+/// yet. A linear material needs no iteration — its plane-stress matrix is
+/// closed-form from `(E, nu)` — so that case is supported and anything else
+/// is an explicit error rather than silently wrong physics.
+#[allow(clippy::type_complexity)]
+pub fn reduced_constitutive(
+    material: &dyn Material,
+    spatial_dim: usize,
+) -> FeaResult<
+    Box<dyn Fn(&nalgebra::DVector<f64>) -> FeaResult<(nalgebra::DVector<f64>, DMatrix<f64>)> + '_>,
+> {
+    match spatial_dim {
+        3 => {
+            let state = material.initialize_state(0.0);
+            Ok(Box::new(move |strain: &nalgebra::DVector<f64>| {
+                let strain6 = Vector6::from_iterator(strain.iter().copied());
+                let response = material.compute_response(&strain6, &state, 0.0)?;
+                if !response.is_valid {
+                    return Err(MaterialError::StateUpdateFailed {
+                        reason: "material reported an invalid response".to_string(),
+                    }
+                    .into());
+                }
+                let stress = nalgebra::DVector::from_iterator(6, response.stress.iter().copied());
+                Ok((stress, response.tangent_matrix))
+            }))
+        }
+        2 => {
+            if !material.is_linear() {
+                return Err(MaterialError::UnsupportedModel {
+                    model: format!(
+                        "plane-stress reduction of nonlinear material '{}'",
+                        material.material_type()
+                    ),
+                }
+                .into());
+            }
+            let properties = material.properties();
+            let (e, nu) = (properties.elastic_modulus, properties.poisson_ratio);
+            let factor = e / (1.0 - nu * nu);
+            let mut d = DMatrix::zeros(3, 3);
+            d[(0, 0)] = factor;
+            d[(1, 1)] = factor;
+            d[(0, 1)] = factor * nu;
+            d[(1, 0)] = factor * nu;
+            d[(2, 2)] = factor * (1.0 - nu) / 2.0;
+            Ok(Box::new(move |strain: &nalgebra::DVector<f64>| {
+                Ok((&d * strain, d.clone()))
+            }))
+        }
+        _ => Err(MaterialError::UnsupportedModel {
+            model: format!("{spatial_dim}-D constitutive reduction"),
+        }
+        .into()),
+    }
+}
+
 pub trait Material: Send + Sync {
     /// Get material properties.
     fn properties(&self) -> &MaterialProperties;

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

@@ -0,0 +1,546 @@
+//! Verification of `NonlinearStaticAnalysis`.
+//!
+//! Two independent instruments, per the crate's verification protocol:
+//!
+//! 1. **Equivalence with the linear path.** With a linear-elastic material
+//!    the Newton loop must land on exactly the solution the linear assembly
+//!    produces — same `B` matrices, same quadrature, same solver — and it
+//!    must get there in one Newton step, because the residual of a linear
+//!    problem after one exact tangent solve is zero. Any disagreement is a
+//!    defect in the nonlinear assembly, since everything else is shared.
+//!
+//! 2. **Manufactured solution with a genuinely nonlinear material.** The
+//!    material `CubicEnergy` below derives from the stored energy
+//!    `W = 1/2 eps' D eps + (alpha/3) I1^3`, so its stress
+//!    `sigma = D eps + alpha I1^2 m` and consistent tangent
+//!    `D_T = D + 2 alpha I1 m m'` (with `m = [1,1,1,0,0,0]'`) are exact by
+//!    construction and the tangent is symmetric. The body force
+//!    `f = -div sigma(eps(u_exact))` is computed by central differences of
+//!    the closed-form stress field, the same trick
+//!    `tests/mms_elastostatics.rs` uses to cross-check its hand-derived
+//!    force. The observed L2 order must be 2 — and it can only get there if
+//!    the nonlinear term is actually solved, because the forcing contains it.
+
+use nalgebra::{DMatrix, DVector, Matrix3, Vector3, Vector6};
+use rtx_fea::analysis::{Analysis, AnalysisConfig, NonlinearConfig, NonlinearStaticAnalysis};
+use rtx_fea::assembly::dof_mapping::{AdvancedDofNumbering, DofComponent, DofMappingStrategy};
+use rtx_fea::boundary::dirichlet::{DirichletBC, DirichletType};
+use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, SpatialFunction};
+use rtx_fea::elements::{ElementMatrixComputer, FiniteElement, StandardFiniteElement};
+use rtx_fea::materials::{
+    LinearElastic, Material, MaterialDatabase, MaterialProperties, MaterialResponse, MaterialState,
+};
+use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+use rtx_fea::solvers::{LinearSolver, LuDirect, SolverOptions};
+
+const E: f64 = 1.0;
+const NU: f64 = 0.3;
+const ALPHA: f64 = 3.0;
+const AMP: f64 = 0.15;
+
+// ---------------------------------------------------------------------------
+// The manufactured field and its gradient, hand-differentiated
+// ---------------------------------------------------------------------------
+
+fn u_exact(p: Vector3<f64>) -> Vector3<f64> {
+    use std::f64::consts::PI;
+    let (x, y, z) = (p.x, p.y, p.z);
+    Vector3::new(
+        AMP * (PI * x).sin() * (PI * y).cos() * (PI * z).cos(),
+        AMP * (PI * x).cos() * (PI * y).sin() * (PI * z).cos(),
+        -2.0 * AMP * (PI * x).cos() * (PI * y).cos() * (PI * z).sin(),
+    )
+}
+
+fn grad_u(p: Vector3<f64>) -> Matrix3<f64> {
+    use std::f64::consts::PI;
+    let (x, y, z) = (p.x, p.y, p.z);
+    let (sx, cx) = ((PI * x).sin(), (PI * x).cos());
+    let (sy, cy) = ((PI * y).sin(), (PI * y).cos());
+    let (sz, cz) = ((PI * z).sin(), (PI * z).cos());
+    let a = AMP * PI;
+    Matrix3::new(
+        a * cx * cy * cz,
+        -a * sx * sy * cz,
+        -a * sx * cy * sz,
+        -a * sx * sy * cz,
+        a * cx * cy * cz,
+        -a * cx * sy * sz,
+        2.0 * a * sx * cy * sz,
+        2.0 * a * cx * sy * sz,
+        -2.0 * a * cx * cy * cz,
+    )
+}
+
+/// Small-strain tensor in Voigt-6 (engineering shear), matching the element
+/// `B` matrix convention.
+fn strain_voigt(p: Vector3<f64>) -> Vector6<f64> {
+    let g = grad_u(p);
+    Vector6::new(
+        g[(0, 0)],
+        g[(1, 1)],
+        g[(2, 2)],
+        g[(1, 2)] + g[(2, 1)],
+        g[(0, 2)] + g[(2, 0)],
+        g[(0, 1)] + g[(1, 0)],
+    )
+}
+
+fn elastic_d() -> DMatrix<f64> {
+    let lambda = E * NU / ((1.0 + NU) * (1.0 - 2.0 * NU));
+    let mu = E / (2.0 * (1.0 + NU));
+    let mut d = DMatrix::zeros(6, 6);
+    for i in 0..3 {
+        for j in 0..3 {
+            d[(i, j)] = lambda;
+        }
+        d[(i, i)] += 2.0 * mu;
+        d[(i + 3, i + 3)] = mu;
+    }
+    d
+}
+
+/// Stress of the `CubicEnergy` material at a point of the exact field.
+fn sigma_exact(p: Vector3<f64>) -> Vector6<f64> {
+    let eps = strain_voigt(p);
+    let d = elastic_d();
+    let i1 = eps[0] + eps[1] + eps[2];
+    let mut sigma = Vector6::zeros();
+    for i in 0..6 {
+        for j in 0..6 {
+            sigma[i] += d[(i, j)] * eps[j];
+        }
+    }
+    for i in 0..3 {
+        sigma[i] += ALPHA * i1 * i1;
+    }
+    sigma
+}
+
+/// Body force `f_i = -d sigma_ij / d x_j`, by central differences of the
+/// closed-form stress. Voigt row (i, j) lookup: the full tensor from Voigt-6.
+fn body_force(p: Vector3<f64>) -> Vector3<f64> {
+    let h = 1e-6;
+    let tensor = |q: Vector3<f64>| -> Matrix3<f64> {
+        let s = sigma_exact(q);
+        Matrix3::new(s[0], s[5], s[4], s[5], s[1], s[3], s[4], s[3], s[2])
+    };
+    let mut f = Vector3::zeros();
+    for j in 0..3 {
+        let mut dq = Vector3::zeros();
+        dq[j] = h;
+        let ds = (tensor(p + dq) - tensor(p - dq)) / (2.0 * h);
+        for i in 0..3 {
+            f[i] -= ds[(i, j)];
+        }
+    }
+    f
+}
+
+// ---------------------------------------------------------------------------
+// The nonlinear test material
+// ---------------------------------------------------------------------------
+
+/// `W = 1/2 eps' D eps + (alpha/3) I1^3`; stress and tangent are exact
+/// derivatives of the energy, so the tangent is consistent by construction.
+struct CubicEnergy {
+    properties: MaterialProperties,
+    d: DMatrix<f64>,
+}
+
+impl CubicEnergy {
+    fn new() -> Self {
+        Self {
+            properties: MaterialProperties::isotropic_elastic(E, NU, 1.0),
+            d: elastic_d(),
+        }
+    }
+}
+
+impl Material for CubicEnergy {
+    fn properties(&self) -> &MaterialProperties {
+        &self.properties
+    }
+
+    fn compute_response(
+        &self,
+        strain: &Vector6<f64>,
+        state: &MaterialState,
+        _dt: f64,
+    ) -> rtx_fea::error::FeaResult<MaterialResponse> {
+        let i1 = strain[0] + strain[1] + strain[2];
+        let mut stress = Vector6::zeros();
+        for i in 0..6 {
+            for j in 0..6 {
+                stress[i] += self.d[(i, j)] * strain[j];
+            }
+        }
+        for i in 0..3 {
+            stress[i] += ALPHA * i1 * i1;
+        }
+
+        let mut tangent = self.d.clone();
+        for i in 0..3 {
+            for j in 0..3 {
+                tangent[(i, j)] += 2.0 * ALPHA * i1;
+            }
+        }
+
+        Ok(MaterialResponse::new(stress, tangent, state.clone()))
+    }
+
+    fn elastic_tangent(&self) -> rtx_fea::error::FeaResult<DMatrix<f64>> {
+        Ok(self.d.clone())
+    }
+
+    fn material_type(&self) -> &'static str {
+        "CubicEnergy"
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Meshing and boundary conditions
+// ---------------------------------------------------------------------------
+
+fn hex8_mesh(n: usize) -> Mesh {
+    let mut mesh = Mesh::new(3).unwrap();
+    let mut grid = vec![vec![vec![NodeId(0); n + 1]; n + 1]; n + 1];
+    for (i, plane) in grid.iter_mut().enumerate() {
+        for (j, column) in plane.iter_mut().enumerate() {
+            for (k, slot) in column.iter_mut().enumerate() {
+                *slot = mesh.add_node(Node::new_3d(
+                    i as f64 / n as f64,
+                    j as f64 / n as f64,
+                    k as f64 / n as f64,
+                ));
+            }
+        }
+    }
+    for i in 0..n {
+        for j in 0..n {
+            for k in 0..n {
+                let nodes = vec![
+                    grid[i][j][k],
+                    grid[i + 1][j][k],
+                    grid[i + 1][j + 1][k],
+                    grid[i][j + 1][k],
+                    grid[i][j][k + 1],
+                    grid[i + 1][j][k + 1],
+                    grid[i + 1][j + 1][k + 1],
+                    grid[i][j + 1][k + 1],
+                ];
+                mesh.add_element(Element::new(ElementType::Hex8, nodes, MaterialId(0)).unwrap())
+                    .unwrap();
+            }
+        }
+    }
+    mesh
+}
+
+fn on_boundary(p: Vector3<f64>) -> bool {
+    (0..3).any(|d| p[d].abs() < 1e-12 || (p[d] - 1.0).abs() < 1e-12)
+}
+
+/// The exact field prescribed on the whole boundary: one spatial Dirichlet
+/// condition per displacement component, over the boundary nodes.
+fn exact_boundary_conditions(mesh: &Mesh) -> BoundaryConditionSet {
+    let boundary_nodes: Vec<NodeId> = mesh
+        .nodes
+        .iter()
+        .filter(|(_, node)| on_boundary(node.position()))
+        .map(|(&id, _)| id)
+        .collect();
+
+    let mut set = BoundaryConditionSet::new();
+    let components = [
+        (DofComponent::DisplacementX, 0usize),
+        (DofComponent::DisplacementY, 1),
+        (DofComponent::DisplacementZ, 2),
+    ];
+    for (component, axis) in components {
+        set.add_condition(BoundaryCondition::Dirichlet(DirichletBC {
+            nodes: boundary_nodes.clone(),
+            components: vec![component],
+            condition_type: DirichletType::Spatial(SpatialFunction(Box::new(move |p| {
+                u_exact(*p)[axis]
+            }))),
+            time_range: None,
+            ramping_factor: 1.0,
+            gradual_enforcement: false,
+        }));
+    }
+    set
+}
+
+/// Quadrature-integrated L2 error of a solved displacement field against the
+/// exact one, mirroring `tests/mms_elastostatics.rs`.
+fn l2_error(mesh: &Mesh, dof_numbering: &AdvancedDofNumbering, solution: &DVector<f64>) -> f64 {
+    let mut squared = 0.0;
+    for element in mesh.elements.values() {
+        let node_coords: Vec<Vector3<f64>> = element
+            .nodes
+            .iter()
+            .map(|id| mesh.get_node(*id).unwrap().position())
+            .collect();
+        let fe = StandardFiniteElement::new(element.element_type, node_coords.clone());
+        let rule = fe.quadrature_rule(None).unwrap();
+        let dofs: Vec<usize> = element
+            .nodes
+            .iter()
+            .flat_map(|node| dof_numbering.get_node_dofs(*node))
+            .collect();
+
+        for point in &rule.points {
+            let shape = fe.shape_functions(&point.coords).unwrap();
+            let jacobian = fe.jacobian(&point.coords, &node_coords).unwrap();
+            let physical = fe.map_to_physical(&point.coords, &node_coords).unwrap();
+
+            let mut uh = Vector3::zeros();
+            for node_index in 0..element.nodes.len() {
+                let n = shape.value(node_index).unwrap();
+                for d in 0..3 {
+                    uh[d] += n * solution[dofs[node_index * 3 + d]];
+                }
+            }
+            let exact = u_exact(physical.coords);
+            squared += (uh - exact).norm_squared() * point.weight * jacobian.determinant().abs();
+        }
+    }
+    squared.sqrt()
+}
+
+fn nonlinear_config() -> NonlinearConfig {
+    NonlinearConfig::default()
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+/// With a linear material, one Newton step must land on exactly the solution
+/// of the directly assembled linear system.
+#[test]
+fn linear_material_reproduces_the_linear_solution_in_one_step() {
+    let mesh = hex8_mesh(3);
+    let bcs = exact_boundary_conditions(&mesh);
+
+    // Nonlinear path.
+    let mut materials = MaterialDatabase::new();
+    materials.add_material(MaterialId(0), LinearElastic::new(E, NU), None);
+    let mut analysis = NonlinearStaticAnalysis::new(
+        mesh.clone(),
+        materials,
+        bcs,
+        nonlinear_config(),
+        AnalysisConfig::default(),
+    );
+    analysis.set_body_force(body_force);
+    let results = analysis.run().unwrap();
+    assert!(results.convergence.converged);
+
+    // A linear problem after one exact tangent solve has zero residual; the
+    // config's single load step should therefore take exactly one iteration.
+    assert!(
+        results.convergence.iterations <= 1,
+        "linear problem took {} Newton iterations",
+        results.convergence.iterations
+    );
+
+    // Direct linear solve with the same numbering, constraints and loads.
+    let mut dof_numbering =
+        AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
+    let mut prescribed = Vec::new();
+    for (&node_id, node) in &mesh.nodes {
+        let p = node.position();
+        if !on_boundary(p) {
+            continue;
+        }
+        let value = u_exact(p);
+        for (i, component) in [
+            DofComponent::DisplacementX,
+            DofComponent::DisplacementY,
+            DofComponent::DisplacementZ,
+        ]
+        .into_iter()
+        .enumerate()
+        {
+            let dof = dof_numbering.get_dof(node_id, component).unwrap();
+            dof_numbering.constrain_dof(dof).unwrap();
+            prescribed.push((dof, value[i]));
+        }
+    }
+    let total = dof_numbering.total_dofs;
+    let mut free_index = vec![None; total];
+    for (i, &dof) in dof_numbering.free_dofs.iter().enumerate() {
+        free_index[dof] = Some(i);
+    }
+    let num_free = dof_numbering.free_dofs.len();
+
+    let mut linear = DVector::zeros(total);
+    for &(dof, value) in &prescribed {
+        linear[dof] = value;
+    }
+
+    // K_ff x = f_f - K_fc u_c, assembled per element.
+    let mut stiffness = rtx_fea::assembly::SparseMatrix::new(num_free, num_free);
+    let mut rhs = DVector::zeros(num_free);
+    for element in mesh.elements.values() {
+        let node_coords: Vec<Vector3<f64>> = element
+            .nodes
+            .iter()
+            .map(|id| mesh.get_node(*id).unwrap().position())
+            .collect();
+        let fe = StandardFiniteElement::new(element.element_type, node_coords.clone());
+        let k = ElementMatrixComputer::compute_stiffness_matrix(&fe, &node_coords, E, NU, None)
+            .unwrap()
+            .matrix;
+        let f =
+            ElementMatrixComputer::compute_body_force_vector(&fe, &node_coords, &body_force, None)
+                .unwrap();
+
+        let dofs: Vec<usize> = element
+            .nodes
+            .iter()
+            .flat_map(|node| dof_numbering.get_node_dofs(*node))
+            .collect();
+        for (row, &dof_row) in dofs.iter().enumerate() {
+            let Some(free_row) = free_index[dof_row] else {
+                continue;
+            };
+            rhs[free_row] += f[row];
+            for (col, &dof_col) in dofs.iter().enumerate() {
+                match free_index[dof_col] {
+                    Some(free_col) => {
+                        if k[(row, col)] != 0.0 {
+                            stiffness
+                                .add_entry(free_row, free_col, k[(row, col)])
+                                .unwrap();
+                        }
+                    }
+                    None => rhs[free_row] -= k[(row, col)] * linear[dof_col],
+                }
+            }
+        }
+    }
+    stiffness.finalize().unwrap();
+    let mut solver = LuDirect::new();
+    let (x, _) = solver
+        .solve(&stiffness, &rhs, &SolverOptions::default())
+        .unwrap();
+    for (i, &dof) in dof_numbering.free_dofs.iter().enumerate() {
+        linear[dof] = x[i];
+    }
+
+    let max_diff = results
+        .displacements
+        .iter()
+        .zip(linear.iter())
+        .map(|(a, b)| (a - b).abs())
+        .fold(0.0_f64, f64::max);
+    assert!(
+        max_diff < 1e-10,
+        "nonlinear path differs from the linear solution by {max_diff:.3e}"
+    );
+}
+
+/// Manufactured solution with the genuinely nonlinear material: the observed
+/// L2 order must be 2, which only happens if the nonlinear term is solved —
+/// the forcing contains it.
+///
+/// Measured (2 -> 4 -> 8 Hex8): L2 error 6.032e-2, 1.780e-2, 4.595e-3 —
+/// observed orders 1.76 and 1.95, climbing to the theoretical 2, with Newton
+/// converging in a handful of iterations at every resolution.
+#[test]
+fn nonlinear_mms_converges_at_second_order() {
+    let resolutions = [2usize, 4, 8];
+    let mut errors = Vec::new();
+
+    for &n in &resolutions {
+        let mesh = hex8_mesh(n);
+        let bcs = exact_boundary_conditions(&mesh);
+        let mut materials = MaterialDatabase::new();
+        materials.add_material(MaterialId(0), CubicEnergy::new(), None);
+
+        let mut analysis = NonlinearStaticAnalysis::new(
+            mesh.clone(),
+            materials,
+            bcs,
+            nonlinear_config(),
+            AnalysisConfig::default(),
+        );
+        analysis.set_body_force(body_force);
+        let results = analysis.run().unwrap();
+        assert!(results.convergence.converged);
+
+        // Full Newton with a consistent tangent converges quadratically: a
+        // handful of iterations per load step, not dozens. A wrong tangent
+        // still creeps to the answer — this is what catches it.
+        let steps = nonlinear_config().max_load_steps.max(1);
+        assert!(
+            results.convergence.iterations <= 8 * steps,
+            "Newton took {} iterations over {steps} load steps — the tangent \
+             is not consistent with the stress",
+            results.convergence.iterations
+        );
+
+        let dof_numbering =
+            AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
+        errors.push(l2_error(&mesh, &dof_numbering, &results.displacements));
+    }
+
+    let rates: Vec<f64> = errors
+        .windows(2)
+        .map(|pair| (pair[0] / pair[1]).log2())
+        .collect();
+    for (i, &n) in resolutions.iter().enumerate() {
+        let rate = if i == 0 {
+            String::from("   -")
+        } else {
+            format!("{:4.2}", rates[i - 1])
+        };
+        println!("  n = {n}   L2 error = {:.6e}   order = {rate}", errors[i]);
+    }
+
+    assert!(errors.windows(2).all(|pair| pair[1] < pair[0]));
+    for &rate in &rates {
+        assert!(
+            (1.7..2.4).contains(&rate),
+            "observed order {rate:.2}, expected 2 for Hex8; errors {errors:?}"
+        );
+    }
+}
+
+/// No load and homogeneous boundary data must produce the zero solution.
+#[test]
+fn zero_problem_stays_zero() {
+    let mesh = hex8_mesh(2);
+    let boundary_nodes: Vec<NodeId> = mesh
+        .nodes
+        .iter()
+        .filter(|(_, node)| on_boundary(node.position()))
+        .map(|(&id, _)| id)
+        .collect();
+    let mut bcs = BoundaryConditionSet::new();
+    bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
+        boundary_nodes,
+        vec![
+            DofComponent::DisplacementX,
+            DofComponent::DisplacementY,
+            DofComponent::DisplacementZ,
+        ],
+        0.0,
+    )));
+
+    let mut materials = MaterialDatabase::new();
+    materials.add_material(MaterialId(0), CubicEnergy::new(), None);
+    let mut analysis = NonlinearStaticAnalysis::new(
+        mesh,
+        materials,
+        bcs,
+        nonlinear_config(),
+        AnalysisConfig::default(),
+    );
+    let results = analysis.run().unwrap();
+    assert!(results.convergence.converged);
+    assert!(results.displacements.norm() < 1e-12);
+}