perf(engine): sparse dynamic solvers — time history and harmonic

engine/src/linalg/sparse.rs

@@ -185,6 +185,61 @@ impl CscMatrix {
         dense
     }
 
+    /// Extract the diagonal (zeros where no explicit diagonal entry exists).
+    pub fn diagonal(&self) -> Vec<f64> {
+        let mut d = vec![0.0; self.n];
+        for j in 0..self.n {
+            // Row indices within a column are sorted ascending and the
+            // diagonal is the smallest possible row, so it comes first.
+            let p = self.col_ptr[j];
+            if p < self.col_ptr[j + 1] && self.row_idx[p] == j {
+                d[j] = self.values[p];
+            }
+        }
+        d
+    }
+
+    /// Linear combination sa·A + sb·B over the union of both sparsity patterns
+    /// (e.g. Newmark's K_eff = a·K + b·M). Both inputs are lower-triangle CSC
+    /// with sorted row indices; the per-column merge is linear time. Entries
+    /// whose combined value rounds to |v| <= 1e-30 are dropped.
+    pub fn linear_combination(sa: f64, a: &CscMatrix, sb: f64, b: &CscMatrix) -> CscMatrix {
+        assert_eq!(a.n, b.n, "linear_combination: dimension mismatch");
+        let n = a.n;
+        let mut col_ptr = vec![0usize; n + 1];
+        let mut row_idx = Vec::with_capacity(a.nnz() + b.nnz());
+        let mut values = Vec::with_capacity(a.nnz() + b.nnz());
+        for j in 0..n {
+            let (mut pa, mut pb) = (a.col_ptr[j], b.col_ptr[j]);
+            let (ea, eb) = (a.col_ptr[j + 1], b.col_ptr[j + 1]);
+            while pa < ea || pb < eb {
+                let (row, v) = if pb >= eb || (pa < ea && a.row_idx[pa] < b.row_idx[pb]) {
+                    let r = a.row_idx[pa];
+                    let v = sa * a.values[pa];
+                    pa += 1;
+                    (r, v)
+                } else if pa >= ea || b.row_idx[pb] < a.row_idx[pa] {
+                    let r = b.row_idx[pb];
+                    let v = sb * b.values[pb];
+                    pb += 1;
+                    (r, v)
+                } else {
+                    let r = a.row_idx[pa];
+                    let v = sa * a.values[pa] + sb * b.values[pb];
+                    pa += 1;
+                    pb += 1;
+                    (r, v)
+                };
+                if v.abs() > 1e-30 {
+                    row_idx.push(row);
+                    values.push(v);
+                }
+            }
+            col_ptr[j + 1] = row_idx.len();
+        }
+        CscMatrix { n, col_ptr, row_idx, values }
+    }
+
     /// Extract principal submatrix for given indices (returns new CscMatrix).
     pub fn extract_principal_submatrix(&self, indices: &[usize]) -> CscMatrix {
         let m = indices.len();
@@ -371,6 +426,39 @@ mod tests {
         assert!((y[1] - 8.0).abs() < 1e-15);  // 2*1 + 3*2
     }
 
+    #[test]
+    fn test_linear_combination() {
+        // A = [[4, 2], [2, 3]], B = [[1, 0], [0, 5]] (diagonal)
+        let a = CscMatrix::from_triplets(
+            2,
+            &[0, 1, 1],
+            &[0, 0, 1],
+            &[4.0, 2.0, 3.0],
+        );
+        let b = CscMatrix::from_triplets(2, &[0, 1], &[0, 1], &[1.0, 5.0]);
+        // 2*A + 0.5*B = [[8.5, 4], [4, 8.5]]
+        let c = CscMatrix::linear_combination(2.0, &a, 0.5, &b);
+        let d = c.to_dense_symmetric();
+        assert!((d[0] - 8.5).abs() < 1e-14);
+        assert!((d[1] - 4.0).abs() < 1e-14);
+        assert!((d[3] - 8.5).abs() < 1e-14);
+
+        // Disjoint patterns: A + 0*B keeps A's pattern with B's zeros dropped
+        let e = CscMatrix::from_triplets(2, &[1], &[0], &[7.0]);
+        let f = CscMatrix::from_triplets(2, &[0, 1], &[0, 1], &[1.0, 1.0]);
+        let g = CscMatrix::linear_combination(1.0, &e, 1.0, &f);
+        assert_eq!(g.nnz(), 3);
+        let dg = g.to_dense_symmetric();
+        assert!((dg[0] - 1.0).abs() < 1e-15);
+        assert!((dg[1] - 7.0).abs() < 1e-15 && (dg[2] - 7.0).abs() < 1e-15);
+        assert!((dg[3] - 1.0).abs() < 1e-15);
+
+        // Exact cancellation drops the entry
+        let h = CscMatrix::linear_combination(1.0, &a, -2.0, &a); // = -A
+        let i = CscMatrix::linear_combination(1.0, &a, 1.0, &h); // A + (-A) = 0
+        assert_eq!(i.nnz(), 0, "A - A should have an empty pattern");
+    }
+
     #[test]
     fn test_sym_mat_vec_vs_dense() {
         // 3×3 SPD: [[10, 2, 1], [2, 8, 3], [1, 3, 6]]

engine/src/solver/assembly.rs

@@ -1818,6 +1818,9 @@ fn assemble_element_loads_3d_mapped(
 /// Sparse assembly result: CSC lower-triangle Kff + dense force vector.
 pub struct SparseAssemblyResult {
     pub k_ff: CscMatrix,
+    /// Full n×n K (all DOFs), only when requested — for reactions and
+    /// prescribed-DOF corrections without a dense assembly.
+    pub k_full: Option<CscMatrix>,
     pub f: Vec<f64>,       // n_total force vector (same as dense)
     pub max_diag_k: f64,
     pub artificial_dofs: Vec<usize>,
@@ -1836,6 +1839,14 @@ pub struct SparseAssemblyResult3D {
 
 /// Assemble sparse Kff for 2D. Returns CSC lower-triangle of the free-DOF block.
 pub fn assemble_sparse_2d(input: &SolverInput, dof_num: &DofNumbering) -> SparseAssemblyResult {
+    assemble_sparse_2d_ex(input, dof_num, false)
+}
+
+/// Sparse 2D assembly with optional full-K (all DOFs) alongside the free block.
+///
+/// Triplets are collected unfiltered (the free×free filter is applied when
+/// building k_ff); with `build_k_full` the full n×n matrix is kept too.
+pub fn assemble_sparse_2d_ex(input: &SolverInput, dof_num: &DofNumbering, build_k_full: bool) -> SparseAssemblyResult {
     let n = dof_num.n_total;
     let nf = dof_num.n_free;
     let mut f_global = vec![0.0; n];
@@ -1876,9 +1887,7 @@ pub fn assemble_sparse_2d(input: &SolverInput, dof_num: &DofNumbering) -> Sparse
                 dof_num.global_dof(elem.node_j, 1).unwrap(),
             ];
             for i in 0..4 {
-                if truss_dofs[i] >= nf { continue; }
                 for j in 0..4 {
-                    if truss_dofs[j] >= nf { continue; }
                     let gi = truss_dofs[i];
                     let gj = truss_dofs[j];
                     if gi >= gj {
@@ -1887,7 +1896,7 @@ pub fn assemble_sparse_2d(input: &SolverInput, dof_num: &DofNumbering) -> Sparse
                         trip_vals.push(k_elem[i * 4 + j]);
                     }
                 }
-                diag_vals[truss_dofs[i]] += k_elem[i * 4 + i];
+                if truss_dofs[i] < nf { diag_vals[truss_dofs[i]] += k_elem[i * 4 + i]; }
             }
 
             // Assemble thermal FEF for 2D truss elements (sparse path)
@@ -1921,9 +1930,7 @@ pub fn assemble_sparse_2d(input: &SolverInput, dof_num: &DofNumbering) -> Sparse
             let ndof = elem_dofs.len();
 
             for i in 0..ndof {
-                if elem_dofs[i] >= nf { continue; }
                 for j in 0..ndof {
-                    if elem_dofs[j] >= nf { continue; }
                     let gi = elem_dofs[i];
                     let gj = elem_dofs[j];
                     if gi >= gj {
@@ -1932,7 +1939,7 @@ pub fn assemble_sparse_2d(input: &SolverInput, dof_num: &DofNumbering) -> Sparse
                         trip_vals.push(k_glob[i * ndof + j]);
                     }
                 }
-                diag_vals[elem_dofs[i]] += k_glob[i * ndof + i];
+                if elem_dofs[i] < nf { diag_vals[elem_dofs[i]] += k_glob[i * ndof + i]; }
             }
 
             let load_refs: Vec<&SolverLoad> = input.loads.iter().collect();
@@ -1956,9 +1963,7 @@ pub fn assemble_sparse_2d(input: &SolverInput, dof_num: &DofNumbering) -> Sparse
         let dofs = dof_num.element_dofs(conn.node_i, conn.node_j);
         let ndof = dofs.len();
         for i in 0..ndof {
-            if dofs[i] >= nf { continue; }
             for j in 0..ndof {
-                if dofs[j] >= nf { continue; }
                 let gi = dofs[i];
                 let gj = dofs[j];
                 if gi >= gj {
@@ -1967,7 +1972,7 @@ pub fn assemble_sparse_2d(input: &SolverInput, dof_num: &DofNumbering) -> Sparse
                     trip_vals.push(ke[i * 6 + j]);
                 }
             }
-            diag_vals[dofs[i]] += ke[i * 6 + i];
+            if dofs[i] < nf { diag_vals[dofs[i]] += ke[i * 6 + i]; }
         }
     }
 
@@ -2040,8 +2045,29 @@ pub fn assemble_sparse_2d(input: &SolverInput, dof_num: &DofNumbering) -> Sparse
         }
     }
 
-    let k_ff = CscMatrix::from_triplets(nf, &trip_rows, &trip_cols, &trip_vals);
-    SparseAssemblyResult { k_ff, f: f_global, max_diag_k: max_diag, artificial_dofs }
+    let k_full = if build_k_full {
+        Some(CscMatrix::from_triplets(n, &trip_rows, &trip_cols, &trip_vals))
+    } else {
+        None
+    };
+    // The triplets cover all DOFs; k_ff keeps only the free×free block.
+    let (ff_rows, ff_cols, ff_vals): (Vec<usize>, Vec<usize>, Vec<f64>) = if n == nf {
+        (trip_rows, trip_cols, trip_vals)
+    } else {
+        let mut fr = Vec::new();
+        let mut fc = Vec::new();
+        let mut fv = Vec::new();
+        for t in 0..trip_rows.len() {
+            if trip_rows[t] < nf && trip_cols[t] < nf {
+                fr.push(trip_rows[t]);
+                fc.push(trip_cols[t]);
+                fv.push(trip_vals[t]);
+            }
+        }
+        (fr, fc, fv)
+    };
+    let k_ff = CscMatrix::from_triplets(nf, &ff_rows, &ff_cols, &ff_vals);
+    SparseAssemblyResult { k_ff, k_full, f: f_global, max_diag_k: max_diag, artificial_dofs }
 }
 
 /// Apply inclined support rotation to COO triplets and force vector.