diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index c344484..ac2947c 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -215,6 +215,9 @@ impl System { if let Some(name) = model.couplings().keys().next() { return Err(Error::NotSupported(format!("couplings (coupling `{name}`) in explicit dynamics"))); } + if let Some(name) = model.ties().keys().next() { + return Err(Error::NotSupported(format!("ties (tie `{name}`) in explicit dynamics"))); + } // Elements by label, with their section's material. let mut materials: Vec = Vec::new(); diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 496b7c3..42804ce 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -341,7 +341,8 @@ pub(crate) struct Setup<'a> { } pub(crate) fn reduction(s: &Setup, values: bool, order: &[usize], budget: Budget, warnings: &mut Vec) -> Result { - let rows = coupling_rows(s.model, s.resolved, s.mesh)?; + 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) } diff --git a/crates/solver/src/contact.rs b/crates/solver/src/contact.rs index 41e7bf1..1e08b2e 100644 --- a/crates/solver/src/contact.rs +++ b/crates/solver/src/contact.rs @@ -128,14 +128,14 @@ const GRAIN: usize = 16; /// A face of a contact surface: its corners (mesh slots) and their translations. #[derive(Debug, Clone)] -struct Face { - nodes: Vec, - dofs: Vec<[usize; 3]>, +pub(crate) struct Face { + pub(crate) nodes: Vec, + pub(crate) dofs: Vec<[usize; 3]>, /// The nodes of its element (never in contact with the face). owner: Vec, /// The penetration beyond which a point has passed through, and its element's thickness. depth: f64, - thickness: f64, + pub(crate) thickness: f64, /// The penalty per length its element sets (before the penalty scale). stiffness: f64, } @@ -303,17 +303,18 @@ fn size(x: &[Vec3]) -> f64 { (0..3).map(|k| b[1][k] - b[0][k]).fold(0.0, f64::max) } -/// The faces of a set, refused unless they are linear faces of analysed solid elements. -fn faces(model: &Model, resolved: &Resolved, mesh: &Mesh, name: &str, set: &str) -> Result> { +/// The faces of set `set` for definition `owner` (``contact `name` ``, ``tie `name` ``), refused +/// unless they are linear faces of analysed solid elements. +pub(crate) fn faces(model: &Model, resolved: &Resolved, mesh: &Mesh, owner: &str, set: &str) -> Result> { let list = resolved.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default(); let mut out = Vec::with_capacity(list.len()); for f in list { - let at = |what: String| Error::Invalid(format!("contact `{name}`: element {} face {}: {what}", f.element, f.face)); + let at = |what: String| Error::Invalid(format!("{owner}: element {} face {}: {what}", f.element, f.face)); let el = mesh.element(f.element).ok_or_else(|| at("the element is not a solid element of the analysis".into()))?; let labels = model.face_nodes(*f).ok_or_else(|| at("no such face".into()))?; if !matches!(labels.len(), 3 | 4) { return Err(Error::NotSupported(format!( - "contact `{name}` on the faces of {} elements (element {}) in implicit steps; it takes the linear faces of tet4, wedge6 and hex8 elements", + "{owner} on the faces of {} elements (element {}) in implicit analyses; it takes the linear faces of tet4, wedge6 and hex8 elements", el.kind.name(), f.element ))); @@ -463,8 +464,9 @@ impl Interface { "general (single-surface) contact `{name}` in implicit steps; define contact pairs (interaction.contact.pair), or solve it in explicit dynamics" ))); }; - let secondary = faces(model, resolved, mesh, name, secondary)?; - let main = faces(model, resolved, mesh, name, main)?; + let owner = format!("contact `{name}`"); + let secondary = faces(model, resolved, mesh, &owner, secondary)?; + let main = faces(model, resolved, mesh, &owner, main)?; let mut nodes: Vec = secondary.iter().flat_map(|f| f.nodes.iter().copied()).collect(); nodes.sort_unstable(); nodes.dedup(); diff --git a/crates/solver/src/heat.rs b/crates/solver/src/heat.rs index 9b676fa..198aed0 100644 --- a/crates/solver/src/heat.rs +++ b/crates/solver/src/heat.rs @@ -414,15 +414,25 @@ pub(crate) fn temperatures(model: &Model, resolved: &Resolved, step: &Step, limi Ok((s.mesh.labels, t)) } +/// Heat transfer does not tie temperatures across a tie yet. +fn refuse_ties(model: &Model) -> Result<()> { + match model.ties().keys().next() { + Some(name) => Err(crate::Error::NotSupported(format!("ties (tie `{name}`) in heat transfer"))), + None => Ok(()), + } +} + /// The steady temperatures of a step's thermal conditions at the nodes `labels` (sorted), `None` /// at the nodes conduction does not reach. pub(crate) fn temperatures_at(model: &Model, resolved: &Resolved, step: &Step, limits: &Limits, labels: &[u32]) -> Result>> { + refuse_ties(model)?; let (solved, t) = temperatures(model, resolved, step, limits)?; Ok(labels.iter().map(|l| solved.binary_search(l).ok().and_then(|i| t.get(i)).copied()).collect()) } /// Solve a heat transfer step: the step named `step`, or the model's only one. pub fn solve_heat(model: &Model, step: Option<&str>, limits: &Limits) -> Result { + refuse_ties(model)?; let resolved = model.resolve()?; let step = model.select_step(step, |p| matches!(p, Procedure::HeatTransfer { .. }), "heat transfer")?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index d55a534..bb46285 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -40,6 +40,7 @@ pub mod nonlinear; pub mod ordering; pub mod sparse; pub mod system; +pub mod tie; pub use analysis::{ ElementStress, FrequencyOptions, FrequencyResult, Location, MassKind, Mode, PointResult, StaticResult, Stats, StepResult, solve_frequency, diff --git a/crates/solver/src/tie.rs b/crates/solver/src/tie.rs new file mode 100644 index 0000000..00d23a5 --- /dev/null +++ b/crates/solver/src/tie.rs @@ -0,0 +1,228 @@ +//! Ties: non-matching meshes bonded by mortar constraints with dual shape functions. +//! +//! Each node A of a tie's secondary surface follows the main surface under it through the +//! mortar condition ∫ ψ_A (u_s − u_m) dA = 0 per direction, integrated over the overlaps of the +//! secondary faces with the main faces in the reference configuration ([`fieldcraft_fem::mortar`]). +//! The multiplier space uses the dual shape functions of each secondary face, ψ_a = Σ_b A_ab N_b +//! with A = D M⁻¹ (M_ab = ∫ N_a N_b over the face, D = diag ∫ N_a), biorthogonal to the face's +//! own shape functions, ∫ ψ_a N_b = δ_ab ∫ N_b (B. I. Wohlmuth, *A mortar finite element method +//! using dual spaces for the Lagrange multiplier*, SIAM J. Numer. Anal. 38(3):989–1012, 2000; +//! A. Popp, B. I. Wohlmuth, M. W. Gee, W. A. Wall, *Dual quadratic mortar finite element methods +//! for 3D finite deformation contact*, SIAM J. Sci. Comput. 34(4):B421–B446, 2012). The mortar +//! matrix of the secondary side is then diagonal, and the condition makes each secondary node a +//! weighted average of the main nodes under it, +//! +//! u_A = Σ_J (M_AJ / Σ_K M_AK) u_J, M_AJ = Σ_overlaps ∫ ψ_A N_J dA, +//! +//! a linear constraint row per direction that the constraint elimination takes like a coupling's +//! ([`crate::constraints`]): the secondary node is the intended dependent. Where a secondary +//! face is wholly covered, Σ_K M_AK = ∫ N_A, the plain mortar condition; dividing by the sum keeps +//! rigid motions exact at the surface's edges, where it is not. A uniform stress crosses a tie on +//! flat faces exactly (the tie patch test), because the dual mortar conditions are satisfied by +//! the linear field on both sides. +//! +//! **Pairing.** Each secondary face takes the main faces whose boxes, grown by the position +//! tolerance, overlap its own, and that it overlaps along its normal with a mean gap within the +//! tolerance (by default half the thickness of the main element). A secondary node with no main +//! 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. + +use std::collections::BTreeMap; + +use fieldcraft_fem::DofKind; +use fieldcraft_fem::contact::{Aabb, Grid, bounds}; +use fieldcraft_fem::mortar::integrate; +use fieldcraft_model::{Model, Resolved}; + +use crate::constraints::Row; +use crate::contact::{Face, faces}; +use crate::system::Mesh; +use crate::{Error, Result}; + +type Vec3 = fieldcraft_fem::frame::Vec3; + +/// Faces per main grid cell budget (the grid's cap is this times the faces). +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]> { + let mut a = *m; + let mut x = b; + for c in 0..n { + let p = (c..n).max_by(|i, j| a[*i][c].abs().total_cmp(&a[*j][c].abs()))?; + if !(a[p][c].abs() > 0.0) { + return None; + } + a.swap(c, p); + x.swap(c, p); + for r in 0..n { + if r != c { + let f = a[r][c] / a[c][c]; + for k in c..n { + a[r][k] -= f * a[c][k]; + } + x[r] -= f * x[c]; + } + } + } + for c in 0..n { + x[c] /= a[c][c]; + } + x.iter().all(|v| v.is_finite()).then_some(x) +} + +/// The dual coefficients A = D M⁻¹ of a face (rows a: ψ_a = Σ_b A_ab N_b). +fn dual(x: &[Vec3]) -> Option<[[f64; 4]; 4]> { + let n = x.len(); + let mirror: Vec = x.iter().rev().copied().collect(); + let mut mass = [[0.0; 4]; 4]; + integrate(x, &mirror, |p| { + for a in 0..4 { + for b in 0..4 { + mass[a][b] += p.weight * p.secondary[a] * p.secondary[b]; + } + } + })?; + // A M = D, M symmetric: row a of A solves M (row a)ᵀ = D_aa e_a. + let mut out = [[0.0; 4]; 4]; + for a in 0..n { + let d: f64 = mass[a].iter().take(n).sum(); + let mut e = [0.0; 4]; + e[a] = d; + out[a] = solve(&mass, n, e)?; + } + 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>> { + 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"))); + } + 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 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. + 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))))?; + for n in &f.nodes { + weights.entry(*n).or_default(); + } + grid.candidates_in(&bounds(&xs, 0.0), &mut candidates); + for &c in &candidates { + let (Some(m), Some(xm)) = (main.get(c as usize), main_x.get(c as usize)) else { continue }; + if m.nodes.iter().any(|n| f.nodes.contains(n)) { + continue; + } + let mut local = [[0.0; 4]; 4]; + let (mut area, mut gap) = (0.0, 0.0); + integrate(&xs, xm, |p| { + area += p.weight; + gap += p.weight * p.gap().abs(); + for a in 0..4 { + let psi: f64 = (0..4).map(|b| coefficients[a][b] * p.secondary[b]).sum(); + for j in 0..4 { + local[a][j] += p.weight * psi * p.main[j]; + } + } + }); + if !(area > 0.0) || gap / area > reach(m) { + continue; + } + for (a, s) in f.nodes.iter().enumerate() { + let row = weights.entry(*s).or_default(); + for (j, mn) in m.nodes.iter().enumerate() { + *row.entry(*mn).or_insert(0.0) += local[a][j]; + } + } + } + } + let mut free = Vec::new(); + 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)); + 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)?; + let mut row = vec![(pivot, 1.0)]; + for (j, w) in &terms { + row.push((mesh.dof(*j, kind)?, -w)); + } + out.push(Row { terms: row, pivot, coupling: name.as_str() }); + } + } + if let Some(first) = free.first() { + warnings.push(format!( + "{owner}: {} secondary node(s) have no main face within the position tolerance and stay free (node {first} first)", + free.len() + )); + } + } + Ok(out) +} + +fn label(mesh: &Mesh, f: &Face) -> u32 { + f.nodes.first().and_then(|n| mesh.labels.get(*n)).copied().unwrap_or(0) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn dual_functions_are_biorthogonal_to_the_face_s_own() { + let quad = [[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [2.2, 1.0, 0.0], [0.1, 1.3, 0.0]]; + for face in [&quad[..], &quad[..3]] { + let a = dual(face).unwrap(); + let mirror: Vec = face.iter().rev().copied().collect(); + let mut check = [[0.0; 4]; 4]; + let mut own = [0.0; 4]; + integrate(face, &mirror, |p| { + for i in 0..4 { + own[i] += p.weight * p.secondary[i]; + let psi: f64 = (0..4).map(|b| a[i][b] * p.secondary[b]).sum(); + for j in 0..4 { + check[i][j] += p.weight * psi * p.secondary[j]; + } + } + }); + for i in 0..face.len() { + for j in 0..face.len() { + let want = if i == j { own[i] } else { 0.0 }; + assert!((check[i][j] - want).abs() < 1e-12, "{i} {j}: {} vs {want}", check[i][j]); + } + } + } + assert!(dual(&[[0.0; 3], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]]).is_none()); + } +} diff --git a/crates/solver/tests/contact.rs b/crates/solver/tests/contact.rs index 9f93692..b9e7294 100644 --- a/crates/solver/tests/contact.rs +++ b/crates/solver/tests/contact.rs @@ -463,3 +463,100 @@ fn coulomb_friction_sticks_below_the_cone_and_slips_on_it() { assert!((t3 - 0.2 * n3).abs() < tol * 0.2 * n3, "finite {finite}: T {t3} vs μk N {}", 0.2 * n3); } } + +/// Two blocks of non-matching meshes tied at z = 1 (`lower_kind` 3 × 3 × 2, graded or distorted, +/// below; `upper_kind` 2 × 2 × 2 above), every outer boundary node moved by the linear field +/// u = G x: the tie patch test. +fn tied(lower_kind: ElementKind, upper_kind: ElementKind, distorted: bool, procedure: Procedure) -> (Blocks, [[f64; 3]; 3]) { + let mut b = Blocks::default(); + let low = if distorted { + b.block(lower_kind, [3, 3, 2], |p| [p[0] + 0.3 * p[0] * (1.0 - p[0]) * (p[1] - 0.5), p[1] + 0.3 * p[1] * (1.0 - p[1]) * (p[0] - 0.4), p[2]]) + } else { + b.block(lower_kind, [3, 3, 2], |p| [p[0].powf(1.4), p[1].powf(0.8), p[2]]) + }; + let up = b.block(upper_kind, [2, 2, 2], |p| [p[0], p[1], 1.0 + p[2]]); + let all: Vec = b.model.elements().keys().copied().collect(); + b.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap(); + b.model.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: 1000.0, poisson: 0.3 }), ..Material::default() }).unwrap(); + b.model + .add_section("solid", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() }) + .unwrap(); + b.model.add_step(Step::new("Step", procedure), None).unwrap(); + b.face_set("lower_top", &low, |x| near(x[2], 1.0)); + b.face_set("upper_bottom", &up, |x| near(x[2], 1.0)); + b.model.add_tie("glue", Tie::new("lower_top", "upper_bottom")).unwrap(); + let g = [[1e-3, 2e-4, -3e-4], [4e-4, -5e-4, 1e-4], [-2e-4, 3e-4, 2e-3]]; + let outer = |x: [f64; 3]| near(x[0], 0.0) || near(x[0], 1.0) || near(x[1], 0.0) || near(x[1], 1.0) || near(x[2], 0.0) || near(x[2], 2.0); + for l in b.nodes(outer) { + let x = b.xyz[&l]; + let u: [f64; 3] = std::array::from_fn(|i| (0..3).map(|j| g[i][j] * x[j]).sum()); + let name = format!("n{l}"); + b.model.add_set(&name, Set { nodes: vec![l], ..Set::default() }).unwrap(); + b.model + .add_boundary_condition( + &name, + BoundaryCondition::Displacement { target: Target::Set(name.clone()), values: [Some(u[0]), Some(u[1]), Some(u[2]), None, None, None] }, + ) + .unwrap(); + b.model.set_boundary_condition_active("Step", &name, true).unwrap(); + } + (b, g) +} + +/// The tie patch test (dual mortar: Wohlmuth 2000): a linear displacement field across a tie +/// of non-matching meshes is reproduced exactly, so every element has the same stress, in linear +/// static, nonlinear static and through the interface nodes the tie moves. +#[test] +fn a_tie_passes_the_patch_test() { + for (lower, upper, distorted, name) in [ + (ElementKind::Hex8, ElementKind::Hex8, false, "hex8 on graded hex8"), + (ElementKind::Hex8, ElementKind::Tet4, false, "tet4 on hex8"), + (ElementKind::Tet4, ElementKind::Hex8, false, "hex8 on tet4"), + (ElementKind::Hex8, ElementKind::Hex8, true, "hex8 on distorted hex8"), + ] { + let (b, g) = tied(lower, upper, distorted, Procedure::Static {}); + let r = fieldcraft_solver::solve_static(&b.model, Some("Step"), &Limits::default()).unwrap(); + // The interface nodes follow the field too. + let tol = if distorted { 1e-6 } else { 1e-11 }; + for l in b.nodes(|x| near(x[2], 1.0)) { + let x = b.xyz[&l]; + let u = r.displacement[r.nodes.binary_search(&l).unwrap()]; + for i in 0..3 { + let want: f64 = (0..3).map(|j| g[i][j] * x[j]).sum(); + assert!((u[i] - want).abs() < tol * 1e-2, "{name}: node {l} u {u:?}"); + } + } + let first = r.elements[0].nodal[0]; + for e in &r.elements { + for s in &e.nodal { + for k in 0..6 { + assert!((s[k] - first[k]).abs() < tol * 10.0, "{name}: element {} {s:?} vs {first:?}", e.label); + } + } + } + assert!(r.warnings.iter().all(|w| !w.contains("stay free")), "{name}: {:?}", r.warnings); + } + // The same in a nonlinear step with small strain. + let (b, _) = tied( + ElementKind::Hex8, + ElementKind::Hex8, + false, + Procedure::NonlinearStatic(NonlinearStatic { nonlinear_geometry: false, ..NonlinearStatic::default() }), + ); + let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap(); + let first = r.elements[0].nodal[0]; + assert!(r.elements.iter().flat_map(|e| &e.nodal).all(|s| (0..6).all(|k| (s[k] - first[k]).abs() < 1e-8))); +} + +#[test] +fn what_ties_refuse() { + let (mut b, _) = tied(ElementKind::Hex8, ElementKind::Hex8, false, Procedure::Static {}); + b.model.add_step(Step::new("Heat", Procedure::HeatTransfer { transient: None }), None).unwrap(); + let e = fieldcraft_solver::solve_heat(&b.model, Some("Heat"), &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("tie `glue`") && e.contains("heat transfer"), "{e}"); + let (mut b, _) = tied(ElementKind::Hex8, ElementKind::Hex8, false, Procedure::Static {}); + // (Touching surfaces are tied within any tolerance.) + b.model.replace_tie("glue", Tie { position_tolerance: Some(1e-9), ..Tie::new("lower_top", "upper_bottom") }).unwrap(); + let r = fieldcraft_solver::solve_static(&b.model, Some("Step"), &Limits::default()).unwrap(); + assert!(r.warnings.iter().all(|w| !w.contains("stay free")), "{:?}", r.warnings); +}