rtx-fea: verify elastostatics by manufactured solution — second order confirmed

Adds the Method of Manufactured Solutions to this crate, and it immediately
paid for itself by finding a bug that every existing test missed.

MMS asserts something stronger than "close enough to a value someone
believed": that the discretisation converges to the exact solution at the
rate the theory predicts. Choose a smooth field, substitute it into the
governing equations, and whatever they fail to balance is the body force
that makes it exact. Solve, refine, read off log2(e_h / e_h/2).

Observed order for Quad4 displacement in L2:

    n =  8   L2 error = 7.816681e-3   order    -
    n = 16   L2 error = 1.962845e-3   order 1.99
    n = 32   L2 error = 4.912786e-4   order 2.00
    (n = 64 reads 2.00 as well, at ten times the cost)

That verifies the whole chain at once -- element matrices, quadrature,
Jacobian, assembly, DOF numbering, constraints and the linear solver --
against a solution none of them can represent exactly. It is the check
that none of the sixteen defects fixed in this crate would have survived.

The manufactured field is u = sin(pi x) sin(pi y), v = x^2(1-x) y(1-y):
smooth, not in the bilinear element space, with the two components
different in form and a non-zero shear strain, so the shear block of the
constitutive matrix is exercised rather than skipped.

Found by it: `compute_shape_functions` inferred how many parametric
coordinates to pass from the coordinate *values* --

    match coords.eta() {
        0.0 if coords.zeta() == 0.0 => vec![coords.xi()],   // 1 component
        ...

-- so any evaluation on an axis was handed a one-component slice, which
every 2-D and 3-D element rejects. That includes the element centre and
the middle point of every odd-order Gauss rule. It survived only because
the default 2-point rule never samples zero; asking for a 3-point rule to
integrate the error was enough to trip it. Dimensionality now comes from
the element, which is where it belongs.

Two prerequisites, both real functional gaps rather than test scaffolding:

  - Consistent body-force integration. `BodyForceBC` distributed load as
    force * volume / num_nodes, which is exact only for a constant force
    on a symmetric element and otherwise first-order -- enough to cap the
    measured order of the whole solver at 1 regardless of the element.
    `ElementMatrixComputer::compute_body_force_vector` now integrates
    the consistent form, taking the force as a closure so a spatially
    varying load can be expressed at all.

  - Non-homogeneous Dirichlet conditions did not exist.
    `StaticLinearAnalysis` read the prescribed value out of the boundary
    condition and discarded it, and `extract_free_system` built the
    reduced right-hand side without the K_fc u_c coupling, so every
    Dirichlet condition behaved as zero whatever the caller asked for.
    `GlobalSystem::set_prescribed_value` and the coupling term close
    that, reusing `extract_submatrix` and `multiply_vector` rather than
    a dof-by-dof loop.

560 tests across the three crates, 0 failing.

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

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

@@ -218,6 +218,10 @@ pub struct GlobalSystem {
     pub dof_numbering: DofNumbering,
     /// Assembly options used
     pub options: AssemblyOptions,
+    /// Prescribed values for constrained DOFs, keyed by global DOF index.
+    ///
+    /// A constrained DOF with no entry here is held at zero.
+    prescribed: HashMap<usize, f64>,
 }
 
 impl GlobalSystem {
@@ -232,9 +236,34 @@ impl GlobalSystem {
             force_vector,
             dof_numbering,
             options,
+            prescribed: HashMap::new(),
         }
     }
 
+    /// Prescribe a non-zero value for a constrained degree of freedom.
+    ///
+    /// Support for this was missing entirely: `StaticLinearAnalysis` read the
+    /// value out of a `DirichletBC` and discarded it, and `extract_free_system`
+    /// built the reduced right-hand side without the coupling term below, so
+    /// every Dirichlet condition behaved as if it were zero regardless of what
+    /// the caller asked for.
+    pub fn set_prescribed_value(&mut self, dof: usize, value: f64) -> FeaResult<()> {
+        if dof >= self.dof_numbering.total_dofs {
+            return Err(AssemblyError::GlobalDofOutOfBounds {
+                index: dof,
+                max_index: self.dof_numbering.total_dofs.saturating_sub(1),
+            }
+            .into());
+        }
+        self.prescribed.insert(dof, value);
+        Ok(())
+    }
+
+    /// Prescribed values, keyed by global DOF index.
+    pub fn prescribed_values(&self) -> &HashMap<usize, f64> {
+        &self.prescribed
+    }
+
     /// Add element contribution to the global system.
     pub fn add_element_contribution(
         &mut self,
@@ -348,12 +377,42 @@ impl GlobalSystem {
             .stiffness_matrix
             .extract_submatrix(free_dofs, free_dofs)?;
 
-        // Extract force vector for free DOFs
+        // Extract force vector for free DOFs, carrying the prescribed values
+        // across to the right-hand side.
+        //
+        // Partitioning `K u = F` into free and constrained blocks gives
+        // `K_ff u_f = F_f - K_fc u_c`. The coupling term was absent, which is
+        // correct only when every prescribed value is zero — so a non-zero
+        // Dirichlet condition was silently solved as a zero one.
         let mut free_force = DVector::zeros(n_free);
         for (i, &global_dof) in free_dofs.iter().enumerate() {
             free_force[i] = self.force_vector[global_dof];
         }
 
+        // Carry the prescribed values across to the right-hand side.
+        //
+        // Partitioning `K u = F` into free and constrained blocks gives
+        // `K_ff u_f = F_f - K_fc u_c`. The coupling term was absent, which is
+        // correct only when every prescribed value is zero — so a non-zero
+        // Dirichlet condition was silently solved as a zero one.
+        let nonzero: Vec<(usize, f64)> = self
+            .prescribed
+            .iter()
+            .filter(|&(_, value)| *value != 0.0)
+            .map(|(&dof, &value)| (dof, value))
+            .collect();
+
+        if !nonzero.is_empty() {
+            let columns: Vec<usize> = nonzero.iter().map(|&(dof, _)| dof).collect();
+            let values =
+                DVector::from_iterator(nonzero.len(), nonzero.iter().map(|&(_, value)| value));
+            let coupling = self
+                .stiffness_matrix
+                .extract_submatrix(free_dofs, &columns)?
+                .multiply_vector(&values)?;
+            free_force -= coupling;
+        }
+
         Ok((free_stiffness, free_force))
     }
 

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

@@ -191,6 +191,59 @@ impl ElementMatrixComputer {
         Self::compute_mass_matrix(element, node_coords, density, quadrature_order)
     }
 
+    /// Consistent element load vector for a distributed body force.
+    ///
+    /// Computes `f_i^a = ∫ N_i(x) f_a(x) dV` over the element, with the force
+    /// evaluated at each quadrature point in *physical* coordinates. Returned
+    /// interleaved by node — `[f_0x, f_0y, f_1x, f_1y, …]` — matching the DOF
+    /// layout the stiffness and mass matrices use.
+    ///
+    /// This is the consistent form, not a lumped one. `BodyForceBC` previously
+    /// distributed a body force as `force * volume / num_nodes`, which is
+    /// exact only for a constant force on a symmetric element and is otherwise
+    /// first-order accurate — enough to cap the observed order of accuracy of
+    /// the whole solver at 1 regardless of the element used.
+    ///
+    /// The force is taken as a closure so a spatially varying load can be
+    /// applied, which constant-per-element handling cannot express and which
+    /// the method of manufactured solutions requires.
+    pub fn compute_body_force_vector(
+        element: &dyn FiniteElement,
+        node_coords: &[Vector3<f64>],
+        body_force: &dyn Fn(Vector3<f64>) -> Vector3<f64>,
+        quadrature_order: Option<usize>,
+    ) -> FeaResult<DVector<f64>> {
+        let quad_rule = element.quadrature_rule(quadrature_order)?;
+        let num_nodes = element.num_nodes();
+        let dim = element.spatial_dimension();
+        let mut load = DVector::zeros(num_nodes * dim);
+
+        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 = element.map_to_physical(&point.coords, node_coords)?;
+            let force = body_force(Vector3::new(physical.x(), physical.y(), physical.z()));
+            let scale = jacobian_eval.determinant().abs() * point.weight;
+
+            for i in 0..num_nodes {
+                let n_i = shape_eval.value(i)?;
+                for d in 0..dim {
+                    load[i * dim + d] += n_i * force[d] * scale;
+                }
+            }
+        }
+
+        Ok(load)
+    }
+
     /// Compute lumped mass matrix (diagonal).
     pub fn compute_lumped_mass_matrix(
         element: &dyn FiniteElement,

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

@@ -62,9 +62,19 @@ pub fn compute_shape_functions(
     };
 
     let shape_fn = create_shape_functions(shape_element_type);
-    let xi = match coords.eta() {
+
-        0.0 if coords.zeta() == 0.0 => vec![coords.xi()],
+    // How many parametric coordinates to pass is a property of the *element*,
-        _ if coords.zeta() == 0.0 => vec![coords.xi(), coords.eta()],
+    // not of the values it is being evaluated at.
+    //
+    // This previously inferred the count from the coordinates themselves,
+    // treating `eta == 0.0` as "this is a one-dimensional point". Any
+    // evaluation on an axis was therefore handed a one-component slice, which
+    // every 2-D and 3-D element rejects — including evaluation at the element
+    // centre, and including the middle point of any odd-order Gauss rule. It
+    // survived only because the default 2-point rule never samples zero.
+    let xi = match shape_fn.dimension() {
+        1 => vec![coords.xi()],
+        2 => vec![coords.xi(), coords.eta()],
         _ => vec![coords.xi(), coords.eta(), coords.zeta()],
     };
 

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

@@ -0,0 +1,301 @@
+//! Code verification by the Method of Manufactured Solutions.
+//!
+//! Every other test in this crate asserts that some quantity is *close enough*
+//! to a value someone believed. MMS asserts something stronger and more
+//! objective: that the discretisation converges to the exact solution at the
+//! rate the theory predicts.
+//!
+//! The method needs no benchmark data. Choose any sufficiently smooth field
+//! `u`, substitute it into the governing equations, and whatever they fail to
+//! balance is the body force `f = -∇·σ(u)` that makes `u` the exact solution of
+//! the problem with that load. Solve with `f` applied and `u` imposed on the
+//! boundary, measure the error, refine, and read off the observed order
+//! `log2(e_h / e_{h/2})`.
+//!
+//! This matters here specifically. Every defect found in this crate returned a
+//! plausible wrong answer and passed the tests written for it — element
+//! matrices of zeros, a quadrature rule with no points, an assembly transposed
+//! per node. None of them survives an order-of-accuracy measurement: a wrong
+//! discretisation either fails to converge or converges at the wrong rate, and
+//! neither can be tuned away.
+//!
+//! # The manufactured solution
+//!
+//! ```text
+//!     u(x, y) = sin(pi x) sin(pi y)
+//!     v(x, y) = x^2 (1 - x) y (1 - y)
+//! ```
+//!
+//! Chosen so that the field is smooth but **not** in the element space — a
+//! bilinear `Quad4` cannot represent either component exactly, so the error is
+//! genuinely discretisation error rather than round-off. The two components are
+//! deliberately different in form, and the shear strain
+//! `du/dy + dv/dx` is non-zero, so the shear block of the constitutive matrix
+//! is exercised rather than silently skipped.
+
+use nalgebra::{DVector, Vector3};
+use rtx_fea::assembly::{
+    AdvancedDofNumbering, AssemblyOptions, DofComponent, DofMappingStrategy, GlobalAssembler,
+};
+use rtx_fea::elements::{
+    ElementMatrixComputer, FiniteElement, NaturalCoords, StandardFiniteElement,
+};
+use rtx_fea::materials::{LinearElastic, MaterialDatabase};
+use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+use rtx_fea::solvers::{SolverFactory, SolverOptions};
+use std::f64::consts::PI;
+
+const E: f64 = 1.0;
+const NU: f64 = 0.3;
+
+/// Exact displacement field.
+fn exact(x: f64, y: f64) -> (f64, f64) {
+    (
+        (PI * x).sin() * (PI * y).sin(),
+        x * x * (1.0 - x) * y * (1.0 - y),
+    )
+}
+
+/// Body force `f = -div(sigma(u))` for plane stress.
+///
+/// Derived by hand from the manufactured field above with
+/// `sigma_xx = c(e_xx + nu e_yy)`, `sigma_yy = c(nu e_xx + e_yy)`,
+/// `sigma_xy = mu gamma_xy`, where `c = E / (1 - nu^2)` and
+/// `mu = E / (2(1 + nu))` — the same plane-stress constitutive matrix
+/// `ElementMatrixComputer::elasticity_matrix` builds.
+///
+/// Strains:
+/// ```text
+///     e_xx = pi cos(pi x) sin(pi y)
+///     e_yy = x^2 (1 - x) (1 - 2y)
+///     gamma_xy = pi sin(pi x) cos(pi y) + (2x - 3x^2) y (1 - y)
+/// ```
+fn body_force(x: f64, y: f64) -> (f64, f64) {
+    let c = E / (1.0 - NU * NU);
+    let mu = E / (2.0 * (1.0 + NU));
+
+    let (sx, cx) = ((PI * x).sin(), (PI * x).cos());
+    let (sy, cy) = ((PI * y).sin(), (PI * y).cos());
+
+    // d(e_xx)/dx, d(e_yy)/dx, d(e_xx)/dy, d(e_yy)/dy
+    let dexx_dx = -PI * PI * sx * sy;
+    let deyy_dx = (2.0 * x - 3.0 * x * x) * (1.0 - 2.0 * y);
+    let dexx_dy = PI * PI * cx * cy;
+    let deyy_dy = -2.0 * x * x * (1.0 - x);
+
+    // d(gamma_xy)/dy and d(gamma_xy)/dx
+    let dgxy_dy = -PI * PI * sx * sy + (2.0 * x - 3.0 * x * x) * (1.0 - 2.0 * y);
+    let dgxy_dx = PI * PI * cx * cy + (2.0 - 6.0 * x) * y * (1.0 - y);
+
+    let dsxx_dx = c * (dexx_dx + NU * deyy_dx);
+    let dsyy_dy = c * (NU * dexx_dy + deyy_dy);
+    let dsxy_dy = mu * dgxy_dy;
+    let dsxy_dx = mu * dgxy_dx;
+
+    (-(dsxx_dx + dsxy_dy), -(dsxy_dx + dsyy_dy))
+}
+
+/// Unit square meshed with `n` by `n` `Quad4` elements.
+fn unit_square_mesh(n: usize) -> (Mesh, Vec<Vec<NodeId>>) {
+    let mut mesh = Mesh::new(2).unwrap();
+    let mut grid = vec![vec![NodeId(0); n + 1]; n + 1];
+    for (i, column) in grid.iter_mut().enumerate() {
+        for (j, slot) in column.iter_mut().enumerate() {
+            let x = i as f64 / n as f64;
+            let y = j as f64 / n as f64;
+            *slot = mesh.add_node(Node::new_2d(x, y));
+        }
+    }
+    for i in 0..n {
+        for j in 0..n {
+            let nodes = vec![
+                grid[i][j],
+                grid[i + 1][j],
+                grid[i + 1][j + 1],
+                grid[i][j + 1],
+            ];
+            mesh.add_element(Element::new(ElementType::Quad4, nodes, MaterialId(0)).unwrap())
+                .unwrap();
+        }
+    }
+    (mesh, grid)
+}
+
+/// Solve the manufactured problem on an `n` by `n` mesh and return the
+/// discrete L2 displacement error.
+fn l2_error(n: usize) -> f64 {
+    let (mesh, grid) = unit_square_mesh(n);
+
+    let mut materials = MaterialDatabase::new();
+    materials.add_material(MaterialId(0), LinearElastic::new(E, NU), None);
+
+    let mut dof_numbering =
+        AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
+
+    // Impose the exact solution on the whole boundary.
+    let mut prescribed = Vec::new();
+    for i in 0..=n {
+        for j in 0..=n {
+            if i != 0 && i != n && j != 0 && j != n {
+                continue;
+            }
+            let node = grid[i][j];
+            let (ux, uy) = exact(i as f64 / n as f64, j as f64 / n as f64);
+            for (component, value) in [
+                (DofComponent::DisplacementX, ux),
+                (DofComponent::DisplacementY, uy),
+            ] {
+                let dof = dof_numbering.get_dof(node, component).unwrap();
+                dof_numbering.constrain_dof(dof).unwrap();
+                prescribed.push((dof, value));
+            }
+        }
+    }
+
+    let simple = dof_numbering.to_dof_numbering();
+    let assembler = GlobalAssembler::new(AssemblyOptions::default(), materials);
+    let mut system = assembler.assemble_system(&mesh, &simple, 0.0).unwrap();
+
+    // Consistent body force: integral of N_i f over each element.
+    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 local = ElementMatrixComputer::compute_body_force_vector(
+            &fe,
+            &node_coords,
+            &|p: Vector3<f64>| {
+                let (fx, fy) = body_force(p.x, p.y);
+                Vector3::new(fx, fy, 0.0)
+            },
+            None,
+        )
+        .unwrap();
+
+        let mut global_dofs = Vec::new();
+        for node in &element.nodes {
+            global_dofs.extend_from_slice(simple.get_node_dofs(*node).unwrap());
+        }
+        for (local_index, &global_dof) in global_dofs.iter().enumerate() {
+            system.force_vector[global_dof] += local[local_index];
+        }
+    }
+
+    for (dof, value) in prescribed {
+        system.set_prescribed_value(dof, value).unwrap();
+    }
+
+    let (free_stiffness, free_force) = system.extract_free_system().unwrap();
+    let properties = SolverFactory::analyze_matrix(&free_stiffness);
+    let options = SolverOptions::default();
+    let mut solver = SolverFactory::create_linear_solver(&properties, &options);
+    let (free_solution, _) = solver
+        .solve(&free_stiffness, &free_force, &options)
+        .unwrap();
+
+    // Expand to the full DOF vector, with the prescribed values in place.
+    let mut solution = DVector::zeros(system.dof_numbering.total_dofs);
+    for (i, &dof) in system.dof_numbering.free_dofs.iter().enumerate() {
+        solution[dof] = free_solution[i];
+    }
+    for (&dof, &value) in system.prescribed_values() {
+        solution[dof] = value;
+    }
+
+    // L2 error by the elements' own quadrature.
+    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(Some(3)).unwrap();
+
+        let dofs: Vec<usize> = element
+            .nodes
+            .iter()
+            .flat_map(|node| simple.get_node_dofs(*node).unwrap().to_vec())
+            .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 = (0.0, 0.0);
+            for (node_index, _) in element.nodes.iter().enumerate() {
+                let n = shape.value(node_index).unwrap();
+                uh.0 += n * solution[dofs[node_index * 2]];
+                uh.1 += n * solution[dofs[node_index * 2 + 1]];
+            }
+
+            let (ux, uy) = exact(physical.x(), physical.y());
+            let weight = point.weight * jacobian.determinant().abs();
+            squared += ((uh.0 - ux).powi(2) + (uh.1 - uy).powi(2)) * weight;
+        }
+    }
+
+    squared.sqrt()
+}
+
+/// The displacement error must fall at second order under refinement.
+///
+/// `Quad4` interpolates bilinearly, so the theoretical L2 displacement rate is
+/// 2. A rate near 1 means something in the chain is only first-order accurate;
+/// a rate near 0 means the solution is not converging to the exact one at all,
+/// which is what every stub found in this crate would produce.
+#[test]
+fn quad4_displacement_converges_at_second_order() {
+    // 8 -> 16 -> 32 gives two independent rate estimates, which is enough to
+    // establish the order. A 64 x 64 rung adds a third at roughly ten times
+    // the cost, and it reads 2.00 as well.
+    let resolutions = [8usize, 16, 32];
+    let errors: Vec<f64> = resolutions.iter().map(|&n| l2_error(n)).collect();
+
+    let rates: Vec<f64> = errors
+        .windows(2)
+        .map(|pair| (pair[0] / pair[1]).log2())
+        .collect();
+
+    // Printed so the measurement is auditable rather than merely asserted.
+    for (i, &n) in resolutions.iter().enumerate() {
+        let rate = if i == 0 {
+            String::from("    -")
+        } else {
+            format!("{:5.2}", rates[i - 1])
+        };
+        println!(
+            "  n = {n:3}   L2 error = {:.6e}   observed order = {rate}",
+            errors[i]
+        );
+    }
+
+    for (i, &rate) in rates.iter().enumerate() {
+        assert!(
+            rate > 1.8 && rate < 2.2,
+            "refinement {} -> {}: observed order {rate:.3}, expected 2 for Quad4. \
+             Errors: {errors:?}",
+            resolutions[i],
+            resolutions[i + 1]
+        );
+    }
+}
+
+/// The error must actually be small, not merely converging at the right rate.
+///
+/// A scheme converging at second order toward the wrong answer would satisfy
+/// the rate test alone.
+#[test]
+fn quad4_error_is_small_on_a_fine_mesh() {
+    let error = l2_error(32);
+    assert!(
+        error < 1e-3,
+        "L2 error {error:.3e} on a 32x32 mesh is too large for a second-order \
+         method on a smooth solution"
+    );
+}