solver: tie weights and constraint reduction independent of the structural mesh (#358)

The tie's mortar rows take any node numbering and degree-of-freedom kinds, and
Reduction takes a DofLayout, so heat transfer can reuse both.

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

crates/solver/src/analysis.rs

@@ -343,7 +343,14 @@ pub(crate) struct Setup<'a> {
 pub(crate) fn reduction(s: &Setup, values: bool, order: &[usize], budget: Budget, warnings: &mut Vec<String>) -> Result<Reduction> {
     let mut rows = coupling_rows(s.model, s.resolved, s.mesh)?;
     rows.extend(crate::tie::rows(s.model, s.resolved, s.mesh, warnings)?);
-    Reduction::new(s.mesh, order, &rows, |active| prescribed(s.model, s.resolved, s.mesh, s.step, values, active, warnings), budget)
+    Reduction::new(
+        &s.mesh.layout,
+        s.mesh.discrete.dofs(),
+        order,
+        &rows,
+        |active| prescribed(s.model, s.resolved, s.mesh, s.step, values, active, warnings),
+        budget,
+    )
 }
 
 /// The reduced matrices' memory budget: the factorization's.

crates/solver/src/constraints.rs

@@ -52,7 +52,7 @@
 use std::collections::BTreeMap;
 
 use fieldcraft_fem::dense::inverse3;
-use fieldcraft_fem::{DofKind, DofSet};
+use fieldcraft_fem::{DofKind, DofLayout, DofSet};
 use fieldcraft_model::{CouplingKind, Model, Resolved};
 
 use crate::sparse::Csr;
@@ -201,11 +201,11 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) ->
 }
 
 /// The degrees of freedom that exist: every mesh node's, and the reference point directions
-/// some coupling row or discrete element involves.
-fn active_dofs(mesh: &Mesh, rows: &[Row]) -> Vec<bool> {
-    let mut active = vec![false; mesh.dofs()];
-    active.iter_mut().take(mesh.mesh_dofs()).for_each(|a| *a = true);
-    for dof in mesh.discrete.dofs() {
+/// some constraint row or `extra` (the discrete elements') involves.
+fn active_dofs(layout: &DofLayout, extra: impl Iterator<Item = usize>, rows: &[Row]) -> Vec<bool> {
+    let mut active = vec![false; layout.dofs()];
+    active.iter_mut().take(layout.node_dofs()).for_each(|a| *a = true);
+    for dof in extra {
         if let Some(a) = active.get_mut(dof) {
             *a = true;
         }
@@ -361,21 +361,25 @@ fn gigabytes(bytes: usize) -> String {
 impl Reduction {
     /// Eliminate the prescribed values (`prescribed(active)` gives them) and the constraint rows.
     /// `order` is the mesh nodes' order for numbering the free equations; reference point
-    /// equations follow them.
+    /// equations follow them. `layout` numbers the degrees of freedom (mesh nodes, then reference
+    /// points) and `extra` lists the reference point ones that exist without a row (discrete
+    /// elements); a heat transfer model passes its temperatures and none.
     ///
     /// A distributing coupling whose reference point is held makes one of its nodes depend on
     /// all the others, which couples every coupled degree of freedom with every other: a dense
     /// (3N)² block in the reduced matrix. Its size is checked against `budget` before anything is
     /// built.
     pub fn new(
-        mesh: &Mesh,
+        layout: &DofLayout,
+        extra: impl Iterator<Item = usize>,
         order: &[usize],
         rows: &[Row],
         prescribed: impl FnOnce(&[bool]) -> Result<Vec<Option<f64>>>,
         budget: Budget,
     ) -> Result<Reduction> {
-        let n = mesh.dofs();
-        let active = active_dofs(mesh, rows);
+        let n = layout.dofs();
+        let mesh_dofs = layout.node_dofs();
+        let active = active_dofs(layout, extra, rows);
         let prescribed = prescribed(&active)?;
         let mut state: Vec<State> = (0..n)
             .map(|d| match (active.get(d).copied().unwrap_or(false), prescribed.get(d).copied().flatten()) {
@@ -389,15 +393,15 @@ impl Reduction {
         // reference point, one block per coupling.
         let mut coupled: BTreeMap<&str, Vec<usize>> = BTreeMap::new();
         for row in rows {
-            if row.pivot >= mesh.mesh_dofs() && state.get(row.pivot) == Some(&State::Fixed) {
-                coupled.entry(row.coupling).or_default().extend(row.terms.iter().map(|t| t.0).filter(|d| *d < mesh.mesh_dofs()));
+            if row.pivot >= mesh_dofs && state.get(row.pivot) == Some(&State::Fixed) {
+                coupled.entry(row.coupling).or_default().extend(row.terms.iter().map(|t| t.0).filter(|d| *d < mesh_dofs));
             }
         }
         let mut blocks: Vec<(&str, usize, usize)> = Vec::new();
         for (name, mut dofs) in coupled {
             dofs.sort_unstable();
             dofs.dedup();
-            let mut nodes: Vec<usize> = dofs.iter().filter_map(|d| mesh.layout.locate(*d).map(|(slot, _)| slot)).collect();
+            let mut nodes: Vec<usize> = dofs.iter().filter_map(|d| layout.locate(*d).map(|(slot, _)| slot)).collect();
             nodes.dedup();
             blocks.push((name, nodes.len(), dofs.len().saturating_mul(dofs.len())));
         }
@@ -477,8 +481,8 @@ impl Reduction {
         // Equations: free mesh degrees of freedom in node order, then reference point ones.
         let mut eq = vec![usize::MAX; n];
         let mut eq_dof = Vec::new();
-        let mesh_order = order.iter().flat_map(|node| mesh.layout.range(*node));
-        for dof in mesh_order.chain(mesh.mesh_dofs()..n) {
+        let mesh_order = order.iter().flat_map(|node| layout.range(*node));
+        for dof in mesh_order.chain(mesh_dofs..n) {
             if state.get(dof) == Some(&State::Free)
                 && let Some(e) = eq.get_mut(dof)
                 && *e == usize::MAX

crates/solver/src/contact.rs

@@ -321,6 +321,12 @@ fn volume(kind: ElementKind, x: &[Vec3]) -> Option<f64> {
     }
 }
 
+/// The depth of element `kind` (corners `element`) under its face or edge `x` of measure `area`:
+/// the element's volume over the area, or the face's size where the volume is unknown.
+pub(crate) fn thickness(kind: ElementKind, element: &[Vec3], x: &[Vec3], area: f64) -> f64 {
+    volume(kind, element).filter(|v| *v > 0.0 && area > 0.0).map_or(size(x), |v| v / area)
+}
+
 /// The overlap integration of a face pair (`extent`: `None`) or an edge pair of a 2D model.
 fn overlap<T: Real>(extent: Option<Extent>, xs: &[[T; 3]], xm: &[[T; 3]], f: impl FnMut(&Sample<T>)) -> Option<()> {
     match extent {
@@ -420,7 +426,7 @@ fn surface_face(mesh: &Mesh, owner: &str, element: u32, labels: Vec<u32>, edge:
         [a, b] if edge => norm([b[0] - a[0], b[1] - a[1], 0.0]),
         _ => area_vector(&x).map_or(0.0, norm),
     };
-    let thickness = volume(el.kind, &el.xyz).filter(|v| *v > 0.0 && area > 0.0).map_or(size(&x), |v| v / area);
+    let thickness = thickness(el.kind, &el.xyz, &x, area);
     if !(thickness > 0.0 && thickness.is_finite()) {
         return Err(at(format!("the {} has no size", if edge { "edge" } else { "face" })));
     }

crates/solver/src/tie.rs

@@ -27,8 +27,14 @@
 //! face under it stays free, with a warning. The constraints are small-displacement (linear)
 //! constraints in nonlinear steps too.
 //!
-//! **Scope.** Linear solid faces (tri3, quad4) in 3D models; heat transfer and explicit dynamics
-//! refuse ties by name.
+//! **Temperature.** The same weights tie each secondary node's temperature to the main nodes'
+//! (T_A = Σ_J w_AJ T_J) in heat transfer and in the heat part of coupled temperature–displacement
+//! steps: the temperature analogue of the displacement tie, so a linear temperature field crosses
+//! a flat tie exactly, and since the constraint is eliminated (no interface resistance), the heat
+//! leaving one side enters the other.
+//!
+//! **Scope.** Linear solid faces (tri3, quad4) in 3D models; explicit dynamics refuses ties by
+//! name.
 
 use std::collections::BTreeMap;
 
@@ -38,7 +44,7 @@ use fieldcraft_fem::mortar::integrate;
 use fieldcraft_model::{Model, Resolved};
 
 use crate::constraints::Row;
-use crate::contact::{Face, faces};
+use crate::contact::faces;
 use crate::system::Mesh;
 use crate::{Error, Result};
 
@@ -50,10 +56,6 @@ const GRID_CELLS_PER_FACE: usize = 4;
 /// A secondary coefficient below this fraction of its node's total is round-off.
 const NEGLIGIBLE: f64 = 1e-12;
 
-fn corners(mesh: &Mesh, f: &Face) -> Vec<Vec3> {
-    f.nodes.iter().map(|n| mesh.xyz.get(*n).copied().unwrap_or_default()).collect()
-}
-
 /// x = M⁻¹ b for a small symmetric positive definite M (Gauss–Jordan with partial pivoting);
 /// `None` when it is singular.
 fn solve(m: &[[f64; 4]; 4], n: usize, b: [f64; 4]) -> Option<[f64; 4]> {
@@ -105,32 +107,53 @@ fn dual(x: &[Vec3]) -> Option<[[f64; 4]; 4]> {
     Some(out)
 }
 
-/// The constraint rows of every tie (ties by name, secondary nodes by label, directions in
-/// order), and warnings for the secondary nodes left free.
-pub fn rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh, warnings: &mut Vec<String>) -> Result<Vec<Row<'m>>> {
+/// One face of a tie surface: its nodes (slots of the caller's numbering, in the face's order)
+/// and the depth of the element under it (the default position tolerance is half of it).
+pub(crate) struct TieFace {
+    pub nodes: Vec<usize>,
+    pub thickness: f64,
+}
+
+/// The nodes a tie works on: their coordinates and labels by slot, and the faces of a set
+/// (`faces(owner, set)`, refusing what the analysis cannot tie).
+pub(crate) struct TieNodes<'a, F> {
+    pub xyz: &'a [Vec3],
+    pub labels: &'a [u32],
+    pub faces: F,
+}
+
+/// The constraint rows of every tie, one per secondary node and kind in `kinds` (ties by name,
+/// secondary nodes by slot, kinds in order); `dof(slot, kind)` numbers the degrees of freedom.
+/// The same mortar weights tie displacements (the translations) and temperatures (one scalar
+/// per node). Warnings name the secondary nodes left free.
+pub(crate) fn rows_of<'m, F>(
+    model: &'m Model,
+    nodes: &TieNodes<'_, F>,
+    kinds: &[DofKind],
+    dof: impl Fn(usize, DofKind) -> Result<usize>,
+    warnings: &mut Vec<String>,
+) -> Result<Vec<Row<'m>>>
+where
+    F: Fn(&str, &str) -> Result<Vec<TieFace>>,
+{
     let mut out = Vec::new();
-    if model.ties().is_empty() {
-        return Ok(out);
-    }
-    if mesh.theory.is_some() {
-        let name = model.ties().keys().next().map_or("", String::as_str);
-        return Err(Error::NotSupported(format!("ties (tie `{name}`) in 2D models")));
-    }
+    let corners = |f: &TieFace| -> Vec<Vec3> { f.nodes.iter().map(|n| nodes.xyz.get(*n).copied().unwrap_or_default()).collect() };
+    let label = |n: usize| nodes.labels.get(n).copied().unwrap_or(0);
     for (name, tie) in model.ties() {
         let owner = format!("tie `{name}`");
-        let secondary = faces(model, resolved, mesh, &owner, &tie.secondary)?;
-        let main = faces(model, resolved, mesh, &owner, &tie.main)?;
-        let main_x: Vec<Vec<Vec3>> = main.iter().map(|f| corners(mesh, f)).collect();
-        let reach = |f: &Face| tie.position_tolerance.unwrap_or(0.5 * f.thickness);
+        let secondary = (nodes.faces)(&owner, &tie.secondary)?;
+        let main = (nodes.faces)(&owner, &tie.main)?;
+        let main_x: Vec<Vec<Vec3>> = main.iter().map(corners).collect();
+        let reach = |f: &TieFace| tie.position_tolerance.unwrap_or(0.5 * f.thickness);
         let boxes: Vec<Aabb> = main.iter().zip(&main_x).map(|(f, x)| bounds(x, reach(f))).collect();
         let grid = Grid::new(&boxes, GRID_CELLS_PER_FACE * boxes.len() + 64);
-        // Per secondary node (mesh slot), its main nodes' weights M_AJ, in face order.
+        // Per secondary node (slot), its main nodes' weights M_AJ, in face order.
         let mut weights: BTreeMap<usize, BTreeMap<usize, f64>> = BTreeMap::new();
         let mut candidates = Vec::new();
         for f in &secondary {
-            let xs = corners(mesh, f);
-            let coefficients =
-                dual(&xs).ok_or_else(|| Error::Invalid(format!("{owner}: a secondary face at node {} has no area", label(mesh, f))))?;
+            let xs = corners(f);
+            let first = f.nodes.first().map_or(0, |n| label(*n));
+            let coefficients = dual(&xs).ok_or_else(|| Error::Invalid(format!("{owner}: a secondary face at node {first} has no area")))?;
             for n in &f.nodes {
                 weights.entry(*n).or_default();
             }
@@ -167,16 +190,15 @@ pub fn rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh, warnings: &m
         for (node, row) in &weights {
             let total: f64 = row.values().sum();
             if !(total.abs() > 0.0 && total.is_finite()) {
-                free.push(mesh.labels.get(*node).copied().unwrap_or(0));
+                free.push(label(*node));
                 continue;
             }
             let terms: Vec<(usize, f64)> = row.iter().map(|(j, w)| (*j, w / total)).filter(|(_, w)| w.abs() > NEGLIGIBLE).collect();
-            for k in 0..3 {
-                let kind = DofKind::translation(k).ok_or_else(|| Error::Internal("axis".into()))?;
-                let pivot = mesh.dof(*node, kind)?;
+            for &kind in kinds {
+                let pivot = dof(*node, kind)?;
                 let mut row = vec![(pivot, 1.0)];
                 for (j, w) in &terms {
-                    row.push((mesh.dof(*j, kind)?, -w));
+                    row.push((dof(*j, kind)?, -w));
                 }
                 out.push(Row { terms: row, pivot, coupling: name.as_str() });
             }
@@ -191,8 +213,25 @@ pub fn rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh, warnings: &m
     Ok(out)
 }
 
-fn label(mesh: &Mesh, f: &Face) -> u32 {
-    f.nodes.first().and_then(|n| mesh.labels.get(*n)).copied().unwrap_or(0)
+/// The displacement constraint rows of every tie in a structural analysis (three translations
+/// per secondary node).
+pub fn rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh, warnings: &mut Vec<String>) -> Result<Vec<Row<'m>>> {
+    if model.ties().is_empty() {
+        return Ok(Vec::new());
+    }
+    if mesh.theory.is_some() {
+        let name = model.ties().keys().next().map_or("", String::as_str);
+        return Err(Error::NotSupported(format!("ties (tie `{name}`) in 2D models")));
+    }
+    let nodes = TieNodes {
+        xyz: &mesh.xyz,
+        labels: &mesh.labels,
+        faces: |owner: &str, set: &str| -> Result<Vec<TieFace>> {
+            Ok(faces(model, resolved, mesh, owner, set)?.into_iter().map(|f| TieFace { nodes: f.nodes, thickness: f.thickness }).collect())
+        },
+    };
+    let translations = [DofKind::Ux, DofKind::Uy, DofKind::Uz];
+    rows_of(model, &nodes, &translations, |slot, kind| mesh.dof(slot, kind), warnings)
 }
 
 #[cfg(test)]