rtx-fea: repair the element library; the crate is now green with no quarantine

Follows the assembly repair. Takes rtx-fea from 21 failures to 253 passing,
0 failing, 0 ignored, with every `#[ignore]` marker gone.

Shape function bugs, all found by one new test asserting two invariants
across the whole element library at once -- partition of unity, and that
the hand-written derivatives sum to zero. The second is the one that gets
skipped, and it is what caught Hexahedron20.

  - Wedge15 summed to 2 at mid-height. Adding a node on a vertical edge
    contributes L_i (1 - t^2) to the sum, so the two corners sharing that
    edge must each give up half of it; the correction was absent. A
    quadratic wedge that doubles every field interpolated through it.

  - Hexahedron20 had sign errors in four hand-written corner
    derivatives -- nodes 3 and 7 in dN/dr, nodes 1 and 5 in dN/ds. The
    values were correct, so partition of unity passed; only the
    derivative-sum invariant exposed it. The strain computed from this
    element was wrong while its interpolation looked right.

  - Quadrilateral9 emitted its shape functions in raw lexicographic
    lattice order while Quad4 and Quad8 use the standard finite-element
    order. A mesh written the usual way paired each node with the wrong
    basis function, which at the element centre made the Jacobian exactly
    singular.

  - Pyramid13 was not a quadratic pyramid basis: it summed to 4 at the
    element centre, and its `derivatives` allocated a 13x3 matrix then
    wrote rows 13 through 15, having been copied from a sixteen-node
    layout, so it panicked before the wrong values could be used. A
    correct 13-node basis is rational, and there is no pyramid quadrature
    rule to integrate it with, so implementing the basis alone would not
    make the element usable. Both now report the gap explicitly rather
    than panicking. Pyramid5 is unaffected and works.

Fixtures corrected rather than tolerances loosened:

  - von Mises stress of an equal biaxial state expected 0, commented "no
    deviatoric stress". Only a hydrostatic state has that. The correct
    value is 100, and expecting 0 would mean a biaxially loaded sheet
    could never yield. The unequal case expected |100-50|; the von Mises
    stress is not a principal difference.
  - A 3-point Gauss rule was required to integrate sin to 1e-10. No
    correct implementation can. Replaced with a convergence assertion,
    which a wrong rule cannot satisfy by luck.
  - MathUtils::SMALL was asserted below EPSILON * 1000, which inverts the
    relationship a practical zero-threshold needs.
  - The Hex20 Jacobian test put all twelve mid-edge nodes at the origin,
    commented "simplified for test". That is not a hexahedron, and its
    mapping is genuinely singular; it only passed because of the
    derivative sign errors above.
  - ElementFactory was required to build every ElementType including
    Point, which has no interpolation and is deliberately rejected.

MemoryInfo displayed decimal GB while its own test constructed binary
GiB, rendering an 8 GiB device as 8.59. Now GiB throughout.

test_mesh_has_real_algorithms searched the *text* of mesh/mod.rs for the
strings "add_node" and "add_element". It broke when those moved into
submodules, but the real problem is that a source-text search cannot tell
a working function from one returning zeros -- it passed throughout the
period when element matrices were a stub and quadrature returned no
points. Replaced with a test that builds a mesh and checks the result.

The crate doc example imported solvers::DirectSolver and
analysis::StaticAnalysis, neither of which has ever existed, so the
doctest never compiled. Replaced with a modal analysis that runs. Also
dropped the "Production Ready: No mocks, stubs, or TODOs - complete
implementation" line, and replaced it with what is actually validated and
what is not.

rtx-fsi unaffected at 26/26.

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

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

@@ -1,13 +1,25 @@
 // Copyright (c) 2024 RustyTorch++ Team
 // Licensed under the Apache License, Version 2.0
 
-//! Modal analysis for eigenvalue problems.
+//! Modal analysis: natural frequencies and mode shapes.
+//!
+//! Solves the generalized eigenvalue problem `K φ = λ M φ` on the free
+//! degrees of freedom, where `λ = ω²`.
+//!
+//! This previously returned `DVector::zeros(num_modes)` without assembling
+//! anything. Both halves it needs already existed —
+//! [`crate::assembly::GlobalAssembler`] for `K` and `M`, and
+//! [`crate::solvers::eigenvalue::EigenvalueSolver`] for the eigenproblem —
+//! they were simply never connected.
 
 use super::{Analysis, AnalysisConfig, AnalysisData, AnalysisResults};
-use crate::error::FeaResult;
+use crate::assembly::{AdvancedDofNumbering, DofMappingStrategy, GlobalAssembler};
+use crate::boundary::{BoundaryCondition, BoundaryConditionSet};
+use crate::error::{AnalysisError, FeaResult};
 use crate::materials::MaterialDatabase;
 use crate::mesh::Mesh;
-use nalgebra::DVector;
+use crate::solvers::eigenvalue::{EigenvalueSolver, ModalResults};
+use nalgebra::{DMatrix, DVector};
 
 /// Modal analysis for computing natural frequencies and mode shapes.
 #[derive(Debug)]
@@ -16,6 +28,8 @@ pub struct ModalAnalysis {
     materials: MaterialDatabase,
     num_modes: usize,
     config: AnalysisConfig,
+    boundary_conditions: BoundaryConditionSet,
+    shift: Option<f64>,
     progress: f64,
     complete: bool,
 }
@@ -32,27 +46,153 @@ impl ModalAnalysis {
             materials,
             num_modes,
             config,
+            boundary_conditions: BoundaryConditionSet::new(),
+            shift: None,
             progress: 0.0,
             complete: false,
         }
     }
+
+    /// Constrain the structure.
+    ///
+    /// Without this the stiffness matrix is singular: an unconstrained body
+    /// has rigid-body modes at zero frequency. That is a physical fact rather
+    /// than a numerical accident, so [`Self::run`] reports it as an error
+    /// naming the remedy instead of returning six meaningless near-zero
+    /// modes.
+    pub fn with_boundary_conditions(mut self, boundary_conditions: BoundaryConditionSet) -> Self {
+        self.boundary_conditions = boundary_conditions;
+        self
+    }
+
+    /// Solve for the modes nearest this frequency-squared value instead of
+    /// the lowest ones.
+    ///
+    /// The one case that genuinely needs it is a free-free structure, where a
+    /// small negative shift moves the rigid-body modes off the singularity.
+    pub fn with_shift(mut self, shift: f64) -> Self {
+        self.shift = Some(shift);
+        self
+    }
+
+    /// Number the degrees of freedom and apply displacement constraints.
+    fn setup_dof_numbering(&self) -> FeaResult<AdvancedDofNumbering> {
+        let mut dof_numbering = AdvancedDofNumbering::displacement_only(
+            &self.mesh,
+            DofMappingStrategy::BandwidthOptimized,
+        )?;
+
+        for bc in self.boundary_conditions.conditions() {
+            // Only Dirichlet conditions remove degrees of freedom. A modal
+            // analysis has no applied load, so Neumann conditions have
+            // nothing to contribute and are ignored rather than silently
+            // mis-applied.
+            if let BoundaryCondition::Dirichlet(dirichlet) = bc {
+                for &node_id in &dirichlet.nodes {
+                    for &component in &dirichlet.components {
+                        if let Some(dof) = dof_numbering.get_dof(node_id, component) {
+                            dof_numbering.constrain_dof(dof)?;
+                        }
+                    }
+                }
+            }
+        }
+
+        Ok(dof_numbering)
+    }
 }
 
 impl Analysis for ModalAnalysis {
     fn run(&mut self) -> FeaResult<AnalysisResults> {
-        // Simplified modal analysis
+        if self.num_modes == 0 {
-        let num_dofs = self.mesh.num_nodes() * 3;
+            return Err(AnalysisError::InvalidConfiguration(
-        let solution = DVector::zeros(num_dofs);
+                "modal analysis requires at least one mode".to_string(),
+            )
+            .into());
+        }
+        if self.mesh.num_elements() == 0 {
+            return Err(AnalysisError::InvalidConfiguration(
+                "mesh contains no elements".to_string(),
+            )
+            .into());
+        }
+
+        let dof_numbering = self.setup_dof_numbering()?;
+        self.progress = 0.2;
+
+        let assembler =
+            GlobalAssembler::new(self.config.assembly_options.clone(), self.materials.clone());
+        let system =
+            assembler.assemble_system(&self.mesh, &dof_numbering.to_dof_numbering(), 0.0)?;
+        self.progress = 0.5;
+
+        let mass = system.mass_matrix.as_ref().ok_or_else(|| {
+            AnalysisError::InvalidConfiguration(
+                "assembly produced no mass matrix; modal analysis needs one".to_string(),
+            )
+        })?;
+
+        let free_dofs = &system.dof_numbering.free_dofs;
+        if free_dofs.is_empty() {
+            return Err(AnalysisError::InvalidConfiguration(
+                "every degree of freedom is constrained; there is nothing to vibrate".to_string(),
+            )
+            .into());
+        }
+        if self.num_modes > free_dofs.len() {
+            return Err(AnalysisError::InvalidConfiguration(format!(
+                "requested {} modes from a system with {} free degrees of freedom",
+                self.num_modes,
+                free_dofs.len()
+            ))
+            .into());
+        }
+
+        let free_stiffness = system
+            .stiffness_matrix
+            .extract_submatrix(free_dofs, free_dofs)?;
+        let free_mass = mass.extract_submatrix(free_dofs, free_dofs)?;
+        self.progress = 0.7;
+
+        let mut solver = EigenvalueSolver::new(self.num_modes);
+        if let Some(shift) = self.shift {
+            solver = solver.with_shift(shift);
+        }
+        let (eigenvalues, free_shapes) = solver.solve(&free_stiffness, &free_mass)?;
+        self.progress = 0.9;
+
+        // Expand each mode shape from the free-DOF space back to the full
+        // DOF vector, leaving constrained entries at zero.
+        let total_dofs = system.dof_numbering.total_dofs;
+        let mut mode_shapes = DMatrix::zeros(total_dofs, self.num_modes);
+        for mode in 0..self.num_modes {
+            for (i, &global_dof) in free_dofs.iter().enumerate() {
+                mode_shapes[(global_dof, mode)] = free_shapes[(i, mode)];
+            }
+        }
+
+        let modal = ModalResults::from_eigenvalues(eigenvalues, free_shapes);
+
+        // The reported solution is the fundamental mode shape, so that a
+        // caller treating this like any other analysis gets something
+        // physical rather than a zero vector.
+        let solution =
+            DVector::from_iterator(total_dofs, (0..total_dofs).map(|d| mode_shapes[(d, 0)]));
 
         let mut results = AnalysisResults::new("Modal".to_string(), solution);
-
+        results.add_data(
-        // Add modal data
+            "eigenvalues".to_string(),
-        let eigenvalues = DVector::zeros(self.num_modes);
+            AnalysisData::Vector(modal.eigenvalues.clone()),
-        let frequencies: DVector<f64> =
+        );
-            eigenvalues.map(|lambda: f64| (lambda / (2.0 * std::f64::consts::PI)).sqrt());
+        results.add_data(
-
+            "frequencies".to_string(),
-        results.add_data("eigenvalues".to_string(), AnalysisData::Vector(eigenvalues));
+            AnalysisData::Vector(modal.frequencies.clone()),
-        results.add_data("frequencies".to_string(), AnalysisData::Vector(frequencies));
+        );
+        results.add_data(
+            "periods".to_string(),
+            AnalysisData::Vector(modal.periods.clone()),
+        );
+        results.add_data("mode_shapes".to_string(), AnalysisData::Matrix(mode_shapes));
 
         self.progress = 1.0;
         self.complete = true;

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

@@ -21,7 +21,12 @@ pub enum DofMappingStrategy {
 }
 
 /// DOF component types for different physics.
-#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
+///
+/// The declaration order is the **canonical DOF order within a node** — the
+/// order element matrices are built in, so `[u, v, w]` before rotations
+/// before scalar fields. [`DofComponent::canonical_index`] exposes it, and
+/// anything that collects a node's degrees of freedom must sort by it.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
 pub enum DofComponent {
     /// Displacement in X direction
     DisplacementX,
@@ -46,6 +51,26 @@ pub enum DofComponent {
 }
 
 impl DofComponent {
+    /// Position of this component in the canonical within-node DOF order.
+    ///
+    /// Element matrices interleave displacement components per node —
+    /// `[u₀, v₀, u₁, v₁, …]` — so any mapping from a node to its global DOFs
+    /// must present them in this order or the assembly is transposed.
+    pub fn canonical_index(self) -> usize {
+        match self {
+            Self::DisplacementX => 0,
+            Self::DisplacementY => 1,
+            Self::DisplacementZ => 2,
+            Self::RotationX => 3,
+            Self::RotationY => 4,
+            Self::RotationZ => 5,
+            Self::Temperature => 6,
+            Self::Pressure => 7,
+            Self::ElectricPotential => 8,
+            Self::MagneticPotential => 9,
+        }
+    }
+
     /// Get all displacement components.
     pub fn displacement_components() -> Vec<Self> {
         vec![
@@ -115,10 +140,22 @@ impl AdvancedDofNumbering {
     }
 
     /// Create simple displacement-only DOF numbering.
+    ///
+    /// The number of displacement components per node follows the **mesh's
+    /// spatial dimension**: two for a planar mesh, three for a solid one.
+    ///
+    /// This previously always allocated all three components. On a 2-D mesh
+    /// that gave a node three degrees of freedom while the element matrices
+    /// — sized from `ElementType::spatial_dimension()` — supplied two, so
+    /// assembly rejected every contribution with a dimension mismatch. A
+    /// planar mesh has no out-of-plane displacement to number, and numbering
+    /// one produces a global stiffness matrix with an entirely zero row and
+    /// column, which is singular by construction.
     pub fn displacement_only(mesh: &Mesh, strategy: DofMappingStrategy) -> FeaResult<Self> {
         let components = DofComponent::displacement_components();
-        let mut node_components = HashMap::new();
+        let components = components[..mesh.spatial_dimension.min(components.len())].to_vec();
 
+        let mut node_components = HashMap::new();
         for &node_id in mesh.nodes.keys() {
             node_components.insert(node_id, components.clone());
         }
@@ -473,16 +510,46 @@ impl AdvancedDofNumbering {
     }
 
     /// Convert to simplified DofNumbering for compatibility with legacy code.
+    /// Flatten into the simple [`super::DofNumbering`] the assembler consumes.
+    ///
+    /// Each node's degrees of freedom are ordered by
+    /// [`DofComponent`]'s canonical order — `[u, v, w]`, then rotations —
+    /// because `GlobalAssembler` maps element local DOF `k` onto
+    /// `get_node_dofs(node)[k]`, and element matrices are built in exactly
+    /// that order.
+    ///
+    /// This previously pushed the DOFs in `HashMap` iteration order, which is
+    /// unspecified. When it came out as `[v, u]` the assembler wrote the
+    /// element's `u` row into the global `v` row: the resulting matrix was
+    /// still symmetric, still had the right rigid-body null space and still
+    /// summed to the right total mass, but described a structure whose axes
+    /// were transposed per node — and `get_dof(node, DisplacementX)` then
+    /// pointed at the wrong row, so constraints were applied to the wrong
+    /// direction too.
     pub fn to_dof_numbering(&self) -> super::DofNumbering {
         use super::DofNumbering;
-        let mut node_to_dofs: HashMap<NodeId, Vec<usize>> = HashMap::new();
+        let mut node_components: HashMap<NodeId, Vec<(DofComponent, usize)>> = HashMap::new();
         let mut dof_to_node: HashMap<usize, NodeId> = HashMap::new();
 
-        for ((node_id, _component), &dof) in &self.node_component_to_dof {
+        for ((node_id, component), &dof) in &self.node_component_to_dof {
-            node_to_dofs.entry(*node_id).or_default().push(dof);
+            node_components
+                .entry(*node_id)
+                .or_default()
+                .push((*component, dof));
             dof_to_node.insert(dof, *node_id);
         }
 
+        let node_to_dofs: HashMap<NodeId, Vec<usize>> = node_components
+            .into_iter()
+            .map(|(node_id, mut components)| {
+                components.sort_by_key(|(component, _)| component.canonical_index());
+                (
+                    node_id,
+                    components.into_iter().map(|(_, dof)| dof).collect(),
+                )
+            })
+            .collect();
+
         DofNumbering {
             total_dofs: self.total_dofs,
             node_to_dofs,

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

@@ -219,13 +219,40 @@ impl JacobianEval {
     }
 
     /// Transform derivatives from natural to physical coordinates.
+    /// Map shape function derivatives from natural to physical coordinates.
+    ///
+    /// `natural_derivatives` is `(num_nodes × parametric_dim)` — the layout
+    /// [`crate::elements::ShapeFunctionEval`] produces — and the result is
+    /// `(num_nodes × spatial_dim)`, which is what
+    /// [`ElementMatrixComputer::strain_displacement_matrix`] indexes as
+    /// `derivatives[(node, direction)]`.
+    ///
+    /// By the chain rule `∂Nᵢ/∂x_k = Σⱼ (∂Nᵢ/∂ξⱼ)(∂ξⱼ/∂x_k)`, and
+    /// `J⁻¹[(j, k)] = ∂ξⱼ/∂x_k`, so the product is
+    /// `natural_derivatives · J⁻¹` — not `J⁻ᵀ · natural_derivatives`, which
+    /// this previously computed. The two agree only when both matrices are
+    /// square and symmetric; for any element with more nodes than parametric
+    /// directions — that is, every element — the old form was a dimension
+    /// mismatch and panicked inside the BLAS call.
     pub fn transform_derivatives(
         &self,
         natural_derivatives: &DMatrix<f64>,
     ) -> FeaResult<DMatrix<f64>> {
-        // Physical derivatives = inverse_jacobian^T * natural_derivatives
+        if natural_derivatives.ncols() != self.inverse_jacobian.nrows() {
-        let physical_derivatives = self.inverse_jacobian.transpose() * natural_derivatives;
+            return Err(ElementError::InvalidGeometry {
-        Ok(physical_derivatives)
+                message: format!(
+                    "shape derivative layout is ({}x{}), expected {} columns to \
+                     match the {}-dimensional parametric space",
+                    natural_derivatives.nrows(),
+                    natural_derivatives.ncols(),
+                    self.inverse_jacobian.nrows(),
+                    self.inverse_jacobian.nrows(),
+                ),
+            }
+            .into());
+        }
+
+        Ok(natural_derivatives * &self.inverse_jacobian)
     }
 
     /// Get the Jacobian determinant (used for integration).

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

@@ -321,25 +321,76 @@ impl StandardFiniteElement {
         }
     }
 
-    /// Compute element matrices (stiffness, mass, force)
+    /// Compute the element stiffness, mass and force matrices.
+    ///
+    /// This is the function [`crate::assembly::GlobalAssembler`] calls for
+    /// every element in the mesh, so every global matrix in the crate is
+    /// downstream of it. It integrates `K = ∫ Bᵀ D B dV` and
+    /// `M = ∫ ρ N ᵀN dV` by the element's own quadrature rule, via
+    /// [`ElementMatrixComputer`].
+    ///
+    /// # Degree-of-freedom layout
+    ///
+    /// Displacement components are interleaved per node — `[u₀, v₀, u₁, v₁,
+    /// …]` in 2-D — matching the strain-displacement matrix assembled by
+    /// [`ElementMatrixComputer::compute_stiffness_matrix`] and the node DOF
+    /// blocks handed out by `DofNumbering`.
+    ///
+    /// `compute_mass_matrix` integrates the *scalar* field mass,
+    /// `num_nodes × num_nodes`, because `∫ ρ Nᵢ Nⱼ dV` does not depend on the
+    /// component. The vector-valued mass matrix is its Kronecker product with
+    /// the spatial identity: each displacement component carries the same
+    /// mass and the components do not couple. Skipping that expansion is how
+    /// a mass matrix ends up dimensionally inconsistent with its stiffness.
+    ///
+    /// # Force vector
+    ///
+    /// Zero. Body forces and surface tractions enter through the boundary
+    /// condition set rather than here, so an element with no body load
+    /// contributes no force — which is different from the previous behaviour
+    /// of returning zero because nothing was computed at all.
     pub fn compute_element_matrices(
         &self,
-        _material: &dyn crate::materials::Material,
+        material: &dyn crate::materials::Material,
         _time: f64,
     ) -> FeaResult<ElementMatrix> {
+        let properties = material.properties();
         let num_nodes = self.num_nodes();
-        let dofs_per_node = self.spatial_dimension();
+        let dim = self.spatial_dimension();
-        let total_dofs = num_nodes * dofs_per_node;
+        let total_dofs = num_nodes * dim;
 
-        // Create placeholder matrices
+        let stiffness_matrix = ElementMatrixComputer::compute_stiffness_matrix(
-        let stiffness_matrix = nalgebra::DMatrix::zeros(total_dofs, total_dofs);
+            self,
-        let force_vector = nalgebra::DVector::zeros(total_dofs);
+            &self.node_coordinates,
-        let mass_matrix = Some(nalgebra::DMatrix::zeros(total_dofs, total_dofs));
+            properties.elastic_modulus,
+            properties.poisson_ratio,
+            None,
+        )?
+        .matrix;
+
+        let scalar_mass = ElementMatrixComputer::compute_consistent_mass_matrix(
+            self,
+            &self.node_coordinates,
+            properties.density,
+            None,
+        )?
+        .matrix;
+
+        // M_vector = M_scalar ⊗ I_dim, interleaved to match the DOF layout.
+        let mut mass_matrix = nalgebra::DMatrix::zeros(total_dofs, total_dofs);
+        for i in 0..num_nodes {
+            for j in 0..num_nodes {
+                let m_ij = scalar_mass[(i, j)];
+                for d in 0..dim {
+                    mass_matrix[(i * dim + d, j * dim + d)] = m_ij;
+                }
+            }
+        }
 
         Ok(ElementMatrix {
             stiffness_matrix,
-            force_vector,
+            force_vector: nalgebra::DVector::zeros(total_dofs),
-            mass_matrix,
+            mass_matrix: Some(mass_matrix),
         })
     }
 
@@ -411,14 +462,72 @@ impl FiniteElement for StandardFiniteElement {
         jacobian::JacobianEval::compute(self, coords, node_coords)
     }
 
+    /// Select the quadrature rule for this element's reference domain.
+    ///
+    /// This previously returned `QuadratureRule::new(vec![], ..)` — a rule
+    /// with **no points**. Every integration loop in the crate iterates over
+    /// `rule.points`, so an empty rule does not fail: it silently skips the
+    /// loop body and yields a zero matrix. That is why the element stiffness
+    /// and mass matrices were zero, and it is the reason an unsupported
+    /// element type below returns an error rather than an empty rule.
+    ///
+    /// The default order is chosen per family to integrate that element's
+    /// stiffness exactly: two Gauss points per direction for tri-/bi-linear
+    /// elements, three for the quadratic ones, and the corresponding
+    /// symmetric rules on simplices.
     fn quadrature_rule(&self, order: Option<usize>) -> FeaResult<QuadratureRule> {
-        let integration_order = order.unwrap_or(2);
+        use ElementType::{
-        // Create quadrature rule based on element type and integration order
+            Hex8, Hex20, Hex27, Line2, Line3, Point, Pyramid5, Pyramid13, Quad4, Quad8, Quad9,
-        Ok(quadrature::QuadratureRule::new(
+            Tet4, Tet10, Tri3, Tri6, Wedge6, Wedge15,
-            vec![],
+        };
-            integration_order,
+
-            self.element_type.spatial_dimension(),
+        let order = order.unwrap_or(match self.element_type {
-        ))
+            Point => 1,
+            Line2 | Tri3 | Tet4 => 1,
+            Line3 | Tri6 | Tet10 | Quad4 | Hex8 | Wedge6 => 2,
+            Quad8 | Quad9 | Hex20 | Hex27 | Wedge15 | Pyramid5 | Pyramid13 => 3,
+        });
+
+        match self.element_type {
+            Point => Ok(QuadratureRule::new(
+                vec![quadrature::QuadraturePoint::new_3d(0.0, 0.0, 0.0, 1.0)],
+                order,
+                0,
+            )),
+            Line2 | Line3 => QuadratureRule::line(order),
+            Tri3 | Tri6 => QuadratureRule::triangle(order),
+            Quad4 | Quad8 | Quad9 => QuadratureRule::quadrilateral(order),
+            Tet4 | Tet10 => QuadratureRule::tetrahedron(order),
+            Hex8 | Hex20 | Hex27 => QuadratureRule::hexahedron(order),
+
+            // A wedge is a triangle extruded along zeta, so its rule is the
+            // tensor product of the triangle rule with a 1-D Gauss rule.
+            Wedge6 | Wedge15 => {
+                let tri = QuadratureRule::triangle(order.min(3))?;
+                let line = QuadratureRule::line(order)?;
+                let mut points = Vec::with_capacity(tri.points.len() * line.points.len());
+                for t in &tri.points {
+                    for l in &line.points {
+                        points.push(quadrature::QuadraturePoint::new_3d(
+                            t.coords.xi(),
+                            t.coords.eta(),
+                            l.coords.xi(),
+                            t.weight * l.weight,
+                        ));
+                    }
+                }
+                Ok(QuadratureRule::new(points, order, 3))
+            }
+
+            // No collapsed-hexahedron rule is implemented for pyramids. An
+            // explicit error beats an empty rule that integrates to zero.
+            Pyramid5 | Pyramid13 => Err(ElementError::UnsupportedQuadratureOrder {
+                element_type: "pyramid".to_string(),
+                requested_order: order,
+                max_order: 0,
+            }
+            .into()),
+        }
     }
 }
 

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

@@ -653,6 +653,32 @@ impl Quadrilateral9 {
     pub fn new() -> Self {
         Self
     }
+
+    /// Which 1-D Lagrange basis pair `(i, j)` each node uses, with
+    /// `i` indexing `r ∈ {-1, 0, +1}` and `j` indexing `s ∈ {-1, 0, +1}`.
+    ///
+    /// The node order is the finite-element convention the rest of this crate
+    /// uses: four corners counter-clockwise, then four mid-edges starting at
+    /// the bottom, then the centre. `Quadrilateral8` and `Quadrilateral4`
+    /// both follow it, and a mesh built for one must be readable by the
+    /// others.
+    ///
+    /// Shape functions were previously emitted in raw lexicographic lattice
+    /// order — `(-1,-1), (0,-1), (1,-1), (-1,0), …` — so a mesh written in
+    /// the usual order paired each node with the wrong basis function. At the
+    /// element centre that made the Jacobian exactly singular, because the
+    /// scrambled ordering mapped the element onto a degenerate shape.
+    const NODE_LATTICE: [(usize, usize); 9] = [
+        (0, 0), // 0: (-1, -1)
+        (2, 0), // 1: (+1, -1)
+        (2, 2), // 2: (+1, +1)
+        (0, 2), // 3: (-1, +1)
+        (1, 0), // 4: ( 0, -1)
+        (2, 1), // 5: (+1,  0)
+        (1, 2), // 6: ( 0, +1)
+        (0, 1), // 7: (-1,  0)
+        (1, 1), // 8: ( 0,  0)
+    ];
 }
 
 impl ShapeFunctions for Quadrilateral9 {
@@ -682,13 +708,9 @@ impl ShapeFunctions for Quadrilateral9 {
             0.5 * s * (s + 1.0), // L₂(s)
         ];
 
-        // 2D tensor product
+        // 2D tensor product, emitted in finite-element node order.
-        let mut idx = 0;
+        for (node, &(i, j)) in Self::NODE_LATTICE.iter().enumerate() {
-        for j in 0..3 {
+            values[node] = l_r[i] * l_s[j];
-            for i in 0..3 {
-                values[idx] = l_r[i] * l_s[j];
-                idx += 1;
-            }
         }
 
         Ok(values)
@@ -715,14 +737,11 @@ impl ShapeFunctions for Quadrilateral9 {
 
         let dl_s = [s - 0.5, -2.0 * s, s + 0.5];
 
-        // Compute derivatives
+        // Same node order as `evaluate`, or the Jacobian would pair each
-        let mut idx = 0;
+        // shape function's derivative with a different node's coordinates.
-        for j in 0..3 {
+        for (node, &(i, j)) in Self::NODE_LATTICE.iter().enumerate() {
-            for i in 0..3 {
+            derivs[(node, 0)] = dl_r[i] * l_s[j]; // dN/dr
-                derivs[(idx, 0)] = dl_r[i] * l_s[j]; // dN/dr
+            derivs[(node, 1)] = l_r[i] * dl_s[j]; // dN/ds
-                derivs[(idx, 1)] = l_r[i] * dl_s[j]; // dN/ds
-                idx += 1;
-            }
         }
 
         Ok(derivs)

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

@@ -593,11 +593,15 @@ impl ShapeFunctions for Hexahedron20 {
         derivs[(0, 0)] = -0.125 * sm * tm * (-r - s - t - 2.0) - 0.125 * rm * sm * tm;
         derivs[(1, 0)] = 0.125 * sm * tm * (r - s - t - 2.0) + 0.125 * rp * sm * tm;
         derivs[(2, 0)] = 0.125 * sp * tm * (r + s - t - 2.0) + 0.125 * rp * sp * tm;
-        derivs[(3, 0)] = -0.125 * sp * tm * (-r + s - t - 2.0) + 0.125 * rm * sp * tm;
+        // N3 = 0.125 rm sp tm (-r + s - t - 2). Both the product rule term
+        // from d(rm)/dr = -1 and the one from d(-r + ...)/dr = -1 are
+        // negative; the second previously carried a `+`.
+        derivs[(3, 0)] = -0.125 * sp * tm * (-r + s - t - 2.0) - 0.125 * rm * sp * tm;
         derivs[(4, 0)] = -0.125 * sm * tp * (-r - s + t - 2.0) - 0.125 * rm * sm * tp;
         derivs[(5, 0)] = 0.125 * sm * tp * (r - s + t - 2.0) + 0.125 * rp * sm * tp;
         derivs[(6, 0)] = 0.125 * sp * tp * (r + s + t - 2.0) + 0.125 * rp * sp * tp;
-        derivs[(7, 0)] = -0.125 * sp * tp * (-r + s + t - 2.0) + 0.125 * rm * sp * tp;
+        // N7 = 0.125 rm sp tp (-r + s + t - 2), same correction as node 3.
+        derivs[(7, 0)] = -0.125 * sp * tp * (-r + s + t - 2.0) - 0.125 * rm * sp * tp;
 
         // dN/dr for mid-edge nodes
         derivs[(8, 0)] = -0.5 * r * sm * tm;
@@ -616,11 +620,15 @@ impl ShapeFunctions for Hexahedron20 {
         // dN/ds (similar pattern)
         // Corner nodes
         derivs[(0, 1)] = -0.125 * rm * tm * (-r - s - t - 2.0) - 0.125 * rm * sm * tm;
-        derivs[(1, 1)] = -0.125 * rp * tm * (r - s - t - 2.0) + 0.125 * rp * sm * tm;
+        // N1 = 0.125 rp sm tm (r - s - t - 2). d(sm)/ds and d(-s + ...)/ds are
+        // both -1, so both product-rule terms are negative; the second
+        // previously carried a `+`. Same correction as node 5 below, and the
+        // mirror of the one applied to nodes 3 and 7 in dN/dr.
+        derivs[(1, 1)] = -0.125 * rp * tm * (r - s - t - 2.0) - 0.125 * rp * sm * tm;
         derivs[(2, 1)] = 0.125 * rp * tm * (r + s - t - 2.0) + 0.125 * rp * sp * tm;
         derivs[(3, 1)] = 0.125 * rm * tm * (-r + s - t - 2.0) + 0.125 * rm * sp * tm;
         derivs[(4, 1)] = -0.125 * rm * tp * (-r - s + t - 2.0) - 0.125 * rm * sm * tp;
-        derivs[(5, 1)] = -0.125 * rp * tp * (r - s + t - 2.0) + 0.125 * rp * sm * tp;
+        derivs[(5, 1)] = -0.125 * rp * tp * (r - s + t - 2.0) - 0.125 * rp * sm * tp;
         derivs[(6, 1)] = 0.125 * rp * tp * (r + s + t - 2.0) + 0.125 * rp * sp * tp;
         derivs[(7, 1)] = 0.125 * rm * tp * (-r + s + t - 2.0) + 0.125 * rm * sp * tp;
 

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

@@ -275,11 +275,25 @@ impl ShapeFunctions for Wedge15 {
             values[6 + i] = tri_funcs[i] * n_plus;
         }
 
-        // Vertical mid-edges (corners only)
+        // Vertical mid-edges. `4 L_i n₋ n₊` is `L_i (1 - t²)`.
         values[12] = 4.0 * l1 * n_minus * n_plus; // Vertical edge at corner 1
         values[13] = 4.0 * l2 * n_minus * n_plus; // Vertical edge at corner 2
         values[14] = 4.0 * l3 * n_minus * n_plus; // Vertical edge at corner 3
 
+        // Corner correction, without which this element does not form a
+        // partition of unity.
+        //
+        // Introducing a mid-side node on a vertical edge adds `L_i (1 - t²)`
+        // to the sum, so the two corners sharing that edge must each give up
+        // half of it. Omitting the correction left the shape functions
+        // summing to `1 + 4 n₋ n₊`, which is 2 at mid-height — a quadratic
+        // wedge that doubled every field interpolated through it.
+        let bubble = 1.0 - t * t;
+        for (corner, l) in [l1, l2, l3].into_iter().enumerate() {
+            values[corner] -= 0.5 * l * bubble;
+            values[6 + corner] -= 0.5 * l * bubble;
+        }
+
         Ok(values)
     }
 
@@ -374,6 +388,24 @@ impl ShapeFunctions for Wedge15 {
         derivs[(13, 2)] = 4.0 * l2 * t_factor; // Vertical edge at corner 2
         derivs[(14, 2)] = 4.0 * l3 * t_factor; // Vertical edge at corner 3
 
+        // Derivative of the corner correction `-0.5 L_i (1 - t²)` applied in
+        // `evaluate`. The derivatives must carry it too, or they would not be
+        // the derivatives of the functions this element actually uses — and
+        // the sum of dN/dr, dN/ds and dN/dt would not vanish as a partition
+        // of unity requires.
+        let bubble = 1.0 - t * t;
+        // d(l1, l2, l3)/dr and /ds for the area coordinates.
+        let dl_dr = [-1.0, 1.0, 0.0];
+        let dl_ds = [-1.0, 0.0, 1.0];
+        for (corner, l) in [l1, l2, l3].into_iter().enumerate() {
+            for face in [0usize, 6] {
+                derivs[(face + corner, 0)] -= 0.5 * dl_dr[corner] * bubble;
+                derivs[(face + corner, 1)] -= 0.5 * dl_ds[corner] * bubble;
+                // d/dt of `-0.5 L (1 - t²)` is `+L t`.
+                derivs[(face + corner, 2)] += l * t;
+            }
+        }
+
         Ok(derivs)
     }
 
@@ -808,155 +840,37 @@ impl FiniteElement for Pyramid13 {
 }
 
 impl ShapeFunctions for Pyramid13 {
-    fn evaluate(&self, xi: &[f64]) -> FeaResult<DVector<f64>> {
+    /// Not implemented.
-        if xi.len() < 3 {
+    ///
-            return Err(crate::error::FeaError::InvalidInput(
+    /// The previous implementation was not a quadratic pyramid basis. Its
-                "Parametric coordinates must have at least 3 components".to_string(),
+    /// thirteen shape functions summed to 4 at the element centre rather than
-            ));
+    /// 1, so any field interpolated through it was scaled by a factor that
+    /// varied over the element, and `derivatives` allocated a 13x3 matrix
+    /// then wrote rows 13, 14 and 15 -- it had been copied from a sixteen-node
+    /// layout -- so it panicked with an out-of-bounds index before the wrong
+    /// values could be used.
+    ///
+    /// A correct 13-node pyramid basis is rational rather than polynomial,
+    /// and this crate has no pyramid quadrature rule to integrate it with
+    /// (see `StandardFiniteElement::quadrature_rule`), so implementing the
+    /// basis alone would not make the element usable. Both gaps are recorded
+    /// in `omni-cortex/docs/solver_status.md`. Returning an error names the
+    /// limitation at the point of use; `Pyramid5` is unaffected and works.
+    fn evaluate(&self, _xi: &[f64]) -> FeaResult<DVector<f64>> {
+        Err(crate::error::FeaError::InvalidInput(
+            "Pyramid13 shape functions are not implemented; use Pyramid5, or \
+             split the pyramid into tetrahedra"
+                .to_string(),
+        ))
     }
 
-        let r = xi[0];
+    /// Not implemented. See [`Pyramid13::evaluate`].
-        let s = xi[1];
+    fn derivatives(&self, _xi: &[f64]) -> FeaResult<DMatrix<f64>> {
-        let t = xi[2];
+        Err(crate::error::FeaError::InvalidInput(
-
+            "Pyramid13 shape function derivatives are not implemented; use \
-        // Pyramid13: quadratic base, linear sides
+             Pyramid5, or split the pyramid into tetrahedra"
-        // Nodes 0-3: base corners
+                .to_string(),
-        // Nodes 4-7: base mid-edges
+        ))
-        // Node 8: base center
-        // Nodes 9-12: vertical mid-edges
-        // Node 13: apex (not used in 13-node version, would be node 12)
-
-        let mut values = DVector::zeros(13);
-
-        // Special handling for degenerate case
-        let eps = 1e-10;
-        let w = if (1.0 - t).abs() < eps { eps } else { 1.0 - t };
-
-        // Modified coordinates for pyramid
-        let _xi_mod = r * w;
-        let _eta_mod = s * w;
-
-        // Base corner nodes
-        values[0] = 0.25 * (1.0 - r) * (1.0 - s) * (1.0 - t) * (-r - s - 2.0 * t + 1.0);
-        values[1] = 0.25 * (1.0 + r) * (1.0 - s) * (1.0 - t) * (r - s - 2.0 * t + 1.0);
-        values[2] = 0.25 * (1.0 + r) * (1.0 + s) * (1.0 - t) * (r + s - 2.0 * t + 1.0);
-        values[3] = 0.25 * (1.0 - r) * (1.0 + s) * (1.0 - t) * (-r + s - 2.0 * t + 1.0);
-
-        // Base mid-edge nodes
-        values[4] = 0.5 * (1.0 - r * r) * (1.0 - s) * (1.0 - t);
-        values[5] = 0.5 * (1.0 + r) * (1.0 - s * s) * (1.0 - t);
-        values[6] = 0.5 * (1.0 - r * r) * (1.0 + s) * (1.0 - t);
-        values[7] = 0.5 * (1.0 - r) * (1.0 - s * s) * (1.0 - t);
-
-        // Base center node
-        values[8] = (1.0 - r * r) * (1.0 - s * s) * (1.0 - t);
-
-        // Vertical mid-edge nodes
-        values[9] = 0.5 * (1.0 - r) * (1.0 - s) * t * (1.0 - t);
-        values[10] = 0.5 * (1.0 + r) * (1.0 - s) * t * (1.0 - t);
-        values[11] = 0.5 * (1.0 + r) * (1.0 + s) * t * (1.0 - t);
-        values[12] = 0.5 * (1.0 - r) * (1.0 + s) * t * (1.0 - t);
-
-        Ok(values)
-    }
-
-    fn derivatives(&self, xi: &[f64]) -> FeaResult<DMatrix<f64>> {
-        if xi.len() < 3 {
-            return Err(crate::error::FeaError::InvalidInput(
-                "Parametric coordinates must have at least 3 components".to_string(),
-            ));
-        }
-
-        let r = xi[0];
-        let s = xi[1];
-        let t = xi[2];
-
-        let mut derivs = DMatrix::zeros(13, 3);
-
-        // Simplified derivatives (full implementation would be more complex)
-        // dN/dr for base corner nodes
-        derivs[(0, 0)] = -0.25 * (1.0 - s) * (1.0 - t) * (-2.0 * r - s - 2.0 * t + 1.0);
-        derivs[(1, 0)] = 0.25 * (1.0 - s) * (1.0 - t) * (2.0 * r - s - 2.0 * t + 1.0);
-        derivs[(2, 0)] = 0.25 * (1.0 + s) * (1.0 - t) * (2.0 * r + s - 2.0 * t + 1.0);
-        derivs[(3, 0)] = -0.25 * (1.0 + s) * (1.0 - t) * (-2.0 * r + s - 2.0 * t + 1.0);
-
-        // dN/dr for base mid-edges
-        derivs[(4, 0)] = -r * (1.0 - s) * (1.0 - t);
-        derivs[(5, 0)] = 0.5 * (1.0 - s * s) * (1.0 - t);
-        derivs[(6, 0)] = -r * (1.0 + s) * (1.0 - t);
-        derivs[(7, 0)] = -0.5 * (1.0 - s * s) * (1.0 - t);
-
-        // dN/ds for base corner nodes
-        derivs[(0, 1)] = -0.25 * (1.0 - r) * (1.0 - t) * (-r - 2.0 * s - 2.0 * t + 1.0);
-        derivs[(1, 1)] = -0.25 * (1.0 + r) * (1.0 - t) * (r - 2.0 * s - 2.0 * t + 1.0);
-        derivs[(2, 1)] = 0.25 * (1.0 + r) * (1.0 - t) * (r + 2.0 * s - 2.0 * t + 1.0);
-        derivs[(3, 1)] = 0.25 * (1.0 - r) * (1.0 - t) * (-r + 2.0 * s - 2.0 * t + 1.0);
-
-        // dN/ds for base mid-edges
-        derivs[(4, 1)] = -0.5 * (1.0 - r * r) * (1.0 - t);
-        derivs[(5, 1)] = -s * (1.0 + r) * (1.0 - t);
-        derivs[(6, 1)] = 0.5 * (1.0 - r * r) * (1.0 - t);
-        derivs[(7, 1)] = -s * (1.0 - r) * (1.0 - t);
-
-        // dN/dt for base corner nodes
-        derivs[(0, 2)] = -0.25 * (1.0 - r) * (1.0 - s) * (-r - s - 2.0 * t + 1.0);
-        derivs[(1, 2)] = -0.25 * (1.0 + r) * (1.0 - s) * (r - s - 2.0 * t + 1.0);
-        derivs[(2, 2)] = -0.25 * (1.0 + r) * (1.0 + s) * (r + s - 2.0 * t + 1.0);
-        derivs[(3, 2)] = -0.25 * (1.0 - r) * (1.0 + s) * (-r + s - 2.0 * t + 1.0);
-
-        // dN/dt for base mid-edges
-        derivs[(4, 2)] = -0.5 * (1.0 - r * r) * (1.0 - s);
-        derivs[(5, 2)] = -0.5 * (1.0 - s * s) * (1.0 + r);
-        derivs[(6, 2)] = -0.5 * (1.0 - r * r) * (1.0 + s);
-        derivs[(7, 2)] = -0.5 * (1.0 - s * s) * (1.0 - r);
-
-        // Top layer nodes (8-15) - similar to base but with (1+t) factor
-        // dN/dr for top corner nodes
-        derivs[(8, 0)] = -0.25 * (1.0 - s) * (1.0 + t) * (-2.0 * r - s + 2.0 * t + 1.0);
-        derivs[(9, 0)] = 0.25 * (1.0 - s) * (1.0 + t) * (2.0 * r - s + 2.0 * t + 1.0);
-        derivs[(10, 0)] = 0.25 * (1.0 + s) * (1.0 + t) * (2.0 * r + s + 2.0 * t + 1.0);
-        derivs[(11, 0)] = -0.25 * (1.0 + s) * (1.0 + t) * (-2.0 * r + s + 2.0 * t + 1.0);
-
-        // dN/dr for top mid-edges
-        derivs[(12, 0)] = -r * (1.0 - s) * (1.0 + t);
-        derivs[(13, 0)] = 0.5 * (1.0 - s * s) * (1.0 + t);
-        derivs[(14, 0)] = -r * (1.0 + s) * (1.0 + t);
-        derivs[(15, 0)] = -0.5 * (1.0 - s * s) * (1.0 + t);
-
-        // dN/ds for top corner nodes
-        derivs[(8, 1)] = -0.25 * (1.0 - r) * (1.0 + t) * (-r - 2.0 * s + 2.0 * t + 1.0);
-        derivs[(9, 1)] = -0.25 * (1.0 + r) * (1.0 + t) * (r - 2.0 * s + 2.0 * t + 1.0);
-        derivs[(10, 1)] = 0.25 * (1.0 + r) * (1.0 + t) * (r + 2.0 * s + 2.0 * t + 1.0);
-        derivs[(11, 1)] = 0.25 * (1.0 - r) * (1.0 + t) * (-r + 2.0 * s + 2.0 * t + 1.0);
-
-        // dN/ds for top mid-edges
-        derivs[(12, 1)] = -0.5 * (1.0 - r * r) * (1.0 + t);
-        derivs[(13, 1)] = -s * (1.0 + r) * (1.0 + t);
-        derivs[(14, 1)] = 0.5 * (1.0 - r * r) * (1.0 + t);
-        derivs[(15, 1)] = -s * (1.0 - r) * (1.0 + t);
-
-        // dN/dt for top corner nodes
-        derivs[(8, 2)] = 0.25 * (1.0 - r) * (1.0 - s) * (-r - s + 2.0 * t + 1.0);
-        derivs[(9, 2)] = 0.25 * (1.0 + r) * (1.0 - s) * (r - s + 2.0 * t + 1.0);
-        derivs[(10, 2)] = 0.25 * (1.0 + r) * (1.0 + s) * (r + s + 2.0 * t + 1.0);
-        derivs[(11, 2)] = 0.25 * (1.0 - r) * (1.0 + s) * (-r + s + 2.0 * t + 1.0);
-
-        // dN/dt for top mid-edges
-        derivs[(12, 2)] = 0.5 * (1.0 - r * r) * (1.0 - s);
-        derivs[(13, 2)] = 0.5 * (1.0 - s * s) * (1.0 + r);
-        derivs[(14, 2)] = 0.5 * (1.0 - r * r) * (1.0 + s);
-        derivs[(15, 2)] = 0.5 * (1.0 - s * s) * (1.0 - r);
-
-        // Vertical mid-edges (16-19)
-        for i in 16..20 {
-            let _base_idx = i - 16; // 0,1,2,3
-            // These are the derivatives for vertical edge shape functions
-            derivs[(i, 0)] = 0.0; // No r-derivative for vertical edges
-            derivs[(i, 1)] = 0.0; // No s-derivative for vertical edges
-            derivs[(i, 2)] = 1.0; // Unit t-derivative for vertical edges
-        }
-
-        Ok(derivs)
     }
 
     fn num_nodes(&self) -> usize {
@@ -1024,12 +938,17 @@ mod tests {
 
     #[test]
     fn test_pyramid13() {
-        let elem = Pyramid13;
+        // Pyramid13 is not implemented and reports so. See the doc comment on
-        let xi = vec![0.0, 0.0, 0.0];
+        // its `evaluate`.
-        let shape = elem.evaluate(&xi).unwrap();
+        //
-
+        // This test previously accepted any positive sum, explaining that the
-        assert_eq!(shape.len(), 13);
+        // "special pyramid formulation" need not sum to 1. No formulation is
-        // Note: Sum may not be exactly 1.0 due to the special pyramid formulation
+        // exempt: shape functions that do not sum to unity cannot reproduce a
-        assert!(shape.sum() > 0.0);
+        // constant field, so the element could not represent a rigid-body
+        // translation. The actual sum was 4 at the element centre.
+        let error = Pyramid13
+            .evaluate(&[0.0, 0.0, 0.0])
+            .expect_err("Pyramid13 must report that it is unimplemented");
+        assert!(error.to_string().to_lowercase().contains("not implemented"));
     }
 }

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

@@ -287,10 +287,14 @@ pub struct MemoryInfo {
 
 impl std::fmt::Display for MemoryInfo {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        // Binary units. Driver APIs report memory in bytes and the sizes are
+        // powers of two, so 1 GiB = 1024^3; dividing by 1e9 instead labelled
+        // an 8 GiB device "8.59 GB".
+        const GIB: f64 = (1024 * 1024 * 1024) as f64;
         writeln!(f, "GPU Memory Information:")?;
-        writeln!(f, "  Total: {:.2} GB", self.total_memory as f64 / 1e9)?;
+        writeln!(f, "  Total: {:.2} GiB", self.total_memory as f64 / GIB)?;
-        writeln!(f, "  Free:  {:.2} GB", self.free_memory as f64 / 1e9)?;
+        writeln!(f, "  Free:  {:.2} GiB", self.free_memory as f64 / GIB)?;
-        writeln!(f, "  Used:  {:.2} GB", self.used_memory as f64 / 1e9)?;
+        writeln!(f, "  Used:  {:.2} GiB", self.used_memory as f64 / GIB)?;
         writeln!(
             f,
             "  Usage: {:.1}%",
@@ -585,11 +589,14 @@ mod tests {
             used_memory: 2 * 1024 * 1024 * 1024,  // 2 GB
         };
 
+        // The fixture builds binary gibibytes, so the display must report
+        // them as such. This previously expected "8.00 GB" from a formatter
+        // that divided by 1e9, which renders 8 GiB as 8.59.
         let display_str = format!("{}", info);
-        assert!(display_str.contains("Total: 8.00 GB"));
+        assert!(display_str.contains("Total: 8.00 GiB"), "{display_str}");
-        assert!(display_str.contains("Free:  6.00 GB"));
+        assert!(display_str.contains("Free:  6.00 GiB"), "{display_str}");
-        assert!(display_str.contains("Used:  2.00 GB"));
+        assert!(display_str.contains("Used:  2.00 GiB"), "{display_str}");
-        assert!(display_str.contains("Usage: 25.0%"));
+        assert!(display_str.contains("Usage: 25.0%"), "{display_str}");
     }
 
     #[test]

crates/specialized/rtx-fea/src/lib.rs

@@ -14,29 +14,68 @@
 //! - **Efficient Assembly**: Sparse matrix assembly with GPU optimization
 //! - **Robust Solvers**: Direct and iterative solvers with GPU acceleration
 //! - **Mesh Management**: Advanced mesh operations including refinement and partitioning
-//! - **Production Ready**: No mocks, stubs, or TODOs - complete implementation
+//!
+//! ## Maturity
+//!
+//! Validated: element stiffness and mass matrices, shape functions across the
+//! element library, global assembly, constraint handling, the generalized
+//! eigensolver, and modal analysis end to end against closed-form bar
+//! frequencies. See `tests/element_matrices_physical.rs`,
+//! `tests/shape_function_invariants.rs`, `tests/eigenvalue_closed_form.rs`
+//! and `tests/modal_closed_form.rs`.
+//!
+//! Not yet: `DynamicAnalysis` and `NonlinearAnalysis` return zeros;
+//! `Pyramid13` and pyramid quadrature are unimplemented and report so. Much
+//! of the crate's test suite is still disabled behind `#[cfg(disabled)]`, so
+//! a green run does not imply coverage of what those modules cover.
 //!
 //! ## Quick Start
 //!
-//! ```rust,no_run
+//! Natural frequencies of a fixed-free bar:
-//! use rtx_fea::{
+//!
-//!     mesh::{Mesh, geometry::Rectangle},
+//! ```rust
-//!     materials::LinearElastic,
+//! use rtx_fea::analysis::{Analysis, AnalysisConfig, AnalysisData, ModalAnalysis};
-//!     solvers::DirectSolver,
+//! use rtx_fea::assembly::DofComponent;
-//!     analysis::StaticAnalysis,
+//! use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, DirichletBC};
+//! use rtx_fea::materials::{LinearElastic, MaterialDatabase};
+//! use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node};
+//!
+//! // A 4 x 1 grid of quadrilaterals spanning a 1.0 x 0.05 strip.
+//! let mut mesh = Mesh::new(2)?;
+//! let mut columns = Vec::new();
+//! for i in 0..=4 {
+//!     let x = f64::from(i) * 0.25;
+//!     columns.push([
+//!         mesh.add_node(Node::new_2d(x, 0.0)),
+//!         mesh.add_node(Node::new_2d(x, 0.05)),
+//!     ]);
+//! }
+//! for i in 0..4 {
+//!     let nodes = vec![columns[i][0], columns[i + 1][0], columns[i + 1][1], columns[i][1]];
+//!     mesh.add_element(Element::new(ElementType::Quad4, nodes, MaterialId(0))?)?;
+//! }
+//!
+//! let mut materials = MaterialDatabase::new();
+//! materials.add_material(MaterialId(0), LinearElastic::new(200e9, 0.3).with_density(8000.0), None);
+//!
+//! // Clamp the left edge. Without constraints the stiffness matrix is
+//! // singular and the analysis reports that rather than returning noise.
+//! let mut boundary_conditions = BoundaryConditionSet::new();
+//! boundary_conditions.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
+//!     columns[0].to_vec(),
+//!     vec![DofComponent::DisplacementX, DofComponent::DisplacementY],
+//!     0.0,
+//! )));
+//!
+//! let mut analysis = ModalAnalysis::new(mesh, materials, 3, AnalysisConfig::default())
+//!     .with_boundary_conditions(boundary_conditions);
+//! let results = analysis.run()?;
+//!
+//! let AnalysisData::Vector(frequencies) = &results.additional_data["frequencies"] else {
+//!     panic!("frequencies should be a vector")
 //! };
-//!
+//! assert_eq!(frequencies.len(), 3);
-//! // Create a mesh
+//! assert!(frequencies.iter().all(|f| f.is_finite() && *f > 0.0));
-//! let rect = Rectangle::new(1.0, 1.0);
-//! let mesh = rect.generate_quad_mesh(10, 10, rtx_fea::mesh::MaterialId(0))?;
-//!
-//! // Set up analysis
-//! let material = LinearElastic::new(200e9, 0.3); // Steel
-//! let mut analysis = StaticAnalysis::new(mesh);
-//! analysis.add_material(material);
-//!
-//! // Apply boundary conditions and solve
-//! // ... (see examples for complete workflow)
 //! # Ok::<(), Box<dyn std::error::Error>>(())
 //! ```
 //!

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

@@ -24,6 +24,7 @@ use cudarc::driver::safe::CudaContext;
 use nalgebra::{DMatrix, Vector3, Vector6};
 use serde::{Deserialize, Serialize};
 use std::collections::HashMap;
+use std::sync::Arc;
 
 pub use composite::*;
 pub use damage::*;
@@ -333,7 +334,7 @@ pub trait Material: Send + Sync {
 /// Material database for managing multiple materials.
 pub struct MaterialDatabase {
     /// Materials storage
-    materials: HashMap<crate::mesh::MaterialId, Box<dyn Material>>,
+    materials: HashMap<crate::mesh::MaterialId, Arc<dyn Material>>,
     /// Material names
     names: HashMap<crate::mesh::MaterialId, String>,
 }
@@ -364,7 +365,7 @@ impl MaterialDatabase {
         name: Option<String>,
     ) {
         let name = name.unwrap_or_else(|| format!("Material_{}", id.as_usize()));
-        self.materials.insert(id, Box::new(material));
+        self.materials.insert(id, Arc::new(material));
         self.names.insert(id, name);
     }
 
@@ -403,16 +404,23 @@ impl MaterialDatabase {
 }
 
 impl Clone for MaterialDatabase {
+    /// Clone the database, sharing the materials themselves.
+    ///
+    /// Materials are immutable once registered, so the `Arc` handles can be
+    /// shared rather than deep-copied and no `clone_box` on the trait object
+    /// is needed.
+    ///
+    /// This previously cloned **names only** and silently dropped every
+    /// material, on the grounds that `Box<dyn Material>` cannot be cloned.
+    /// Because `GlobalAssembler` is constructed with `materials.clone()`,
+    /// that meant every assembler ever built received an empty database and
+    /// every analysis failed with `MaterialNotFound` — for a mesh that was
+    /// correctly specified.
     fn clone(&self) -> Self {
-        // For now, create an empty database with same names
+        Self {
-        // Material cloning would require dyn trait object cloning support
+            materials: self.materials.clone(),
-        let mut new_db = Self::new();
+            names: self.names.clone(),
-        for (id, name) in &self.names {
-            new_db.names.insert(*id, name.clone());
         }
-        // Note: materials are not cloned, this is a shallow clone for names only
-        // Full implementation would require Material trait to have clone_box method
-        new_db
     }
 }
 

crates/specialized/rtx-fea/src/solvers/eigenvalue.rs

@@ -102,6 +102,36 @@ impl EigenvalueSolver {
         let x = solve_lower(&l, &shifted)?;
         let b = solve_lower(&l, &x.transpose())?.transpose();
 
+        // Reject a rank-deficient problem explicitly.
+        //
+        // `try_inverse` does not fail on a matrix that is singular only to
+        // working precision, and an unconstrained structure is exactly that
+        // case: its rigid-body modes are zero eigenvalues that round to
+        // ~1e-16 of the largest. Inverting anyway succeeds and returns
+        // enormous, meaningless entries, which emerge as a handful of
+        // near-zero frequencies that look like real low-frequency modes.
+        // Better to say what is wrong and name the remedy.
+        //
+        // Only meaningful without a shift: `K - σM` is deliberately
+        // indefinite when a shift is used to bracket interior modes.
+        if self.shift.is_none() {
+            let spectrum = b.clone().symmetric_eigenvalues();
+            let largest = spectrum.iter().fold(0.0_f64, |acc, v| acc.max(v.abs()));
+            let smallest = spectrum
+                .iter()
+                .fold(f64::INFINITY, |acc: f64, v| acc.min(v.abs()));
+
+            if largest <= 0.0 || smallest / largest < 1e-12 {
+                return Err(SolverError::SolveError {
+                    reason: "stiffness matrix is singular — an unconstrained structure has \
+                             rigid-body modes at zero frequency. Constrain it, or pass an \
+                             explicit shift via `with_shift` to compute the modes around one."
+                        .to_string(),
+                }
+                .into());
+            }
+        }
+
         // A = B⁻¹. Symmetrized to shed the asymmetry that the inverse of a
         // numerically-symmetric matrix picks up at rounding level; Lanczos is
         // sensitive to it.

crates/specialized/rtx-fea/src/utils/math_tests.rs

@@ -12,7 +12,20 @@ mod tests {
     fn test_math_utils_constants() {
         assert!(MathUtils::EPSILON > 0.0);
         assert!(MathUtils::SMALL > 0.0);
-        assert!(MathUtils::SMALL < MathUtils::EPSILON * 1000.0);
+
+        // SMALL is the practical "is this zero" threshold and must sit well
+        // above machine epsilon, or comparisons against it would be decided
+        // by rounding noise.
+        //
+        // This previously asserted `SMALL < EPSILON * 1000`, which is false
+        // and inverts the relationship: EPSILON is f64::EPSILON (2.2e-16), so
+        // that bound is 2.2e-13, while SMALL is 1e-12. A threshold below a
+        // thousand machine epsilons would defeat its own purpose.
+        assert!(MathUtils::SMALL > MathUtils::EPSILON);
+
+        // ...and far below any length, stress or stiffness a model uses, so
+        // it never swallows a physically meaningful quantity.
+        assert!(MathUtils::SMALL < 1e-6);
     }
 
     #[test]
@@ -198,16 +211,31 @@ mod tests {
 
     #[test]
     fn test_von_mises_stress_biaxial() {
-        // Equal biaxial stress
+        // Equal biaxial stress: sigma_x = sigma_y = 100, sigma_z = 0.
+        //
+        // sigma_vm = sqrt(1/2 [(sx-sy)^2 + (sy-sz)^2 + (sz-sx)^2])
+        //          = sqrt(1/2 [0 + 100^2 + 100^2]) = 100.
+        //
+        // This previously expected 0, commented "no deviatoric stress". That
+        // is only true of a *hydrostatic* state, where all three principal
+        // stresses are equal. An equal biaxial state with a free third
+        // direction is strongly deviatoric -- and expecting zero would mean a
+        // biaxially loaded sheet could never yield.
         let stress = DVector::from_vec(vec![100.0, 100.0, 0.0, 0.0, 0.0, 0.0]);
         let vm_stress = VoigtUtils::von_mises_stress(&stress);
-        assert!(MathUtils::approx_eq(vm_stress, 0.0, 1e-10)); // No deviatoric stress
+        assert!(MathUtils::approx_eq(vm_stress, 100.0, 1e-10));
 
-        // Unequal biaxial stress
+        // Genuine hydrostatic stress, which does have zero von Mises stress.
+        // This is the case the comment above was reaching for.
+        let stress = DVector::from_vec(vec![100.0, 100.0, 100.0, 0.0, 0.0, 0.0]);
+        let vm_stress = VoigtUtils::von_mises_stress(&stress);
+        assert!(MathUtils::approx_eq(vm_stress, 0.0, 1e-10));
+
+        // Unequal biaxial stress: sqrt(1/2 [50^2 + 50^2 + 100^2]) = sqrt(7500).
+        // Not |100 - 50|; the von Mises stress is not a principal difference.
         let stress = DVector::from_vec(vec![100.0, 50.0, 0.0, 0.0, 0.0, 0.0]);
         let vm_stress = VoigtUtils::von_mises_stress(&stress);
-        let expected = 50.0; // |100 - 50| = 50
+        assert!(MathUtils::approx_eq(vm_stress, 7500.0_f64.sqrt(), 1e-10));
-        assert!(MathUtils::approx_eq(vm_stress, expected, 1e-10));
     }
 
     #[test]
@@ -281,13 +309,43 @@ mod tests {
     fn test_integrate_1d_trigonometric() {
         use std::f64::consts::PI;
 
-        // Test ∫ sin(x) dx from 0 to π = 2
+        // Gauss-Legendre with n points is exact for polynomials of degree
-        let result = IntegrationUtils::integrate_1d(|x| x.sin(), 0.0, PI, 3).unwrap();
+        // 2n-1, and sin is not a polynomial. Three points therefore give an
-        assert!(MathUtils::approx_eq(result, 2.0, 1e-10));
+        // approximation, not an exact answer -- the error here is ~1e-3.
+        //
+        // This previously demanded 1e-10 from a 3-point rule, which no
+        // correct implementation could satisfy. Rather than relax the bound
+        // to whatever happens to pass, assert what actually pins the
+        // quadrature: that the error *shrinks* as points are added. A wrong
+        // rule can land near the right answer once; it will not converge.
+        let exact = 2.0;
+        let errors: Vec<f64> = (1..=3)
+            .map(|n| {
+                let result = IntegrationUtils::integrate_1d(|x| x.sin(), 0.0, PI, n).unwrap();
+                (result - exact).abs()
+            })
+            .collect();
 
-        // Test ∫ cos(x) dx from 0 to π/2 = 1
+        for pair in errors.windows(2) {
+            assert!(
+                pair[1] < pair[0],
+                "adding a Gauss point did not reduce the error: {errors:?}"
+            );
+        }
+        assert!(
+            errors[2] < 2e-3,
+            "3-point Gauss on sin over [0, pi] should be accurate to ~1e-3, got {}",
+            errors[2]
+        );
+
+        // Same for cos over a quarter period.
         let result = IntegrationUtils::integrate_1d(|x| x.cos(), 0.0, PI / 2.0, 3).unwrap();
-        assert!(MathUtils::approx_eq(result, 1.0, 1e-10));
+        assert!(MathUtils::approx_eq(result, 1.0, 1e-3));
+
+        // A cubic *is* within the exactness of a 2-point rule, so there the
+        // tight tolerance is the right expectation.
+        let result = IntegrationUtils::integrate_1d(|x| x * x * x, 0.0, 2.0, 2).unwrap();
+        assert!(MathUtils::approx_eq(result, 4.0, 1e-12));
     }
 
     #[test]

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

@@ -0,0 +1,290 @@
+//! Physical invariants of the element stiffness and mass matrices.
+//!
+//! These are the cheapest decisive checks that element matrix computation is
+//! real. None of them needs a reference table: each is a property the exact
+//! matrices must satisfy for any correct implementation, so a failure
+//! localises to the quadrature, the Jacobian, the shape function derivatives
+//! or the constitutive matrix rather than to an accuracy budget.
+//!
+//! They exist because `StandardFiniteElement::compute_element_matrices`
+//! returned `DMatrix::zeros(..)` for every element in the mesh, which made
+//! every global matrix `GlobalAssembler` produced zero, and every analysis
+//! built on it a solve of a null system. A zero matrix passes a symmetry
+//! check and a "does not crash" check, so the invariants below are chosen to
+//! be ones a zero matrix fails.
+
+use approx::assert_relative_eq;
+use nalgebra::{DVector, Vector3};
+use rtx_fea::elements::StandardFiniteElement;
+use rtx_fea::materials::LinearElastic;
+use rtx_fea::mesh::ElementType;
+
+const E: f64 = 210e9;
+const NU: f64 = 0.3;
+const RHO: f64 = 7850.0;
+
+fn steel() -> LinearElastic {
+    LinearElastic::new(E, NU).with_density(RHO)
+}
+
+/// Unit square Quad4, counter-clockwise from the origin. Area = 1.
+fn unit_square() -> StandardFiniteElement {
+    StandardFiniteElement::new(
+        ElementType::Quad4,
+        vec![
+            Vector3::new(0.0, 0.0, 0.0),
+            Vector3::new(1.0, 0.0, 0.0),
+            Vector3::new(1.0, 1.0, 0.0),
+            Vector3::new(0.0, 1.0, 0.0),
+        ],
+    )
+}
+
+/// A 2 x 3 rectangle, to catch anything that only works when `det J == 1`.
+fn rectangle_2x3() -> StandardFiniteElement {
+    StandardFiniteElement::new(
+        ElementType::Quad4,
+        vec![
+            Vector3::new(0.0, 0.0, 0.0),
+            Vector3::new(2.0, 0.0, 0.0),
+            Vector3::new(2.0, 3.0, 0.0),
+            Vector3::new(0.0, 3.0, 0.0),
+        ],
+    )
+}
+
+/// The stub this suite exists to catch returned zeros everywhere.
+#[test]
+fn element_matrices_are_not_zero() {
+    let m = unit_square()
+        .compute_element_matrices(&steel(), 0.0)
+        .unwrap();
+
+    assert!(
+        m.stiffness_matrix.amax() > 0.0,
+        "stiffness matrix is entirely zero — element matrix computation is a stub"
+    );
+    let mass = m
+        .mass_matrix
+        .as_ref()
+        .expect("a mass matrix is required for modal and dynamic analysis");
+    assert!(
+        mass.amax() > 0.0,
+        "mass matrix is entirely zero — element matrix computation is a stub"
+    );
+}
+
+/// `K` must be symmetric: it is `∫ Bᵀ D B dV` with `D` symmetric.
+#[test]
+fn stiffness_is_symmetric() {
+    let m = rectangle_2x3()
+        .compute_element_matrices(&steel(), 0.0)
+        .unwrap();
+    let k = &m.stiffness_matrix;
+
+    for i in 0..k.nrows() {
+        for j in 0..k.ncols() {
+            assert_relative_eq!(k[(i, j)], k[(j, i)], epsilon = 1e-6 * k.amax());
+        }
+    }
+}
+
+/// A rigid-body translation stores no strain energy, so `K t = 0`.
+///
+/// This is the single most informative check on the strain-displacement
+/// matrix `B`: it fails for a sign error, a mis-transformed derivative, or a
+/// wrong DOF ordering, none of which a symmetry check detects. It is the
+/// element-level half of the patch test.
+#[test]
+fn rigid_translation_produces_no_internal_force() {
+    let element = rectangle_2x3();
+    let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
+    let k = &m.stiffness_matrix;
+
+    let num_nodes = 4;
+    let dim = 2;
+    assert_eq!(k.nrows(), num_nodes * dim);
+
+    // One translation per spatial direction.
+    for d in 0..dim {
+        let mut t = DVector::zeros(num_nodes * dim);
+        for node in 0..num_nodes {
+            t[node * dim + d] = 1.0;
+        }
+
+        let f = k * &t;
+        assert!(
+            f.amax() < 1e-6 * k.amax(),
+            "translating the element in direction {d} produced internal force {:.3e} \
+             against a stiffness scale of {:.3e}",
+            f.amax(),
+            k.amax()
+        );
+    }
+}
+
+/// An unconstrained plane element has exactly three rigid-body modes — two
+/// translations and one rotation — so `K` has a three-dimensional null space.
+///
+/// Fewer means the element is spuriously stiff; more means it is rank
+/// deficient and admits a zero-energy deformation (hourglassing), which shows
+/// up in a real analysis as a mode shape that is pure noise.
+#[test]
+fn stiffness_has_exactly_three_rigid_body_modes_in_2d() {
+    let m = rectangle_2x3()
+        .compute_element_matrices(&steel(), 0.0)
+        .unwrap();
+    let k = m.stiffness_matrix;
+
+    let eigenvalues = k.clone().symmetric_eigenvalues();
+    let scale = k.amax();
+    let num_zero = eigenvalues
+        .iter()
+        .filter(|&&e| e.abs() < 1e-9 * scale)
+        .count();
+
+    assert_eq!(
+        num_zero,
+        3,
+        "expected 3 rigid-body modes (2 translations + 1 rotation), found {num_zero}; \
+         eigenvalues (scaled): {:?}",
+        eigenvalues.iter().map(|e| e / scale).collect::<Vec<_>>()
+    );
+}
+
+/// The consistent mass matrix integrates to the element's total mass.
+///
+/// `Σᵢⱼ Mᵢⱼ = ∫ ρ (Σᵢ Nᵢ)(Σⱼ Nⱼ) dV = ∫ ρ dV = ρV`, using the partition of
+/// unity. This is decisive for the quadrature rule and the Jacobian
+/// determinant together: get either wrong and the total mass is wrong by
+/// exactly the factor of the error.
+#[test]
+fn consistent_mass_integrates_to_rho_times_volume() {
+    for (element, area) in [(unit_square(), 1.0), (rectangle_2x3(), 6.0)] {
+        let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
+        let mass = m.mass_matrix.expect("mass matrix required");
+
+        // Vector-valued mass: each spatial direction carries the full mass,
+        // so the whole matrix sums to `dim * rho * V`.
+        let dim = 2;
+        let total: f64 = mass.iter().sum();
+        assert_relative_eq!(total, dim as f64 * RHO * area, max_relative = 1e-9);
+    }
+}
+
+/// The mass matrix must be positive definite.
+///
+/// It is the metric in `K φ = λ M φ`; if it is not positive definite the
+/// eigenproblem has no real spectrum and the Cholesky reduction in the
+/// eigensolver fails outright.
+#[test]
+fn mass_matrix_is_positive_definite() {
+    let m = rectangle_2x3()
+        .compute_element_matrices(&steel(), 0.0)
+        .unwrap();
+    let mass = m.mass_matrix.expect("mass matrix required");
+
+    let eigenvalues = mass.clone().symmetric_eigenvalues();
+    let min = eigenvalues.iter().cloned().fold(f64::INFINITY, f64::min);
+    assert!(
+        min > 0.0,
+        "mass matrix has a non-positive eigenvalue {min:.3e}; \
+         eigenvalues: {eigenvalues:?}"
+    );
+}
+
+/// Mass scales linearly with density, stiffness with elastic modulus, and
+/// neither picks up the other's property.
+///
+/// A single hard-coded default leaking into the computation would break this,
+/// and would otherwise be invisible in a suite that only ever uses one
+/// material.
+#[test]
+fn matrices_scale_with_the_material_they_are_given() {
+    let base = steel();
+    let stiffer = LinearElastic::new(2.0 * E, NU).with_density(RHO);
+    let denser = LinearElastic::new(E, NU).with_density(3.0 * RHO);
+
+    let element = rectangle_2x3();
+    let m0 = element.compute_element_matrices(&base, 0.0).unwrap();
+    let m_stiff = element.compute_element_matrices(&stiffer, 0.0).unwrap();
+    let m_dense = element.compute_element_matrices(&denser, 0.0).unwrap();
+
+    assert_relative_eq!(
+        m_stiff.stiffness_matrix.amax(),
+        2.0 * m0.stiffness_matrix.amax(),
+        max_relative = 1e-12
+    );
+    assert_relative_eq!(
+        m_dense.mass_matrix.as_ref().unwrap().amax(),
+        3.0 * m0.mass_matrix.as_ref().unwrap().amax(),
+        max_relative = 1e-12
+    );
+
+    // Doubling E must not change the mass; tripling ρ must not change K.
+    assert_relative_eq!(
+        m_stiff.mass_matrix.as_ref().unwrap().amax(),
+        m0.mass_matrix.as_ref().unwrap().amax(),
+        max_relative = 1e-12
+    );
+    assert_relative_eq!(
+        m_dense.stiffness_matrix.amax(),
+        m0.stiffness_matrix.amax(),
+        max_relative = 1e-12
+    );
+}
+
+/// The same invariants must hold in 3D, where the null space is six
+/// dimensional — three translations and three rotations.
+#[test]
+fn hexahedron_satisfies_the_same_invariants() {
+    let element = StandardFiniteElement::new(
+        ElementType::Hex8,
+        vec![
+            Vector3::new(0.0, 0.0, 0.0),
+            Vector3::new(2.0, 0.0, 0.0),
+            Vector3::new(2.0, 1.0, 0.0),
+            Vector3::new(0.0, 1.0, 0.0),
+            Vector3::new(0.0, 0.0, 3.0),
+            Vector3::new(2.0, 0.0, 3.0),
+            Vector3::new(2.0, 1.0, 3.0),
+            Vector3::new(0.0, 1.0, 3.0),
+        ],
+    );
+    let volume = 2.0 * 1.0 * 3.0;
+
+    let m = element.compute_element_matrices(&steel(), 0.0).unwrap();
+    let k = m.stiffness_matrix;
+    let mass = m.mass_matrix.expect("mass matrix required");
+
+    assert_eq!(k.nrows(), 8 * 3);
+
+    // Rigid translations.
+    for d in 0..3 {
+        let mut t = DVector::zeros(24);
+        for node in 0..8 {
+            t[node * 3 + d] = 1.0;
+        }
+        let f = &k * &t;
+        assert!(
+            f.amax() < 1e-6 * k.amax(),
+            "3-D rigid translation in direction {d} produced force {:.3e}",
+            f.amax()
+        );
+    }
+
+    let scale = k.amax();
+    let num_zero = k
+        .clone()
+        .symmetric_eigenvalues()
+        .iter()
+        .filter(|&&e| e.abs() < 1e-9 * scale)
+        .count();
+    assert_eq!(
+        num_zero, 6,
+        "expected 6 rigid-body modes in 3-D, found {num_zero}"
+    );
+
+    let total: f64 = mass.iter().sum();
+    assert_relative_eq!(total, 3.0 * RHO * volume, max_relative = 1e-9);
+}

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

@@ -58,32 +58,52 @@ mod hexahedron20_tests {
         // RED: Test Jacobian computation
         let hex20 = Hexahedron20::new();
 
-        // Create a regular hex with side length 2
+        // The reference cube spanning [-1, 1] in each direction, in the node
-        let mut node_coords = Vec::new();
+        // order the shape functions are written in: four bottom corners
-
+        // counter-clockwise, four top corners, then the bottom, top and
-        // Corner nodes (8)
+        // vertical mid-edges.
-        for k in 0..2 {
+        //
-            for j in 0..2 {
+        // This previously listed the corners in lexicographic order and put
-                for i in 0..2 {
+        // *all twelve* mid-edge nodes at the origin, commented "simplified
-                    node_coords.push(Vector3::new(
+        // for test". That is not a degenerate hexahedron so much as not a
-                        2.0 * i as f64 - 1.0,
+        // hexahedron: collapsing the mid-side nodes to a point makes the
-                        2.0 * j as f64 - 1.0,
+        // mapping genuinely singular. It only passed because the dN/dr and
-                        2.0 * k as f64 - 1.0,
+        // dN/ds derivatives carried sign errors at nodes 1, 3, 5 and 7, which
-                    ));
+        // produced a non-zero determinant for a geometry that has none.
-                }
+        let node_coords = vec![
-            }
+            // Bottom corners
-        }
+            Vector3::new(-1.0, -1.0, -1.0),
-
+            Vector3::new(1.0, -1.0, -1.0),
-        // Mid-edge nodes (12) - simplified for test
+            Vector3::new(1.0, 1.0, -1.0),
-        for _ in 8..20 {
+            Vector3::new(-1.0, 1.0, -1.0),
-            node_coords.push(Vector3::new(0.0, 0.0, 0.0));
+            // Top corners
-        }
+            Vector3::new(-1.0, -1.0, 1.0),
+            Vector3::new(1.0, -1.0, 1.0),
+            Vector3::new(1.0, 1.0, 1.0),
+            Vector3::new(-1.0, 1.0, 1.0),
+            // Bottom mid-edges
+            Vector3::new(0.0, -1.0, -1.0),
+            Vector3::new(1.0, 0.0, -1.0),
+            Vector3::new(0.0, 1.0, -1.0),
+            Vector3::new(-1.0, 0.0, -1.0),
+            // Top mid-edges
+            Vector3::new(0.0, -1.0, 1.0),
+            Vector3::new(1.0, 0.0, 1.0),
+            Vector3::new(0.0, 1.0, 1.0),
+            Vector3::new(-1.0, 0.0, 1.0),
+            // Vertical mid-edges
+            Vector3::new(-1.0, -1.0, 0.0),
+            Vector3::new(1.0, -1.0, 0.0),
+            Vector3::new(1.0, 1.0, 0.0),
+            Vector3::new(-1.0, 1.0, 0.0),
+        ];
 
         let coords = NaturalCoords::new_3d(0.0, 0.0, 0.0);
         let jac = hex20.jacobian(&coords, &node_coords).unwrap();
 
-        // GREEN: For regular hex, Jacobian at center should be diagonal
+        // The element occupies its own reference domain, so the mapping is the
-        assert!(jac.determinant().abs() > 1e-10); // Non-singular
+        // identity and the Jacobian determinant is exactly 1.
+        assert!((jac.determinant() - 1.0).abs() < 1e-10);
     }
 
     #[test]

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

@@ -69,13 +69,29 @@ fn test_material_database() {
 
 #[test]
 fn test_element_factory() {
-    // Test element creation for all supported types
+    // `Point`, `Line2` and `Line3` have no interpolation over an area or
+    // volume, so the factory deliberately rejects them. This previously
+    // required every variant of `ElementType::all()` to construct, which
+    // could only pass if the factory stopped making that distinction.
+    let unsupported = [ElementType::Point, ElementType::Line2, ElementType::Line3];
+
     for element_type in ElementType::all() {
         let element = ElementFactory::create(element_type);
-        assert!(element.is_ok(), "Failed to create {:?}", element_type);
 
-        let elem = element.unwrap();
+        if unsupported.contains(&element_type) {
-        assert!(elem.num_nodes() > 0);
+            assert!(
+                element.is_err(),
+                "{element_type:?} has no area or volume interpolation and should be rejected"
+            );
+            continue;
+        }
+
+        let elem = element.unwrap_or_else(|e| panic!("failed to create {element_type:?}: {e}"));
+        assert_eq!(
+            elem.num_nodes(),
+            element_type.node_count(),
+            "{element_type:?} reported the wrong node count"
+        );
         assert!(elem.spatial_dimension() >= 1 && elem.spatial_dimension() <= 3);
     }
 }
@@ -87,15 +103,25 @@ fn test_dof_numbering() {
     let dof_numbering =
         AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::Sequential).unwrap();
 
-    assert_eq!(dof_numbering.total_dofs, 27); // 9 nodes × 3 displacement components
+    // 9 nodes x 2 displacement components: `generate_rectangle` builds a
-    assert_eq!(dof_numbering.num_free_dofs(), 27);
+    // planar mesh, which has no out-of-plane displacement to number.
+    //
+    // This previously expected 27, from numbering all three components on a
+    // 2-D mesh. That gave every node a degree of freedom the element
+    // matrices never supply, so assembly rejected every element with a
+    // dimension mismatch — and `comprehensive_tdd_tests::test_dof_numbering`
+    // asserts `num_nodes * 2` for the same situation, so the two tests
+    // contradicted each other.
+    assert_eq!(dof_numbering.total_dofs, 18);
+    assert_eq!(dof_numbering.num_free_dofs(), 18);
     assert_eq!(dof_numbering.num_constrained_dofs(), 0);
 
-    // Test bandwidth optimized numbering
+    // Test bandwidth optimized numbering. Reordering DOFs must not change how
+    // many there are.
     let optimized_numbering =
         AdvancedDofNumbering::displacement_only(&mesh, DofMappingStrategy::BandwidthOptimized)
             .unwrap();
-    assert_eq!(optimized_numbering.total_dofs, 27);
+    assert_eq!(optimized_numbering.total_dofs, 18);
 }
 
 #[test]

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

@@ -0,0 +1,331 @@
+//! Modal analysis against closed-form natural frequencies.
+//!
+//! This is the end-to-end check that mesh, DOF numbering, constraints,
+//! element matrices, global assembly and the eigensolver are all correct
+//! *together*. Each has its own unit tests; none of those would catch a
+//! mismatch between them, such as a mass matrix assembled in a different DOF
+//! order than its stiffness.
+//!
+//! # Why axial modes and not a cantilever
+//!
+//! The obvious benchmark is the bending frequency of a cantilever,
+//! `β₁L = 1.8751`. It is the wrong first test here. Bilinear `Quad4` elements
+//! suffer **shear locking** in bending: their assumed displacement field
+//! cannot represent pure bending without spurious shear strain, so a coarse
+//! mesh is far too stiff and reports frequencies well above the true value.
+//! A cantilever test would fail for a reason that has nothing to do with
+//! whether the code under test is correct, and tuning the tolerance until it
+//! passed would destroy its value as evidence.
+//!
+//! Longitudinal (axial) vibration has no such problem. The exact solution of
+//! the 1-D wave equation for a fixed-free bar is
+//!
+//! ```text
+//!     f_n = (2n - 1) / (4L) * sqrt(E / rho),    n = 1, 2, 3, ...
+//! ```
+//!
+//! and a plane-stress mesh with transverse motion suppressed reduces to
+//! exactly that problem. Linear elements with a consistent mass matrix
+//! converge to it from above at `O(h²)`, so a modest mesh lands within a
+//! fraction of a percent — tight enough that a real error cannot hide.
+//!
+//! Bending is still checked below, but as a *convergence* statement rather
+//! than a single tolerance, which is the honest way to assert on an element
+//! that is known to lock.
+
+use rtx_fea::analysis::{Analysis, AnalysisConfig, AnalysisData, ModalAnalysis};
+use rtx_fea::assembly::DofComponent;
+use rtx_fea::boundary::{BoundaryCondition, BoundaryConditionSet, DirichletBC};
+use rtx_fea::materials::{LinearElastic, MaterialDatabase};
+use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node, NodeId};
+
+const E: f64 = 200e9;
+const RHO: f64 = 8000.0;
+const LENGTH: f64 = 1.0;
+const HEIGHT: f64 = 0.05;
+
+/// A rectangular `nx` by `ny` grid of `Quad4` elements spanning
+/// `[0, LENGTH] x [0, HEIGHT]`, returned with its node grid so tests can pick
+/// out edges to constrain.
+fn bar_mesh(nx: usize, ny: usize) -> (Mesh, Vec<Vec<NodeId>>) {
+    let mut mesh = Mesh::new(2).unwrap();
+
+    let mut grid = vec![vec![NodeId(0); ny + 1]; nx + 1];
+    for (i, column) in grid.iter_mut().enumerate() {
+        for (j, slot) in column.iter_mut().enumerate() {
+            let x = LENGTH * i as f64 / nx as f64;
+            let y = HEIGHT * j as f64 / ny as f64;
+            *slot = mesh.add_node(Node::new_2d(x, y));
+        }
+    }
+
+    for i in 0..nx {
+        for j in 0..ny {
+            // Counter-clockwise, so the Jacobian determinant is positive.
+            let nodes = vec![
+                grid[i][j],
+                grid[i + 1][j],
+                grid[i + 1][j + 1],
+                grid[i][j + 1],
+            ];
+            let element = Element::new(ElementType::Quad4, nodes, MaterialId(0)).unwrap();
+            mesh.add_element(element).unwrap();
+        }
+    }
+
+    (mesh, grid)
+}
+
+/// Poisson's ratio is zero throughout.
+///
+/// This is a deliberate modelling choice, not a convenience: with `nu = 0`
+/// the axial and transverse responses decouple exactly, so the plane-stress
+/// model reduces to the 1-D bar the closed form describes. A non-zero
+/// Poisson's ratio would introduce a real physical difference between the two
+/// and the comparison would no longer be exact.
+fn steel_no_poisson() -> MaterialDatabase {
+    let mut materials = MaterialDatabase::new();
+    materials.add_material(
+        MaterialId(0),
+        LinearElastic::new(E, 0.0).with_density(RHO),
+        Some("steel".to_string()),
+    );
+    materials
+}
+
+fn frequencies_of(results: &rtx_fea::analysis::AnalysisResults) -> Vec<f64> {
+    match results
+        .additional_data
+        .get("frequencies")
+        .expect("frequencies missing")
+    {
+        AnalysisData::Vector(v) => v.iter().copied().collect(),
+        other => panic!("frequencies had unexpected type {other:?}"),
+    }
+}
+
+/// Longitudinal modes of a fixed-free bar against `f_n = (2n-1)/(4L)·√(E/ρ)`.
+#[test]
+fn axial_modes_match_the_closed_form_bar() {
+    let nx = 24;
+    let ny = 2;
+    let (mesh, grid) = bar_mesh(nx, ny);
+
+    let mut bcs = BoundaryConditionSet::new();
+
+    // Clamp the left edge axially.
+    let left_edge: Vec<NodeId> = grid[0].clone();
+    bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
+        left_edge,
+        vec![DofComponent::DisplacementX],
+        0.0,
+    )));
+
+    // Suppress transverse motion everywhere, reducing the plane-stress model
+    // to the 1-D bar the closed form describes. Without this the spectrum is
+    // interleaved with bending modes and the comparison is meaningless.
+    let all_nodes: Vec<NodeId> = grid.iter().flatten().copied().collect();
+    bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
+        all_nodes,
+        vec![DofComponent::DisplacementY],
+        0.0,
+    )));
+
+    let num_modes = 3;
+    let mut analysis = ModalAnalysis::new(
+        mesh,
+        steel_no_poisson(),
+        num_modes,
+        AnalysisConfig::default(),
+    )
+    .with_boundary_conditions(bcs);
+
+    let results = analysis.run().expect("modal analysis failed");
+    let computed = frequencies_of(&results);
+
+    let wave_speed = (E / RHO).sqrt();
+    for n in 1..=num_modes {
+        let exact = (2 * n - 1) as f64 / (4.0 * LENGTH) * wave_speed;
+        let got = computed[n - 1];
+        let relative_error = (got - exact).abs() / exact;
+
+        assert!(
+            relative_error < 0.01,
+            "mode {n}: computed {got:.4} Hz against exact {exact:.4} Hz \
+             ({:.3}% error)",
+            relative_error * 100.0
+        );
+
+        // Linear elements with a consistent mass matrix are stiffer than the
+        // continuum, so the discrete frequency must come in high. Landing
+        // below the exact value means something is wrong even if the
+        // magnitude looks plausible.
+        assert!(
+            got >= exact * (1.0 - 1e-9),
+            "mode {n}: computed {got:.4} Hz is below the exact {exact:.4} Hz; \
+             a consistent-mass discretisation cannot be softer than the continuum"
+        );
+    }
+}
+
+/// Refining the mesh must drive the axial error down, and at the expected
+/// second-order rate.
+///
+/// A single tolerance check can be satisfied by a wrong formula with a
+/// compensating error. A convergence *rate* cannot: it pins the
+/// discretisation itself.
+#[test]
+fn axial_frequency_converges_at_second_order() {
+    let wave_speed = (E / RHO).sqrt();
+    let exact = wave_speed / (4.0 * LENGTH);
+
+    let mut errors = Vec::new();
+    for nx in [4usize, 8, 16] {
+        let (mesh, grid) = bar_mesh(nx, 1);
+
+        let mut bcs = BoundaryConditionSet::new();
+        bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
+            grid[0].clone(),
+            vec![DofComponent::DisplacementX],
+            0.0,
+        )));
+        bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
+            grid.iter().flatten().copied().collect(),
+            vec![DofComponent::DisplacementY],
+            0.0,
+        )));
+
+        let mut analysis =
+            ModalAnalysis::new(mesh, steel_no_poisson(), 1, AnalysisConfig::default())
+                .with_boundary_conditions(bcs);
+        let results = analysis.run().expect("modal analysis failed");
+        errors.push((frequencies_of(&results)[0] - exact).abs() / exact);
+    }
+
+    for window in errors.windows(2) {
+        let rate = (window[0] / window[1]).log2();
+        assert!(
+            rate > 1.7,
+            "halving the element size reduced the error by only 2^{rate:.2}; \
+             expected close to second order. errors: {errors:?}"
+        );
+    }
+}
+
+/// An unconstrained structure has rigid-body modes, so `K` is singular.
+///
+/// The failure must be a clear error rather than a set of near-zero
+/// eigenvalues that look like real low-frequency modes.
+#[test]
+fn unconstrained_structure_is_rejected_rather_than_silently_wrong() {
+    let (mesh, _) = bar_mesh(4, 1);
+    let mut analysis = ModalAnalysis::new(mesh, steel_no_poisson(), 2, AnalysisConfig::default());
+
+    let error = analysis
+        .run()
+        .expect_err("an unconstrained structure must not yield frequencies");
+    let message = error.to_string().to_lowercase();
+    assert!(
+        message.contains("singular") || message.contains("shift"),
+        "error should name the singular stiffness or the shift remedy, got: {error}"
+    );
+}
+
+/// Every reported natural frequency must be real and positive.
+///
+/// A constrained, positive-definite structure has no zero-frequency mode. A
+/// zero or NaN here means the constraints did not reach the assembled system
+/// or the eigenvalues came back negative.
+#[test]
+fn frequencies_are_real_and_positive() {
+    let (mesh, grid) = bar_mesh(6, 2);
+
+    let mut bcs = BoundaryConditionSet::new();
+    bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
+        grid[0].clone(),
+        vec![DofComponent::DisplacementX, DofComponent::DisplacementY],
+        0.0,
+    )));
+
+    let mut analysis = ModalAnalysis::new(mesh, steel_no_poisson(), 4, AnalysisConfig::default())
+        .with_boundary_conditions(bcs);
+    let results = analysis.run().expect("modal analysis failed");
+
+    let frequencies = frequencies_of(&results);
+    assert_eq!(frequencies.len(), 4);
+    for (i, f) in frequencies.iter().enumerate() {
+        assert!(
+            f.is_finite() && *f > 0.0,
+            "mode {} frequency is {f}, which is not a physical frequency",
+            i + 1
+        );
+    }
+
+    // Ascending, since modes are named by index.
+    for pair in frequencies.windows(2) {
+        assert!(
+            pair[1] >= pair[0],
+            "frequencies are not ascending: {frequencies:?}"
+        );
+    }
+}
+
+/// Cantilever bending, asserted as convergence rather than as a tolerance.
+///
+/// `Quad4` locks in bending, so the coarse-mesh frequency is far too high.
+/// What must still hold is that refinement moves it monotonically *towards*
+/// the Euler-Bernoulli value `f₁ = (β₁L)²/(2πL²)·√(EI/ρA)` with
+/// `β₁L = 1.8751` — and that it approaches from above, which is the signature
+/// of locking rather than of a bug.
+#[test]
+fn cantilever_bending_converges_towards_euler_bernoulli_from_above() {
+    let beta_l: f64 = 1.8751;
+    // Plane stress with unit thickness: A = h, I = h³/12.
+    let area = HEIGHT;
+    let second_moment = HEIGHT.powi(3) / 12.0;
+    let exact = beta_l.powi(2) / (2.0 * std::f64::consts::PI * LENGTH.powi(2))
+        * (E * second_moment / (RHO * area)).sqrt();
+
+    let mut computed = Vec::new();
+    for (nx, ny) in [(8usize, 2usize), (16, 4), (32, 8)] {
+        let (mesh, grid) = bar_mesh(nx, ny);
+
+        let mut bcs = BoundaryConditionSet::new();
+        bcs.add_condition(BoundaryCondition::Dirichlet(DirichletBC::fixed(
+            grid[0].clone(),
+            vec![DofComponent::DisplacementX, DofComponent::DisplacementY],
+            0.0,
+        )));
+
+        let mut analysis =
+            ModalAnalysis::new(mesh, steel_no_poisson(), 1, AnalysisConfig::default())
+                .with_boundary_conditions(bcs);
+        let results = analysis.run().expect("modal analysis failed");
+        computed.push(frequencies_of(&results)[0]);
+    }
+
+    for (i, f) in computed.iter().enumerate() {
+        assert!(
+            *f > exact * 0.95,
+            "mesh {i}: {f:.3} Hz is below the Euler-Bernoulli value {exact:.3} Hz; \
+             a locking element cannot be softer than the beam theory it approximates"
+        );
+    }
+
+    for pair in computed.windows(2) {
+        assert!(
+            pair[1] <= pair[0] * 1.001,
+            "refining the mesh increased the bending frequency ({:?}); \
+             locking must relax under refinement, not worsen",
+            computed
+        );
+    }
+
+    let coarse_error = (computed[0] - exact).abs() / exact;
+    let fine_error = (computed[computed.len() - 1] - exact).abs() / exact;
+    assert!(
+        fine_error < coarse_error,
+        "refinement did not reduce the bending error: {coarse_error:.4} -> {fine_error:.4} \
+         against exact {exact:.3} Hz, computed {computed:?}"
+    );
+}

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

@@ -100,21 +100,23 @@ mod pyramid13_tests {
 
     #[test]
     fn test_pyramid13_partition_of_unity() {
-        // RED: Test that shape functions sum to 1 everywhere
+        // Pyramid13 is not implemented; it reports that rather than returning
+        // a basis that does not form a partition of unity.
+        //
+        // The previous implementation summed to 4 at the base centre, and its
+        // `derivatives` allocated a 13x3 matrix then wrote rows 13 to 15 --
+        // copied from a sixteen-node layout -- so this test panicked on an
+        // out-of-bounds index rather than on the assertion below.
+        //
+        // Pyramid5 is unaffected; see `test_pyramid5_partition_of_unity`.
         let pyramid13 = Pyramid13::new();
 
-        let test_points = vec![
+        for (x, y, z) in [(0.0, 0.0, 0.0), (0.25, 0.25, 0.25)] {
-            (0.0, 0.0, 0.0),    // Base center
-            (0.25, 0.25, 0.25), // Quarter pyramid
-        ];
-
-        for (x, y, z) in test_points {
             let coords = NaturalCoords::new_3d(x, y, z);
-            let shape = pyramid13.shape_functions(&coords).unwrap();
+            let error = pyramid13
-
+                .shape_functions(&coords)
-            // GREEN: Sum of all shape functions should be 1
+                .expect_err("Pyramid13 must report that it is unimplemented");
-            let sum: f64 = (0..13).map(|i| shape.value(i).unwrap()).sum();
+            assert!(error.to_string().to_lowercase().contains("not implemented"));
-            assert!((sum - 1.0).abs() < 1e-9); // Slightly more tolerance for complex element
         }
     }
 }

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

@@ -0,0 +1,157 @@
+//! Invariants every finite element shape function set must satisfy.
+//!
+//! These hold for any correct element regardless of order or geometry, so a
+//! failure localises to the basis itself rather than to an accuracy budget.
+//! Checking them across the whole element library at once is what catches an
+//! element that was written but never exercised — `Wedge15` summed to 2 at
+//! mid-height and `Quadrilateral9` emitted its functions in a node order no
+//! other element in this crate uses.
+//!
+//! Two invariants, and the second is the one that gets skipped:
+//!
+//! 1. **Partition of unity**: `Σ Nᵢ(ξ) = 1` everywhere. An interpolation that
+//!    does not reproduce a constant cannot represent rigid-body motion.
+//! 2. **Vanishing derivative sum**: `Σ ∂Nᵢ/∂ξⱼ = 0` everywhere. This follows
+//!    from the first by differentiation, but the derivatives are usually
+//!    written out by hand separately from the values, so they can drift apart
+//!    — and when they do, the shape functions are consistent while the strain
+//!    they produce is not.
+
+use rtx_fea::elements::shape_functions::{
+    Hexahedron8, Hexahedron20, Pyramid5, Pyramid13, Quadrilateral4, Quadrilateral8, Quadrilateral9,
+    ShapeFunctions, Tetrahedron4, Tetrahedron10, Triangle3, Triangle6, Wedge6, Wedge15,
+};
+
+/// Sample points inside the reference domain of a quadrilateral or hexahedron,
+/// which spans `[-1, 1]` in each direction.
+const CUBE_SAMPLES: &[[f64; 3]] = &[
+    [0.0, 0.0, 0.0],
+    [0.5, -0.25, 0.75],
+    [-0.9, 0.9, -0.3],
+    [0.33, 0.67, 0.1],
+    [1.0, -1.0, 1.0],
+];
+
+/// Sample points inside a simplex reference domain, where the coordinates are
+/// non-negative and sum to at most one.
+const SIMPLEX_SAMPLES: &[[f64; 3]] = &[
+    [0.25, 0.25, 0.25],
+    [0.1, 0.2, 0.3],
+    [0.0, 0.0, 0.0],
+    [1.0 / 3.0, 1.0 / 3.0, 0.0],
+    [0.5, 0.25, 0.125],
+];
+
+/// Wedges are a triangle in `(r, s)` extruded over `t ∈ [-1, 1]`, so they need
+/// their own sample set. `t = 0` is included deliberately: it is where a
+/// missing vertical mid-edge correction shows up most strongly.
+const WEDGE_SAMPLES: &[[f64; 3]] = &[
+    [1.0 / 3.0, 1.0 / 3.0, 0.0],
+    [0.25, 0.25, 0.0],
+    [0.2, 0.3, -0.6],
+    [0.5, 0.1, 0.8],
+    [0.0, 0.0, -1.0],
+];
+
+fn check_element<S: ShapeFunctions>(name: &str, element: &S, samples: &[[f64; 3]]) {
+    let expected_nodes = element.num_nodes();
+
+    for point in samples {
+        let xi = &point[..];
+
+        let values = element
+            .evaluate(xi)
+            .unwrap_or_else(|e| panic!("{name}: evaluate failed at {point:?}: {e}"));
+        assert_eq!(
+            values.len(),
+            expected_nodes,
+            "{name}: returned {} values for {expected_nodes} nodes",
+            values.len()
+        );
+
+        let sum: f64 = values.iter().sum();
+        assert!(
+            (sum - 1.0).abs() < 1e-10,
+            "{name}: shape functions sum to {sum} at {point:?}, not 1 — \
+             the element cannot reproduce a constant field"
+        );
+
+        let derivatives = element
+            .derivatives(xi)
+            .unwrap_or_else(|e| panic!("{name}: derivatives failed at {point:?}: {e}"));
+        assert_eq!(
+            derivatives.nrows(),
+            expected_nodes,
+            "{name}: derivative matrix has {} rows for {expected_nodes} nodes",
+            derivatives.nrows()
+        );
+
+        for direction in 0..derivatives.ncols() {
+            let column_sum: f64 = derivatives.column(direction).iter().sum();
+            assert!(
+                column_sum.abs() < 1e-10,
+                "{name}: d/d(xi_{direction}) of the shape functions sums to \
+                 {column_sum} at {point:?}, not 0 — the derivatives are not \
+                 those of the values"
+            );
+        }
+    }
+}
+
+#[test]
+fn triangles_satisfy_the_invariants() {
+    check_element("Triangle3", &Triangle3, SIMPLEX_SAMPLES);
+    check_element("Triangle6", &Triangle6, SIMPLEX_SAMPLES);
+}
+
+#[test]
+fn quadrilaterals_satisfy_the_invariants() {
+    check_element("Quadrilateral4", &Quadrilateral4, CUBE_SAMPLES);
+    check_element("Quadrilateral8", &Quadrilateral8, CUBE_SAMPLES);
+    check_element("Quadrilateral9", &Quadrilateral9, CUBE_SAMPLES);
+}
+
+#[test]
+fn tetrahedra_satisfy_the_invariants() {
+    check_element("Tetrahedron4", &Tetrahedron4, SIMPLEX_SAMPLES);
+    check_element("Tetrahedron10", &Tetrahedron10, SIMPLEX_SAMPLES);
+}
+
+#[test]
+fn hexahedra_satisfy_the_invariants() {
+    check_element("Hexahedron8", &Hexahedron8, CUBE_SAMPLES);
+    check_element("Hexahedron20", &Hexahedron20, CUBE_SAMPLES);
+}
+
+#[test]
+fn wedges_satisfy_the_invariants() {
+    check_element("Wedge6", &Wedge6, WEDGE_SAMPLES);
+    check_element("Wedge15", &Wedge15, WEDGE_SAMPLES);
+}
+
+#[test]
+fn pyramid5_satisfies_the_invariants() {
+    check_element("Pyramid5", &Pyramid5, SIMPLEX_SAMPLES);
+}
+
+/// `Pyramid13` is not implemented and must say so rather than return values.
+///
+/// It previously returned a basis summing to 4 at the element centre, and its
+/// derivative routine indexed past the end of the matrix it had allocated.
+/// Reporting the gap is the honest behaviour until a correct rational basis
+/// and a matching pyramid quadrature rule both exist.
+#[test]
+fn pyramid13_reports_that_it_is_unimplemented() {
+    let error = Pyramid13
+        .evaluate(&[0.0, 0.0, 0.0])
+        .expect_err("Pyramid13 must not return shape function values");
+    assert!(
+        error.to_string().to_lowercase().contains("not implemented"),
+        "error should say the element is unimplemented, got: {error}"
+    );
+
+    let error = Pyramid13
+        .derivatives(&[0.0, 0.0, 0.0])
+        .expect_err("Pyramid13 must not return shape function derivatives");
+    assert!(error.to_string().to_lowercase().contains("not implemented"));
+}

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

@@ -31,21 +31,45 @@ mod standalone_tests {
         );
     }
 
+    /// Exercise the mesh operations rather than grepping for their names.
+    ///
+    /// This previously searched the text of `src/mesh/mod.rs` for the strings
+    /// "add_node", "add_element" and "generate_rectangle". It broke when the
+    /// implementations moved into submodules, which is the smaller problem:
+    /// the larger one is that a source-text search cannot distinguish a
+    /// working function from one that returns zeros. Every such check in this
+    /// file passed for the entire period during which element matrix
+    /// computation was a stub returning `DMatrix::zeros`, quadrature returned
+    /// no points at all, and cloning the material database silently dropped
+    /// every material.
     #[test]
     fn test_mesh_has_real_algorithms() {
-        // Verify mesh module has real implementations
+        use rtx_fea::mesh::{Element, ElementType, MaterialId, Mesh, Node};
-        let mesh_path = include_str!("../src/mesh/mod.rs");
 
-        // Check for real mesh operations
+        let mut mesh = Mesh::new(2).unwrap();
-        assert!(mesh_path.contains("add_node"), "Missing node addition");
+
-        assert!(
+        let n0 = mesh.add_node(Node::new_2d(0.0, 0.0));
-            mesh_path.contains("add_element"),
+        let n1 = mesh.add_node(Node::new_2d(1.0, 0.0));
-            "Missing element addition"
+        let n2 = mesh.add_node(Node::new_2d(1.0, 1.0));
-        );
+        let n3 = mesh.add_node(Node::new_2d(0.0, 1.0));
-        assert!(
+        assert_eq!(mesh.num_nodes(), 4);
-            mesh_path.contains("generate_rectangle"),
+
-            "Missing mesh generation"
+        mesh.add_element(
-        );
+            Element::new(ElementType::Quad4, vec![n0, n1, n2, n3], MaterialId(0)).unwrap(),
+        )
+        .unwrap();
+        assert_eq!(mesh.num_elements(), 1);
+
+        // Nodes come back at the coordinates they went in at.
+        let position = mesh.get_node(n2).expect("node 2 should exist").position();
+        assert!((position.x - 1.0).abs() < 1e-12);
+        assert!((position.y - 1.0).abs() < 1e-12);
+
+        // `generate_rectangle` takes node counts per direction, so a 3 x 3
+        // grid of nodes yields 9 nodes and 2 x 2 = 4 quadrilaterals.
+        let generated = Mesh::generate_rectangle(1.0, 1.0, 3, 3).unwrap();
+        assert_eq!(generated.num_nodes(), 9);
+        assert_eq!(generated.num_elements(), 4);
     }
 
     #[test]