From 0ec188269717edafe3c74796c5641b7bd7f395c4 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Thu, 15 Oct 2026 08:25:50 -0400 Subject: [PATCH] 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 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/solver/src/analysis.rs | 9 ++- crates/solver/src/constraints.rs | 34 +++++----- crates/solver/src/contact.rs | 8 ++- crates/solver/src/tie.rs | 103 +++++++++++++++++++++---------- 4 files changed, 105 insertions(+), 49 deletions(-) diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 42804ce..40722e1 100644 --- a/crates/solver/src/analysis.rs +++ b/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) -> Result { 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. diff --git a/crates/solver/src/constraints.rs b/crates/solver/src/constraints.rs index 504eee4..3375185 100644 --- a/crates/solver/src/constraints.rs +++ b/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 { - 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, rows: &[Row]) -> Vec { + 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, order: &[usize], rows: &[Row], prescribed: impl FnOnce(&[bool]) -> Result>>, budget: Budget, ) -> Result { - 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 = (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> = 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 = dofs.iter().filter_map(|d| mesh.layout.locate(*d).map(|(slot, _)| slot)).collect(); + let mut nodes: Vec = 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 diff --git a/crates/solver/src/contact.rs b/crates/solver/src/contact.rs index 512a949..6102670 100644 --- a/crates/solver/src/contact.rs +++ b/crates/solver/src/contact.rs @@ -321,6 +321,12 @@ fn volume(kind: ElementKind, x: &[Vec3]) -> Option { } } +/// 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(extent: Option, xs: &[[T; 3]], xm: &[[T; 3]], f: impl FnMut(&Sample)) -> Option<()> { match extent { @@ -420,7 +426,7 @@ fn surface_face(mesh: &Mesh, owner: &str, element: u32, labels: Vec, 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" }))); } diff --git a/crates/solver/src/tie.rs b/crates/solver/src/tie.rs index 00d23a5..fba428c 100644 --- a/crates/solver/src/tie.rs +++ b/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 { - 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) -> Result>> { +/// 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, + 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, + warnings: &mut Vec, +) -> Result>> +where + F: Fn(&str, &str) -> Result>, +{ 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 { 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> = 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> = main.iter().map(corners).collect(); + let reach = |f: &TieFace| tie.position_tolerance.unwrap_or(0.5 * f.thickness); let boxes: Vec = 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> = 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) -> Result>> { + 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> { + 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)]