rtx-fea: make the analysis stack produce physics, validated against closed form

The census found rtx-fea could not produce a non-zero answer for any
analysis type. Six defects sat between a correctly specified mesh and a
natural frequency, each of which alone was fatal. Every one was found by
writing the closed-form test first and confirming red.

1. Element matrices were a stub. StandardFiniteElement::
   compute_element_matrices returned DMatrix::zeros for stiffness, force
   and mass -- and it is what GlobalAssembler calls for every element, so
   every global matrix in the crate was zero. Real quadrature-based
   stiffness and mass already existed in ElementMatrixComputer; nothing
   called them. Now wired, with the scalar mass matrix expanded by a
   Kronecker product with the spatial identity to match the interleaved
   per-node DOF layout its stiffness uses.

2. Quadrature returned no points. quadrature_rule built
   QuadratureRule::new(vec![], ..). Every integration loop iterates over
   rule.points, so an empty rule does not fail -- it skips the loop and
   yields a zero matrix. Real Gauss rules for line, triangle, quad, tet
   and hex existed unused; now dispatched by element type, with wedges as
   the triangle-line tensor product and pyramids an explicit error rather
   than an empty rule.

3. transform_derivatives computed J^-T * dN where dN is
   (num_nodes x param_dim). By the chain rule it is dN * J^-1. The two
   agree only when both are square and symmetric; for any element with
   more nodes than parametric directions -- every element -- the old form
   was a dimension mismatch that panicked inside BLAS.

4. MaterialDatabase::clone silently dropped every material, cloning
   names only, because Box<dyn Material> is not Clone. GlobalAssembler is
   constructed with materials.clone(), so every assembler ever built got
   an empty database and every analysis failed MaterialNotFound on a
   correctly specified mesh. Materials are immutable once registered, so
   the map now holds Arc and cloning shares them.

5. displacement_only numbered three displacement components on a 2-D
   mesh. Elements supply two, so assembly rejected every contribution.

6. to_dof_numbering pushed each node's DOFs in HashMap iteration order.
   When that came out [v, u] the assembler wrote the element's u row into
   the global v row. The result was still symmetric, still had the right
   rigid-body null space and still summed to the right total mass -- it
   simply described a structure with its axes transposed per node, and
   get_dof(node, DisplacementX) then pointed at the wrong row so
   constraints were applied to the wrong direction too. DofComponent now
   carries a canonical_index and the DOFs are sorted by it.

ModalAnalysis is wired to real assembly and the repaired eigensolver, and
takes boundary conditions, which it previously had no way to accept. The
eigensolver now rejects a singular stiffness explicitly: try_inverse does
not fail on a matrix singular only to working precision, so an
unconstrained structure used to return rigid-body noise dressed up as
low-frequency modes.

Validation, 18 tests:

  - Element matrices: rigid translation stores no energy, exactly 3
    rigid-body modes in 2-D and 6 in 3-D, consistent mass integrates to
    rho*V, mass positive definite, and K and M each scale only with the
    property they depend on. A zero matrix passes symmetry and
    does-not-crash checks, so these are chosen to be ones it fails.
  - Modal, end to end: longitudinal modes of a fixed-free bar against
    f_n = (2n-1)/(4L) sqrt(E/rho), within 1% on the first three, and
    second-order convergence under refinement. Axial rather than
    cantilever bending on purpose: Quad4 shear-locks, so a bending
    tolerance would fail for a reason unrelated to correctness. Bending
    is asserted as convergence from above instead, which is the honest
    claim for a locking element.

Two fixtures corrected rather than tolerances loosened: integration_tests
expected 27 DOFs for a 9-node planar mesh (3 components per node), which
encoded defect 5 and contradicted comprehensive_tdd_tests asserting
num_nodes * 2 for the same situation.

rtx-fsi stays 26/26. No new failures; the rtx-cfd quarantine is
unchanged.

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
-        let num_dofs = self.mesh.num_nodes() * 3;
-        let solution = DVector::zeros(num_dofs);
+        if self.num_modes == 0 {
+            return Err(AnalysisError::InvalidConfiguration(
+                "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);
-
-        // Add modal data
-        let eigenvalues = DVector::zeros(self.num_modes);
-        let frequencies: DVector<f64> =
-            eigenvalues.map(|lambda: f64| (lambda / (2.0 * std::f64::consts::PI)).sqrt());
-
-        results.add_data("eigenvalues".to_string(), AnalysisData::Vector(eigenvalues));
-        results.add_data("frequencies".to_string(), AnalysisData::Vector(frequencies));
+        results.add_data(
+            "eigenvalues".to_string(),
+            AnalysisData::Vector(modal.eigenvalues.clone()),
+        );
+        results.add_data(
+            "frequencies".to_string(),
+            AnalysisData::Vector(modal.frequencies.clone()),
+        );
+        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;