rtx-fea: re-enable the remaining CPU test modules; fix three real defects they caught

All 33 remaining #[cfg(disabled)] test modules outside the GPU cluster are
now enabled: assembly (dof_mapping, constraints, global assembly), boundary
(mod + dirichlet/neumann/robin/thermal/contact), analysis (mod + static),
materials (mod, linear_elastic, hyperelastic, plasticity), elements (mod,
element_matrices, isoparametric, jacobian, quadrature), mesh (element_types,
connectivity, topology, topology_repair), solvers (mod, direct, iterative,
nonlinear) and lib.rs. Lib tests 117 -> 335, stable across repeated runs.
Only gpu_solver_tests and the GpuMeshData fixture stay disabled — they need
CUDA hardware and belong to the GPU tranche.

Three real defects found by the newly-compiling tests, each fixed:

- Direct solvers reused factorizations keyed on matrix SIZE alone.
  In a Newton loop the Jacobian changes every iteration but never its
  dimension, so LuDirect/CholeskyDirect/LdltDirect silently solved with the
  first iteration's factorization forever — Newton on x^2-4 crawled to
  x=1.955 in 1000 iterations instead of converging in 5. Invisible in
  single-solve linear analysis, which is why every green test passed over
  it. solve() now factorizes the matrix it is given.

- AdaptiveQuadrature's refinement re-integrated the WHOLE domain once per
  subdomain, so each level multiplied the estimate by the subdomain count:
  integrating e^x over [-1,1] at tolerance 1e-10 returned ~75 instead of
  2.35. The recursion now descends into each sub-box with its share of the
  error budget.

- compute_skewness read Jacobian columns as coordinate-line tangents, but
  the trait's jacobian() stores tangents in ROWS: on a sheared
  parallelogram whose tangents meet at 14 degrees it reported skewness 0.43
  instead of 0.84 — measuring per-component gradients, not mesh skew.

Fixtures corrected rather than the code where the fixture was wrong:
sigma_yy ~ 0 asserted uniaxial-stress physics on a uniaxial-strain state
(exact Lame values now asserted); an "unstable" orthotropic parameter set
that satisfies the determinant stability condition (delta = 0.187 > 0); a
unit-cube hex Jacobian of 1.0 that assumed a unit reference element (it is
0.125 from [-1,1]^3); a "distorted" quad whose centre Jacobian is exactly
orthogonal, asserted as skewed (flattening and shearing now tested
separately); a quality score below the implementation's own calibration;
Rayleigh damping fed the scalar-field mass (now expanded via the Kronecker
identity, with C = alpha*M + beta*K asserted entry-wise); an element
factory required to construct Point/Line types that have no implementation;
and DOF counts that encoded the repaired 3-DOFs-per-node-on-2-D defect.

MaterialDatabase::add_material call sites updated to the (id, material,
name) signature; ConnectivityInfo::build takes elements only;
TopologyRepair::triangle_quality (normalized 4*sqrt(3)*A/sum(a^2)) added
for the repair tests; create_subdomain_rule_* widened to pub(super) for the
quadrature tests.

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

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

@@ -12,7 +12,7 @@ pub mod quadrature_types;
 pub use quadrature_adaptive::*;
 pub use quadrature_types::*;
 
-#[cfg(disabled)]
+#[cfg(test)]
 mod tests {
     use super::*;
 
@@ -247,9 +247,19 @@ mod tests {
             .integrate(|coords| coords.xi().exp())
             .unwrap();
 
-        // Tight tolerance should be more accurate
-        let expected = (1.0_f64.exp() - (-1.0_f64).exp()); // e - e^(-1)
-        assert!((result_tight - expected).abs() < (result_loose - expected).abs());
+        // Both must actually approximate the integral — the broken recursion
+        // this test caught returned 2^levels times the estimate — and the
+        // tight tolerance must not be less accurate than the loose one.
+        let expected = 1.0_f64.exp() - (-1.0_f64).exp(); // e - e^(-1)
+        assert!(
+            (result_tight - expected).abs() < 1e-7,
+            "tight: expected {expected}, got {result_tight}"
+        );
+        assert!(
+            (result_loose - expected).abs() < 1e-2,
+            "loose: expected {expected}, got {result_loose}"
+        );
+        assert!((result_tight - expected).abs() <= (result_loose - expected).abs());
     }
 
     #[test]

crates/specialized/rtx-fea/src/elements/quadrature/quadrature_adaptive.rs

@@ -29,111 +29,97 @@ impl AdaptiveQuadrature {
     }
 
     /// Integrate with adaptive refinement.
+    ///
+    /// Bisects every axis of any sub-box whose one-level Richardson error
+    /// estimate exceeds its share of the tolerance, and recurses **into that
+    /// sub-box**. (The previous implementation re-integrated the whole
+    /// original domain once per subdomain, so each refinement level
+    /// *multiplied* the estimate by the subdomain count instead of improving
+    /// it.)
     pub fn integrate<F>(&self, function: F) -> FeaResult<f64>
     where
         F: Fn(&NaturalCoords) -> f64,
     {
-        self.integrate_recursive(&function, 0)
+        let dim = self.base_rule.dimension;
+        if !(1..=3).contains(&dim) {
+            return Err(ElementError::InvalidQuadrature {
+                reason: format!("Unsupported quadrature dimension: {dim}"),
+            }
+            .into());
+        }
+        let bounds = vec![(-1.0, 1.0); dim];
+        self.integrate_over(&function, &bounds, self.tolerance, 0)
     }
 
-    /// Recursive integration with refinement.
-    fn integrate_recursive<F>(&self, function: &F, level: usize) -> FeaResult<f64>
+    /// Recursive integration over one sub-box of the reference domain.
+    fn integrate_over<F>(
+        &self,
+        function: &F,
+        bounds: &[(f64, f64)],
+        tolerance: f64,
+        level: usize,
+    ) -> FeaResult<f64>
     where
         F: Fn(&NaturalCoords) -> f64,
     {
-        if level >= self.max_levels {
-            // Maximum refinement reached, use base rule
-            return Ok(self.base_rule.integrate(function));
+        let coarse = self.rule_for(bounds)?.integrate(function);
+
+        let children = Self::bisect(bounds);
+        let mut fine = 0.0;
+        for child in &children {
+            fine += self.rule_for(child)?.integrate(function);
         }
 
-        // Compute integral with current rule
-        let integral = self.base_rule.integrate(function);
+        if (coarse - fine).abs() < tolerance || level + 1 >= self.max_levels {
+            return Ok(fine);
+        }
 
-        // Estimate error by subdividing domain
-        let subdomain_integrals = self.integrate_subdomains(function)?;
-        let refined_integral: f64 = subdomain_integrals.iter().sum();
+        // Each child gets an equal share of the error budget.
+        let child_tolerance = tolerance / children.len() as f64;
+        let mut total = 0.0;
+        for child in &children {
+            total += self.integrate_over(function, child, child_tolerance, level + 1)?;
+        }
+        Ok(total)
+    }
 
-        let error = (integral - refined_integral).abs();
+    /// The 2^d sub-boxes obtained by bisecting every axis.
+    fn bisect(bounds: &[(f64, f64)]) -> Vec<Vec<(f64, f64)>> {
+        let mut boxes: Vec<Vec<(f64, f64)>> = vec![Vec::new()];
+        for &(a, b) in bounds {
+            let mid = f64::midpoint(a, b);
+            boxes = boxes
+                .into_iter()
+                .flat_map(|prefix| {
+                    [(a, mid), (mid, b)].into_iter().map(move |segment| {
+                        let mut extended = prefix.clone();
+                        extended.push(segment);
+                        extended
+                    })
+                })
+                .collect();
+        }
+        boxes
+    }
 
-        if error < self.tolerance {
-            // Converged
-            Ok(refined_integral)
-        } else {
-            // Need more refinement
-            let mut total = 0.0;
-            for _subdomain_integral in subdomain_integrals {
-                total += self.integrate_recursive(function, level + 1)?;
+    /// The base rule mapped onto the given sub-box.
+    fn rule_for(&self, bounds: &[(f64, f64)]) -> FeaResult<QuadratureRule> {
+        match bounds {
+            [(a, b)] => self.create_subdomain_rule_1d(*a, *b),
+            [(a, b), (c, d)] => self.create_subdomain_rule_2d(*a, *b, *c, *d),
+            [(a, b), (c, d), (e, f)] => self.create_subdomain_rule_3d(*a, *b, *c, *d, *e, *f),
+            _ => Err(ElementError::InvalidQuadrature {
+                reason: format!("Unsupported quadrature dimension: {}", bounds.len()),
             }
-            Ok(total)
+            .into()),
         }
     }
 
-    /// Integrate over subdomains for error estimation.
-    fn integrate_subdomains<F>(&self, function: &F) -> FeaResult<Vec<f64>>
-    where
-        F: Fn(&NaturalCoords) -> f64,
-    {
-        let mut subdomain_integrals = Vec::new();
-
-        match self.base_rule.dimension {
-            1 => {
-                // 1D: split at midpoint
-                let subdomain_rules = vec![
-                    self.create_subdomain_rule_1d(-1.0, 0.0)?,
-                    self.create_subdomain_rule_1d(0.0, 1.0)?,
-                ];
-                for rule in subdomain_rules {
-                    subdomain_integrals.push(rule.integrate(function));
-                }
-            }
-            2 => {
-                // 2D: split into 4 quadrants
-                let subdomains = vec![
-                    (-1.0, 0.0, -1.0, 0.0), // Bottom-left
-                    (0.0, 1.0, -1.0, 0.0),  // Bottom-right
-                    (-1.0, 0.0, 0.0, 1.0),  // Top-left
-                    (0.0, 1.0, 0.0, 1.0),   // Top-right
-                ];
-                for (xi_min, xi_max, eta_min, eta_max) in subdomains {
-                    let subdomain_rule =
-                        self.create_subdomain_rule_2d(xi_min, xi_max, eta_min, eta_max)?;
-                    subdomain_integrals.push(subdomain_rule.integrate(function));
-                }
-            }
-            3 => {
-                // 3D: split into 8 sub-cubes
-                let subdomains = vec![
-                    (-1.0, 0.0, -1.0, 0.0, -1.0, 0.0), // Bottom-left-front
-                    (0.0, 1.0, -1.0, 0.0, -1.0, 0.0),  // Bottom-right-front
-                    (-1.0, 0.0, 0.0, 1.0, -1.0, 0.0),  // Top-left-front
-                    (0.0, 1.0, 0.0, 1.0, -1.0, 0.0),   // Top-right-front
-                    (-1.0, 0.0, -1.0, 0.0, 0.0, 1.0),  // Bottom-left-back
-                    (0.0, 1.0, -1.0, 0.0, 0.0, 1.0),   // Bottom-right-back
-                    (-1.0, 0.0, 0.0, 1.0, 0.0, 1.0),   // Top-left-back
-                    (0.0, 1.0, 0.0, 1.0, 0.0, 1.0),    // Top-right-back
-                ];
-                for (xi_min, xi_max, eta_min, eta_max, zeta_min, zeta_max) in subdomains {
-                    let subdomain_rule = self.create_subdomain_rule_3d(
-                        xi_min, xi_max, eta_min, eta_max, zeta_min, zeta_max,
-                    )?;
-                    subdomain_integrals.push(subdomain_rule.integrate(function));
-                }
-            }
-            _ => {
-                return Err(ElementError::InvalidQuadrature {
-                    reason: format!(
-                        "Unsupported quadrature dimension: {}",
-                        self.base_rule.dimension
-                    ),
-                }
-                .into());
-            }
-        }
-
-        Ok(subdomain_integrals)
-    }
-
-    fn create_subdomain_rule_1d(&self, start: f64, end: f64) -> FeaResult<QuadratureRule> {
+    pub(super) fn create_subdomain_rule_1d(
+        &self,
+        start: f64,
+        end: f64,
+    ) -> FeaResult<QuadratureRule> {
         let scale = (end - start) / 2.0;
         let shift = f64::midpoint(start, end);
 
@@ -150,7 +136,7 @@ impl AdaptiveQuadrature {
         Ok(QuadratureRule::new(points, self.base_rule.order, 1))
     }
 
-    fn create_subdomain_rule_2d(
+    pub(super) fn create_subdomain_rule_2d(
         &self,
         xi_min: f64,
         xi_max: f64,
@@ -177,7 +163,7 @@ impl AdaptiveQuadrature {
         Ok(QuadratureRule::new(points, self.base_rule.order, 2))
     }
 
-    fn create_subdomain_rule_3d(
+    pub(super) fn create_subdomain_rule_3d(
         &self,
         xi_min: f64,
         xi_max: f64,