diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index ec57169..496b7c3 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -113,6 +113,8 @@ pub struct StaticResult { /// The increments of a nonlinear static step (its stresses are Cauchy stresses and its /// strains logarithmic strains). pub nonlinear: Option, + /// Each contact's state at the end of a nonlinear static step, by name. + pub contacts: Vec, } /// One natural mode. @@ -362,6 +364,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul |p| matches!(p, Procedure::Static {} | Procedure::CoupledTemperatureDisplacement {}), "static or coupled temperature–displacement", )?; + crate::system::refuse_contact(model, "linear static")?; let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).chain(mesh.warnings.iter().cloned()).collect(); warnings.extend(inelastic_warnings(model)); @@ -497,7 +500,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul } let stress: Vec<[f64; 6]> = sum.iter().zip(&count).map(|(s, c)| s.map(|v| if *c > 0 { v / f64::from(*c) } else { 0.0 })).collect(); let von_mises = stress.iter().map(von_mises).collect(); - let scale = applied_total.iter().chain(&reaction_total).fold(0.0f64, |m, v| m.max(v.abs())); + let scale = crate::system::balance_scale(&applied_total, &reaction_total, &reaction); let ground = mesh.discrete.ground_force(&u, &mesh); let imbalance = (0..3).map(|c| (applied_total[c] + reaction_total[c] + ground[c]).abs()).fold(0.0, f64::max); if scale > 0.0 && imbalance > 1e-6 * scale { @@ -522,6 +525,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul warnings, temperature, nonlinear: None, + contacts: Vec::new(), }) } @@ -568,6 +572,7 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt { return Err(Error::Invalid(format!("the shift must be a finite number, got {s}"))); } + crate::system::refuse_contact(model, "frequency")?; let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).chain(mesh.warnings.iter().cloned()).collect(); warnings.extend(inelastic_warnings(model)); diff --git a/crates/solver/src/buckle.rs b/crates/solver/src/buckle.rs index d5d21dd..52a921c 100644 --- a/crates/solver/src/buckle.rs +++ b/crates/solver/src/buckle.rs @@ -113,6 +113,7 @@ pub fn solve_buckle(model: &Model, step: Option<&str>, limits: &Limits) -> Resul }; step.procedure.validate()?; let resolved = model.resolve()?; + crate::system::refuse_contact(model, "buckling")?; let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; if let Some(t) = mesh.trusses.first() { return Err(Error::NotSupported(format!( diff --git a/crates/solver/src/contact.rs b/crates/solver/src/contact.rs new file mode 100644 index 0000000..d6678f0 --- /dev/null +++ b/crates/solver/src/contact.rs @@ -0,0 +1,751 @@ +//! Contact in nonlinear static steps: mortar (segment-to-segment) normal contact enforced by an +//! augmented Lagrangian, with its consistent tangent. +//! +//! **Discretisation: mortar, one pass.** A contact pair's `secondary` surface is the non-mortar +//! side, its `main` surface the mortar side (M. A. Puso, T. A. Laursen, *A mortar +//! segment-to-segment contact method for large deformation solid mechanics*, CMAME 193:601–629, +//! 2004). Each secondary node A carries one normal pressure p_A and one weighted gap +//! +//! g̃_A = Σ_segments ∫ Φ_A (x_m − x_s) · n dA, +//! +//! integrated over the overlaps of its faces with the main faces they face +//! ([`fieldcraft_fem::mortar`]: projection along the secondary face's normal n, clipping, +//! seven-point rule), with Φ_A the secondary face's own shape functions. The pressure field is +//! interpolated like the displacements, so a uniform pressure across non-matching meshes is +//! transmitted exactly: both sides receive the consistent nodal forces ∫ p N dA of it, which is +//! the contact patch test that node-to-surface contact fails (R. L. Taylor, P. Papadopoulos, +//! *On a patch test for contact problems in two dimensions*, 1991; T. W. McDevitt, T. A. +//! Laursen, *A mortar-finite element formulation for frictional contact problems*, IJNME +//! 48:1525–1547, 2000). A node's gap is normalised by its reference area A_A = ∫ Φ_A dA over its +//! own faces: ĝ_A = g̃_A / A_A, a length (negative: penetration). +//! +//! **Enforcement: augmented Lagrangian** (J. C. Simo, T. A. Laursen, *An augmented Lagrangian +//! treatment of contact problems involving friction*, Comput. Struct. 42(1):97–116, 1992; +//! P. Alart, A. Curnier, *A mixed formulation for frictional contact problems prone to Newton +//! like solution methods*, CMAME 92:353–375, 1991). With the node's multiplier λ_A held fixed, +//! the pressure is p_A = ⟨λ_A − ε_A ĝ_A⟩ (⟨x⟩ = max(x, 0)), the derivative of the potential +//! Π = Σ_A A_A/(2 ε_A) ⟨λ_A − ε_A ĝ_A⟩², whose forces −p_A ∂g̃_A/∂x enter the internal forces and +//! whose Hessian +//! +//! K_c = Σ_{A closed} (ε_A / A_A) ∂g̃_A ⊗ ∂g̃_A − p_A ∂²g̃_A +//! +//! the tangent: the consistent linearisation of the segmentation, the normal and the projection +//! (the second term is the geometric stiffness of the contact, which matters on curved and +//! sliding surfaces), both from the automatic differentiation of the mortar integrals +//! ([`fieldcraft_fem::dual`]). The penalty only regularises: once Newton has converged an +//! increment, the multipliers take the pressures (Uzawa's update λ_A ← p_A) and the increment is +//! solved again from its solution, until every closed node's gap is within [`GAP_TOLERANCE`] of +//! the thickness of the element under it and no open node keeps a multiplier +//! ([`Interface::augmentation_error`]; at most [`MAX_AUGMENTATIONS`] times, then with a warning). +//! The penalty per length is ε_A = s · 10 · 3Ḡ / h +//! over the secondary faces at the node (Ḡ the mean shear modulus of the element under the face, +//! h its thickness V/A, s the contact's `penalty_scale`): about ten times the stiffness of the +//! element under the node, stiff enough that a few augmentations reach the tolerance and soft +//! enough to leave the tangent well conditioned (3Ḡ rather than the constrained modulus keeps it +//! so for nearly incompressible rubber). +//! +//! **Active set.** The pressure's max() is the active set: each evaluation decides which nodes +//! are closed (p_A > 0) from the current gaps and multipliers, and Newton runs on the +//! semismooth system (Alart & Curnier 1991; A. Popp, M. W. Gee, W. A. Wall, *A finite +//! deformation mortar contact formulation using a primal–dual active set strategy*, IJNME +//! 79:1354–1391, 2009). The nonlinear solver cuts the increment back when the set alternates +//! between two states ([`CHATTER_LIMIT`] times: chattering). +//! +//! **Small and finite sliding.** Without nonlinear geometry the face pairs (the main faces within +//! a main element's thickness of each secondary face) and their mortar integrals +//! D_ai = ∫ Φ_a N_i and M_aj = ∫ Φ_a N_j are found once in the reference configuration, so that +//! g̃_a = n · (Σ_j M_aj x_j − Σ_i D_ai x_i) is linear in the positions (small sliding: exact, and +//! consistent with the small-strain elements); with nonlinear geometry they are searched and +//! integrated on the current configuration at every evaluation (finite sliding). +//! +//! **Search.** Every evaluation of finite sliding: the main faces' boxes, grown by their depth limit, go into a +//! bucket grid ([`fieldcraft_fem::contact::Grid`]); each secondary face takes the main faces in +//! the cells its box overlaps that face it (unit normals more than [`FACING`] from parallel) and +//! share no node with its element. A pair whose overlap reaches deeper than [`DEPTH_FRACTION`] +//! of the main element's thickness has passed through and is left out. The pairs are integrated +//! in parallel and summed in face order, so results do not depend on the thread count. +//! +//! **Initial overclosure.** Before the step, nodes whose gap starts negative are reported: +//! `adjust` (the default) takes the initial overclosure out of their gap (an offset of the +//! weighted gap: strain-free, the nodes are not moved), `refuse` refuses the step naming the +//! first node and its depth. +//! +//! **Output.** Per secondary node the pressure as the node's normal contact force over the area +//! it carries (∫ Φ_A over the overlaps: smoother than the multipliers, which are the coefficients +//! of the pressure field and ring at the contact's edge), the multiplier, the gap and the status; +//! per contact the force on the secondary surface and its normal part. +//! +//! **Scope.** Pairs of surfaces of linear solid faces (tri3 and quad4 faces: tet4, wedge6 and +//! hex8 elements) in 3D models, frictionless. General (single-surface) contact, quadratic faces, +//! 2D models and arc-length steps are refused by name. + +use std::collections::BTreeMap; + +use fieldcraft_fem::contact::{Aabb, Grid, bounds}; +use fieldcraft_fem::dual::{Dual, Hyper, Real, VARS}; +use fieldcraft_fem::frame::{Vec3, dot, norm, unit}; +use fieldcraft_fem::mortar::{area_vector, integrate}; +use fieldcraft_model::{Contact, ContactSurfaces, ElementKind, InitialPenetration, Model, Resolved}; +use fieldcraft_parallel as parallel; +use serde::Serialize; + +use crate::system::{Mesh, translations}; +use crate::{Error, Result}; + +/// Facing faces: the dot product of their unit normals is below minus this. +pub const FACING: f64 = 0.1; + +/// A pair deeper than this fraction of the main element's thickness has passed through. +pub const DEPTH_FRACTION: f64 = 0.5; + +/// ε = s · PENALTY_FACTOR · 3Ḡ / h (see the module docs). +pub const PENALTY_FACTOR: f64 = 10.0; + +/// Augmentations stop when every closed node's gap is within this fraction of the thickness of +/// the element under it (and no open node keeps a multiplier). +pub const GAP_TOLERANCE: f64 = 1e-5; + +/// Most augmentations of one increment. +pub const MAX_AUGMENTATIONS: usize = 25; + +/// An increment is cut back when the active set alternates between two states this many times. +pub const CHATTER_LIMIT: usize = 3; + +/// Most faces over the contacts of a model. +pub const MAX_CONTACT_FACES: usize = 2_000_000; + +/// Secondary faces per parallel task. +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]>, + /// 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, + /// The penalty per length its element sets (before the penalty scale). + stiffness: f64, +} + +/// A contact pair as the solver runs it. +#[derive(Debug, Clone)] +pub struct Interface { + pub name: String, + secondary: Vec, + main: Vec, + /// The secondary nodes (mesh slots), sorted, and per node: its reference area A_A, its + /// penalty ε_A, the thinnest element under it and the offset taken out of its weighted gap + /// (initial overclosure). + pub nodes: Vec, + area: Vec, + penalty: Vec, + thickness: Vec, + offset: Vec, + /// Small sliding (steps without nonlinear geometry): the face pairs and their integrals, + /// fixed in the reference configuration; `None`: finite sliding, searched and integrated + /// on the current configuration at every evaluation. + fixed: Option>>, +} + +/// Whether a secondary node presses on the main surface. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)] +#[serde(rename_all = "snake_case")] +pub enum Status { + Open, + Closed, +} + +impl ContactResult { + /// The closed nodes. + pub fn closed(&self) -> usize { + self.status.iter().filter(|s| **s == Status::Closed).count() + } +} + +/// One evaluation of a contact: per secondary node (in [`Interface::nodes`] order). +#[derive(Debug, Clone, PartialEq, Default)] +pub struct State { + /// ĝ_A, `None` where no main face lies over the node. + pub gap: Vec>, + pub pressure: Vec, + /// ∫ Φ_A over the overlaps (the area the pressure acts on). + pub support: Vec, + /// The force the contact puts on each secondary node, and the node's unit outward normal + /// (its faces' area vectors summed). + pub nodal_force: Vec, + pub normal: Vec, +} + +impl State { + /// A hash of which nodes are closed (for the chattering check). + pub fn active_set(&self) -> u64 { + // FNV-1a over the closed nodes' indices: deterministic on every platform. + let mut h: u64 = 0xcbf2_9ce4_8422_2325; + for (i, _) in self.pressure.iter().enumerate().filter(|(_, p)| **p > 0.0) { + for b in (i as u64).to_le_bytes() { + h ^= u64::from(b); + h = h.wrapping_mul(0x0100_0000_01b3); + } + } + h + } +} + +/// A contact's result at the end of a step, per secondary node. +#[derive(Debug, Clone, PartialEq, Serialize)] +pub struct ContactResult { + pub name: String, + /// The secondary surface's node labels, sorted. + pub nodes: Vec, + /// The contact pressure: the normal contact force on the node over the area it carries + /// (∫ Φ_A over the overlaps), so that Σ p_A ∫Φ_A is the force across the contact. + pub pressure: Vec, + /// The mortar multipliers (the coefficients of the pressure field, Φ_A-interpolated). + pub multiplier: Vec, + /// The normalised gap ĝ (negative: penetration), `None` where no main face is over the node. + pub gap: Vec>, + pub status: Vec, + /// The force the contact puts on the secondary surface, and its normal part Σ p_A ∫Φ_A. + pub force: Vec3, + pub normal_force: f64, + /// The deepest penetration of a closed node. + pub penetration: f64, + /// Augmentations over the step. + pub augmentations: usize, +} + +/// The contribution of one face pair to the forces and the tangent. +struct PairTerms { + /// Mesh slots of the pair's corners (secondary first) and their degrees of freedom. + nodes: Vec, + dofs: Vec<[usize; 3]>, + /// Secondary node indices (in [`Interface::nodes`]) of the secondary corners. + corners: Vec, + /// ∂g̃_a/∂x per secondary corner (24 variables; corner c of the pair at 3c, main corners from 12). + gradient: Vec<[f64; VARS]>, + /// Σ_a p_a ∂²g̃_a, packed (tangent evaluations only). + hessian: Option>, +} + +/// The tangent blocks of a contact: per block its mesh slots and the dense matrix over their +/// translations (row-major, node by node). +pub type Blocks = Vec<(Vec, Vec)>; + +fn position(mesh: &Mesh, u: &[f64], node: usize, dofs: &[usize; 3]) -> Vec3 { + let x = mesh.xyz.get(node).copied().unwrap_or_default(); + [0, 1, 2].map(|k| x[k] + u.get(dofs[k]).copied().unwrap_or(0.0)) +} + +/// The element's volume (hex8, tet4) for the thickness under a face. +fn volume(kind: ElementKind, x: &[Vec3]) -> Option { + match kind { + ElementKind::Hex8 => <[Vec3; 8]>::try_from(x).ok().and_then(|a| fieldcraft_fem::uniform::Hex8::hex8(&a).ok()).map(|k| k.volume), + ElementKind::Tet4 => <[Vec3; 4]>::try_from(x).ok().and_then(|a| fieldcraft_fem::uniform::Tet4::tet4(&a).ok()).map(|k| k.volume), + _ => None, + } +} + +fn size(x: &[Vec3]) -> f64 { + let b = bounds(x, 0.0); + (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> { + 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 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", + el.kind.name(), + f.element + ))); + } + let nodes: Vec = labels.iter().map(|l| mesh.node(*l).ok_or_else(|| at(format!("node {l} is not analysed")))).collect::>()?; + let dofs = nodes + .iter() + .map(|n| translations(mesh, *n).and_then(|d| <[usize; 3]>::try_from(d).map_err(|_| Error::Internal("solid node translations".into())))) + .collect::>>()?; + let x: Vec = nodes.iter().map(|n| mesh.xyz.get(*n).copied().unwrap_or_default()).collect(); + let area = 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); + if !(thickness > 0.0 && thickness.is_finite()) { + return Err(at("the face has no size".into())); + } + let shear = (el.d[3][3] + el.d[4][4] + el.d[5][5]) / 3.0; + let mut owner = el.nodes.clone(); + owner.sort_unstable(); + out.push(Face { nodes, dofs, owner, depth: DEPTH_FRACTION * thickness, thickness, stiffness: PENALTY_FACTOR * 3.0 * shear / thickness }); + } + Ok(out) +} + +/// A face pair of small-sliding contact: its mortar integrals in the reference configuration, +/// D_ai = ∫ Φ_a N_i (secondary corners i) and M_aj = ∫ Φ_a N_j (main corners j), and the +/// secondary face's normal there, so that g̃_a = n · (Σ_j M_aj x_j − Σ_i D_ai x_i). +#[derive(Debug, Clone)] +struct Fixed { + main: usize, + d: [[f64; 4]; 4], + m: [[f64; 4]; 4], + normal: Vec3, + support: [f64; 4], +} + +impl Fixed { + /// ∂g̃_a/∂x over the pair's variables (secondary corners from 0, main corners from 12). + fn gradient(&self, a: usize) -> [f64; VARS] { + let mut g = [0.0; VARS]; + let (Some(d), Some(m)) = (self.d.get(a), self.m.get(a)) else { return g }; + for k in 0..3 { + for i in 0..4 { + g[3 * i + k] = -d[i] * self.normal[k]; + g[12 + 3 * i + k] = m[i] * self.normal[k]; + } + } + g + } +} + +/// What a secondary face finds against the main surface: per main face, ∫ Φ_a g and ∫ Φ_a for +/// each corner a. +struct Found { + main: usize, + gap: [f64; 4], + support: [f64; 4], +} + +/// The weighted gaps g̃_a of one face pair (corners `xs`, `xm`) over any scalar. +fn weighted_gaps(xs: &[[T; 3]], xm: &[[T; 3]]) -> Option<[T; 4]> { + let mut g = [T::cst(0.0); 4]; + integrate(xs, xm, |p| { + let w = p.weight * p.gap(); + for a in 0..4 { + g[a] = g[a] + w * p.secondary[a]; + } + })?; + Some(g) +} + +/// The variables of a face pair: secondary corners at 0, main corners from 12. +fn variables(xs: &[Vec3], xm: &[Vec3], var: impl Fn(f64, usize) -> T) -> (Vec<[T; 3]>, Vec<[T; 3]>) { + let s = xs.iter().enumerate().map(|(i, x)| [0, 1, 2].map(|k| var(x[k], 3 * i + k))).collect(); + let m = xm.iter().enumerate().map(|(j, x)| [0, 1, 2].map(|k| var(x[k], 12 + 3 * j + k))).collect(); + (s, m) +} + +impl Interface { + /// The contact pairs of a model, ready for the nonlinear solver (initial overclosures + /// adjusted or refused, and reported in `warnings`). + pub fn all(model: &Model, resolved: &Resolved, mesh: &Mesh, finite: bool, warnings: &mut Vec) -> Result> { + let mut out = Vec::new(); + let mut total = 0usize; + for (name, contact) in model.contacts() { + let mut interface = Interface::new(model, resolved, mesh, name, contact)?; + if !finite { + interface.fixed = Some(interface.fixed_pairs(mesh)?); + } + total = total.saturating_add(interface.secondary.len() + interface.main.len()); + if total > MAX_CONTACT_FACES { + return Err(Error::Limit(format!("the contacts have more than {MAX_CONTACT_FACES} faces"))); + } + out.push(interface.overclosure(mesh, contact.initial_penetration, warnings)?); + } + Ok(out) + } + + fn new(model: &Model, resolved: &Resolved, mesh: &Mesh, name: &str, contact: &Contact) -> Result { + let ContactSurfaces::Pair { main, secondary } = &contact.surfaces else { + return Err(Error::NotSupported(format!( + "general (single-surface) contact `{name}` in implicit steps; define contact pairs (interaction.contact.pair), or solve it in explicit dynamics" + ))); + }; + if !contact.friction.is_frictionless() { + return Err(Error::NotSupported(format!("friction (contact `{name}`) in implicit steps"))); + } + let secondary = faces(model, resolved, mesh, name, secondary)?; + let main = faces(model, resolved, mesh, name, main)?; + let mut nodes: Vec = secondary.iter().flat_map(|f| f.nodes.iter().copied()).collect(); + nodes.sort_unstable(); + nodes.dedup(); + let (mut area, mut penalty, mut thickness) = (vec![0.0; nodes.len()], vec![0.0f64; nodes.len()], vec![f64::INFINITY; nodes.len()]); + for f in &secondary { + let x: Vec = f.nodes.iter().map(|n| mesh.xyz.get(*n).copied().unwrap_or_default()).collect(); + // ∫ Φ_a over the face itself (the face against its own mirror image). + let mirror: Vec = x.iter().rev().copied().collect(); + let mut own = [0.0; 4]; + integrate(&x, &mirror, |p| { + for a in 0..4 { + own[a] += p.weight * p.secondary[a]; + } + }); + for (n, share) in f.nodes.iter().zip(own) { + if let Ok(i) = nodes.binary_search(n) { + if let Some(a) = area.get_mut(i) { + *a += share; + } + if let Some(e) = penalty.get_mut(i) { + *e = e.max(contact.penalty_scale * f.stiffness); + } + if let Some(h) = thickness.get_mut(i) { + *h = h.min(f.thickness); + } + } + } + } + if let Some(i) = area.iter().position(|a| !(*a > 0.0 && a.is_finite())) { + let label = nodes.get(i).and_then(|n| mesh.labels.get(*n)).copied().unwrap_or(0); + return Err(Error::Invalid(format!("contact `{name}`: the secondary faces at node {label} have no area"))); + } + let offset = vec![0.0; nodes.len()]; + Ok(Interface { name: name.to_string(), secondary, main, nodes, area, penalty, thickness, offset, fixed: None }) + } + + /// Adjust (offset) or refuse the nodes that start behind the main surface. + fn overclosure(mut self, mesh: &Mesh, mode: InitialPenetration, warnings: &mut Vec) -> Result { + let zero = vec![0.0; mesh.dofs()]; + let (gaps, _) = self.gaps(mesh, &zero)?; + let scale = mesh.xyz.iter().flatten().fold(1.0f64, |m, c| m.max(c.abs())); + let (mut count, mut deepest) = (0usize, (0.0f64, 0u32)); + for (i, g) in gaps.iter().enumerate() { + let (Some(g), Some(a)) = (g, self.area.get(i)) else { continue }; + let depth = -g / a; + // Coordinates are exact to a few units in the last place: a node a hair behind a + // face it was meshed on is touching. + if depth <= 1e-9 * scale { + continue; + } + let label = self.nodes.get(i).and_then(|n| mesh.labels.get(*n)).copied().unwrap_or(0); + if mode == InitialPenetration::Refuse { + return Err(Error::Invalid(format!( + "contact `{}`: node {label} starts {depth:.4e} inside the main surface (move it out, or let the contact adjust it: initial_penetration \"adjust\")", + self.name + ))); + } + if let Some(o) = self.offset.get_mut(i) { + *o = *g; + } + count += 1; + if depth > deepest.0 { + deepest = (depth, label); + } + } + if count > 0 { + warnings.push(format!( + "contact `{}`: {count} node(s) start inside the main surface (deepest node {}, {:.4e}); their initial overclosure is taken out of the gap", + self.name, deepest.1, deepest.0 + )); + } + Ok(self) + } + + fn corners(mesh: &Mesh, u: &[f64], f: &Face) -> Vec { + f.nodes.iter().zip(&f.dofs).map(|(n, d)| position(mesh, u, *n, d)).collect() + } + + /// The face pairs: per secondary face, the main faces it faces and their integrals. + fn pairs(&self, mesh: &Mesh, u: &[f64]) -> Result>> { + let Some(fixed) = &self.fixed else { return self.search(mesh, u, 1.0, |_, _, found| found) }; + Ok(parallel::map_grained(self.secondary.len(), GRAIN, |s| { + let (Some(f), Some(list)) = (self.secondary.get(s), fixed.get(s)) else { return Vec::new() }; + let xs = Interface::corners(mesh, u, f); + list.iter() + .map(|p| { + let xm = self.main.get(p.main).map(|m| Interface::corners(mesh, u, m)).unwrap_or_default(); + let gap = std::array::from_fn(|a| { + let (Some(d), Some(m)) = (p.d.get(a), p.m.get(a)) else { return 0.0 }; + let mut v = [0.0; 3]; + for (x, c) in xm.iter().zip(m) { + (0..3).for_each(|k| v[k] += c * x[k]); + } + for (x, c) in xs.iter().zip(d) { + (0..3).for_each(|k| v[k] -= c * x[k]); + } + dot(v, p.normal) + }); + Found { main: p.main, gap, support: p.support } + }) + .collect() + })) + } + + /// The small-sliding pairs: found and integrated once, in the reference configuration. + fn fixed_pairs(&self, mesh: &Mesh) -> Result>> { + let zero = vec![0.0; mesh.dofs()]; + // Faces within a main element's thickness may touch later. + self.search(mesh, &zero, 2.0, |xs, xm, found| { + found + .into_iter() + .filter_map(|f| { + let xm = xm.get(f.main)?; + let (mut d, mut m, mut normal) = ([[0.0; 4]; 4], [[0.0; 4]; 4], [0.0; 3]); + integrate(xs, xm, |p| { + normal = p.normal; + for a in 0..4 { + for i in 0..4 { + d[a][i] += p.weight * p.secondary[a] * p.secondary[i]; + m[a][i] += p.weight * p.secondary[a] * p.main[i]; + } + } + })?; + Some(Fixed { main: f.main, d, m, normal, support: f.support }) + }) + .collect() + }) + } + + /// The global and local search at displacements `u`: per secondary face (corners `xs`), the + /// main faces (corners `xm`, all of them) within `reach` times their depth limit that it + /// faces, with their integrals, passed through `keep`. + fn search(&self, mesh: &Mesh, u: &[f64], reach: f64, keep: impl Fn(&[Vec3], &[Vec], Vec) -> T + Sync) -> Result> { + let main_x: Vec> = self.main.iter().map(|f| Interface::corners(mesh, u, f)).collect(); + let boxes: Vec = self.main.iter().zip(&main_x).map(|(f, x)| bounds(x, reach * f.depth)).collect(); + let grid = Grid::new(&boxes, 4 * boxes.len() + 64); + let main_normals: Vec> = main_x.iter().map(|x| area_vector(x).filter(|a| norm(*a) > 0.0)).collect(); + Ok(parallel::map_grained(self.secondary.len(), GRAIN, |i| { + let mut found = Vec::new(); + let Some(f) = self.secondary.get(i) else { return keep(&[], &main_x, found) }; + let xs = Interface::corners(mesh, u, f); + let Some(ns) = area_vector(&xs).filter(|a| norm(*a) > 0.0) else { return keep(&xs, &main_x, found) }; + let mut candidates = Vec::new(); + grid.candidates_in(&bounds(&xs, 0.0), &mut candidates); + for &c in &candidates { + let c = c as usize; + let (Some(m), Some(xm), Some(Some(nm))) = (self.main.get(c), main_x.get(c), main_normals.get(c)) else { continue }; + if m.nodes.iter().any(|n| f.owner.binary_search(n).is_ok()) || f.nodes.iter().any(|n| m.owner.binary_search(n).is_ok()) { + continue; + } + if dot(ns, *nm) >= -FACING * norm(ns) * norm(*nm) { + continue; + } + let (mut gap, mut support, mut deepest, mut any) = ([0.0; 4], [0.0; 4], f64::INFINITY, false); + integrate(&xs, xm, |p| { + any = true; + let g = p.gap(); + deepest = deepest.min(g); + for a in 0..4 { + gap[a] += p.weight * g * p.secondary[a]; + support[a] += p.weight * p.secondary[a]; + } + }); + if any && deepest >= -m.depth && gap.iter().chain(&support).all(|v| v.is_finite()) { + found.push(Found { main: c, gap, support }); + } + } + keep(&xs, &main_x, found) + })) + } + + /// The weighted gaps g̃_A less the offsets (`None` where no main face lies over the node), and + /// the supports ∫ Φ_A, at displacements `u`. + fn gaps(&self, mesh: &Mesh, u: &[f64]) -> Result<(Vec>, Vec)> { + let pairs = self.pairs(mesh, u)?; + let (gaps, support, _) = self.sums(&pairs); + Ok((gaps, support)) + } + + /// Node sums of the pairs' integrals, in face order. + fn sums(&self, pairs: &[Vec]) -> (Vec>, Vec, Vec>) { + let n = self.nodes.len(); + let (mut gap, mut support) = (vec![None; n], vec![0.0; n]); + let mut index = Vec::with_capacity(self.secondary.len()); + for (f, found) in self.secondary.iter().zip(pairs) { + let corners: Vec = f.nodes.iter().filter_map(|x| self.nodes.binary_search(x).ok()).collect(); + for p in found { + for (a, i) in corners.iter().enumerate() { + if let (Some(g), Some(s)) = (gap.get_mut(*i), support.get_mut(*i)) { + *g = Some(g.unwrap_or(0.0) + p.gap.get(a).copied().unwrap_or(0.0)); + *s += p.support.get(a).copied().unwrap_or(0.0); + } + } + } + index.push(corners); + } + for (g, o) in gap.iter_mut().zip(&self.offset) { + if let Some(g) = g { + *g -= o; + } + } + (gap, support, index) + } + + /// The contact's forces at displacements `u` (added to `internal`), its tangent blocks when + /// `tangent`, and its state, with the multipliers `multipliers` (per node). + pub fn evaluate(&self, mesh: &Mesh, u: &[f64], multipliers: &[f64], tangent: bool, internal: &mut [f64], blocks: &mut Blocks) -> Result { + let pairs = self.pairs(mesh, u)?; + let (weighted, support, corner_index) = self.sums(&pairs); + let n = self.nodes.len(); + let mut state = State { gap: vec![None; n], pressure: vec![0.0; n], support, nodal_force: vec![[0.0; 3]; n], normal: vec![[0.0; 3]; n] }; + for (f, corners) in self.secondary.iter().zip(&corner_index) { + let Some(a) = area_vector(&Interface::corners(mesh, u, f)) else { continue }; + for i in corners { + if let Some(v) = state.normal.get_mut(*i) { + (0..3).for_each(|k| v[k] += a[k]); + } + } + } + for v in &mut state.normal { + *v = unit(*v).unwrap_or(*v); + } + for i in 0..self.nodes.len() { + let (Some(Some(g)), Some(a), Some(e)) = (weighted.get(i), self.area.get(i), self.penalty.get(i)) else { continue }; + let gap = g / a; + let lambda = multipliers.get(i).copied().unwrap_or(0.0); + if let (Some(slot), Some(p)) = (state.gap.get_mut(i), state.pressure.get_mut(i)) { + *slot = Some(gap); + *p = (lambda - e * gap).max(0.0); + } + } + // The pairs with a closed secondary corner, linearised. + let mut active: Vec<(usize, usize)> = Vec::new(); + for (s, found) in pairs.iter().enumerate() { + let closed = corner_index.get(s).is_some_and(|c| c.iter().any(|i| state.pressure.get(*i).is_some_and(|p| *p > 0.0))); + if closed { + active.extend((0..found.len()).map(|k| (s, k))); + } + } + let terms = parallel::map_grained(active.len(), GRAIN, |t| -> Result> { + let Some(&(s, k)) = active.get(t) else { return Ok(None) }; + let (Some(f), Some(found), Some(corners)) = (self.secondary.get(s), pairs.get(s).and_then(|p| p.get(k)), corner_index.get(s)) else { + return Ok(None); + }; + let Some(m) = self.main.get(found.main) else { return Ok(None) }; + let (xs, xm) = (Interface::corners(mesh, u, f), Interface::corners(mesh, u, m)); + let pressure: Vec = corners.iter().map(|i| state.pressure.get(*i).copied().unwrap_or(0.0)).collect(); + let fixed = self.fixed.as_ref().and_then(|f| f.get(s)).and_then(|l| l.get(k)); + let (gradient, hessian) = if let Some(fixed) = fixed { + // Small sliding: linear in the positions. + ((0..corners.len()).map(|a| fixed.gradient(a)).collect(), None) + } else if tangent { + let (vs, vm) = variables(&xs, &xm, Hyper::var); + let Some(g) = weighted_gaps(&vs, &vm) else { return Ok(None) }; + let mut h = vec![0.0; fieldcraft_fem::dual::PACKED]; + for (ga, p) in g.iter().zip(&pressure) { + h.iter_mut().zip(&ga.h).for_each(|(x, y)| *x += p * y); + } + (g.iter().take(corners.len()).map(|ga| ga.g).collect::>(), Some(h)) + } else { + let (vs, vm) = variables(&xs, &xm, Dual::var); + let Some(g) = weighted_gaps(&vs, &vm) else { return Ok(None) }; + (g.iter().take(corners.len()).map(|ga| ga.g).collect(), None) + }; + let nodes = f.nodes.iter().chain(&m.nodes).copied().collect(); + let dofs = f.dofs.iter().chain(&m.dofs).copied().collect(); + Ok(Some(PairTerms { nodes, dofs, corners: corners.clone(), gradient, hessian })) + }); + // Forces −p_A ∂g̃_A/∂x, and per closed node its gradient over the nodes it reaches. + let mut node_gradients: BTreeMap> = BTreeMap::new(); + for t in terms { + let Some(t) = t? else { continue }; + for (a, i) in t.corners.iter().enumerate() { + let p = state.pressure.get(*i).copied().unwrap_or(0.0); + let Some(g) = t.gradient.get(a) else { continue }; + for (c, (node, dofs)) in t.nodes.iter().zip(&t.dofs).enumerate() { + let var = if c < t.corners.len() { 3 * c } else { 12 + 3 * (c - t.corners.len()) }; + let gv = [0, 1, 2].map(|k| g.get(var + k).copied().unwrap_or(0.0)); + let own = t.corners.get(c).and_then(|j| state.nodal_force.get_mut(*j)); + let mut on = [0.0; 3]; + for k in 0..3 { + let f = -p * gv[k]; + if let Some(x) = internal.get_mut(dofs[k]) { + *x += f; + } + on[k] = -f; + } + if let Some(own) = own { + (0..3).for_each(|k| own[k] += on[k]); + } + if p > 0.0 && tangent { + let entry = node_gradients.entry(*i).or_default().entry(*node).or_insert([0.0; 3]); + for k in 0..3 { + entry[k] += gv[k]; + } + } + } + } + if let Some(h) = &t.hessian { + // −Σ_a p_a ∂²g̃_a over the pair's nodes. + let count = t.nodes.len(); + let s = t.corners.len(); + let var = |c: usize, k: usize| if c < s { 3 * c + k } else { 12 + 3 * (c - s) + k }; + let width = 3 * count; + let mut data = vec![0.0; width * width]; + for ci in 0..count { + for ki in 0..3 { + for cj in 0..count { + for kj in 0..3 { + let at = fieldcraft_fem::dual::packed(var(ci, ki), var(cj, kj)); + if let Some(x) = data.get_mut((3 * ci + ki) * width + 3 * cj + kj) { + *x = -h.get(at).copied().unwrap_or(0.0); + } + } + } + } + } + blocks.push((t.nodes.clone(), data)); + } + } + // (ε_A / A_A) ∂g̃_A ⊗ ∂g̃_A per closed node. + for (i, gradient) in node_gradients { + let (Some(e), Some(a)) = (self.penalty.get(i), self.area.get(i)) else { continue }; + let c = e / a; + let nodes: Vec = gradient.keys().copied().collect(); + let g: Vec = gradient.values().flatten().copied().collect(); + let width = g.len(); + let mut data = vec![0.0; width * width]; + for (r, gr) in g.iter().enumerate() { + for (s, gs) in g.iter().enumerate() { + if let Some(x) = data.get_mut(r * width + s) { + *x = c * gr * gs; + } + } + } + blocks.push((nodes, data)); + } + Ok(state) + } + + /// The result of a converged state. + pub fn result(&self, mesh: &Mesh, state: &State, augmentations: usize) -> ContactResult { + let status: Vec = state.pressure.iter().map(|p| if *p > 0.0 { Status::Closed } else { Status::Open }).collect(); + let penetration = state.gap.iter().zip(&status).filter(|(_, s)| **s == Status::Closed).filter_map(|(g, _)| *g).fold(0.0f64, |m, g| m.max(-g)); + // The force on the secondary nodes pushes them back along −n. + let normal: Vec = state.nodal_force.iter().zip(&state.normal).map(|(f, n)| -dot(*f, *n)).collect(); + let pressure = normal.iter().zip(&state.support).map(|(f, a)| if *a > 0.0 { f / a } else { 0.0 }).collect(); + let force = state.nodal_force.iter().fold([0.0; 3], |s, f| [s[0] + f[0], s[1] + f[1], s[2] + f[2]]); + ContactResult { + name: self.name.clone(), + nodes: self.nodes.iter().map(|n| mesh.labels.get(*n).copied().unwrap_or(0)).collect(), + pressure, + multiplier: state.pressure.clone(), + gap: state.gap.clone(), + status, + force, + normal_force: normal.iter().sum(), + penetration, + augmentations, + } + } + + /// How far the augmentation is from converged at `state`, computed with `multipliers`: the + /// largest gap of a closed node, or the gap equivalent λ_A / ε_A of the multiplier an open + /// node keeps, over the thickness of the element under the node. + pub fn augmentation_error(&self, multipliers: &[f64], state: &State) -> f64 { + let mut worst = 0.0f64; + for i in 0..self.nodes.len() { + let (Some(p), Some(e), Some(h)) = (state.pressure.get(i), self.penalty.get(i), self.thickness.get(i)) else { continue }; + let lambda = multipliers.get(i).copied().unwrap_or(0.0); + let gap = if *p > 0.0 { state.gap.get(i).copied().flatten().unwrap_or(0.0).abs() } else { lambda.max(0.0) / e }; + worst = worst.max(gap / h); + } + worst + } +} diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index 3a777a3..d55a534 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -31,6 +31,7 @@ pub mod analysis; pub mod buckle; pub mod cholesky; pub mod constraints; +pub mod contact; pub mod discrete; pub mod eigen; pub mod heat; diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs index 3f9cc31..fedb1f3 100644 --- a/crates/solver/src/nonlinear.rs +++ b/crates/solver/src/nonlinear.rs @@ -88,6 +88,15 @@ //! strain at the integration points. //! Couplings stay the linear (small-rotation) constraints of the linear solver, with a warning. //! +//! **Contact** ([`crate::contact`]): the model's contact pairs add their mortar forces to the +//! internal forces and their consistent tangent to the stiffness, whose pattern takes the blocks +//! of the closed nodes' face pairs at every evaluation. Each converged increment is solved again +//! with updated multipliers until the gaps are tight (the augmented Lagrangian's Uzawa loop, the +//! increment's `augmentations`, each iteration tagged with its `augmentation`); an increment whose +//! active set chatters is cut back. With contact the force scale of the convergence test is that +//! of the converged states, the loads and the current state, not the largest force of every +//! trial: a trial state can press a body deep into another before the active set settles. +//! //! Element forces and tangents are computed in parallel and added in label order, so results are //! the same bit for bit on any number of threads. @@ -111,8 +120,9 @@ use crate::analysis::{ }; use crate::cholesky::FactorError; use crate::constraints::Reduction; +use crate::contact::{self, ContactResult, Interface}; use crate::ldlt::{Kind, Symmetric}; -use crate::sparse::{Csr, add_block}; +use crate::sparse::{Csr, NodeGraph, add_block, add_entry}; use crate::system::{Mesh, add_element_thermal_loads, loads, temperature_loads, translations}; use crate::{Error, Limits, Result}; @@ -131,6 +141,9 @@ pub struct Iteration { /// tangent is factored by LDLᵀ, [`crate::ldlt`]). #[serde(skip_serializing_if = "is_zero")] pub negative_pivots: usize, + /// Contact: the augmentation it belongs to (0: the increment's first solution). + #[serde(skip_serializing_if = "is_zero")] + pub augmentation: usize, } fn is_zero(n: &usize) -> bool { @@ -155,6 +168,9 @@ pub struct Increment { /// factor it for the next predictor). #[serde(skip_serializing_if = "is_zero")] pub negative_pivots: usize, + /// Contact: the times the increment was solved again with updated multipliers. + #[serde(skip_serializing_if = "is_zero")] + pub augmentations: usize, } /// Where the tangent's inertia changes along an arc-length path. @@ -259,8 +275,12 @@ struct Problem<'a> { /// The step's temperature field at λ = 1, per mesh node (`None` where the field leaves a /// node at its material's reference temperature); `None` without one. field: Option>>, - /// The tangent's pattern, zero. + /// The tangent's pattern, zero, and the mesh's node graph it comes from (contact adds the + /// blocks of its face pairs to it). pattern: Csr, + graph: NodeGraph, + /// The contact pairs ([`crate::contact`]). + contacts: Vec, } /// The internal and external forces at a state, and the tangent when asked for. @@ -271,6 +291,8 @@ struct Evaluation { /// Every element's material history at this state (committed when the increment /// converges). history: Vec>, + /// Each contact's state. + contact: Vec, } /// Why an attempt at an increment failed, when that calls for a cutback rather than an error. @@ -291,6 +313,11 @@ impl From for Failure { } } +/// One hash of the contacts' active sets. +fn active_sets(e: &Evaluation) -> u64 { + e.contact.iter().fold(0u64, |h, s| h.rotate_left(17) ^ s.active_set()) +} + fn norm_inf(v: &[f64]) -> f64 { v.iter().fold(0.0, |m, x| m.max(x.abs())) } @@ -387,11 +414,30 @@ impl Problem<'_> { /// The forces at displacements `u` and load factor λ, and the tangent (symmetric part of the /// load stiffness included) when `tangent`. - fn evaluate(&self, u: &[f64], lambda: f64, tangent: bool, committed: &[Vec]) -> Result { + /// `multipliers` are the contacts' augmented Lagrange multipliers per secondary node (none: + /// zero). + fn evaluate(&self, u: &[f64], lambda: f64, tangent: bool, committed: &[Vec], multipliers: &[Vec]) -> Result { let layout = &self.mesh.layout; let n = u.len(); let mut internal = vec![0.0; n]; - let mut k = tangent.then(|| self.pattern.clone()); + // Contact first: its face pairs extend the tangent's pattern. + let mut blocks = contact::Blocks::new(); + let mut states = Vec::with_capacity(self.contacts.len()); + for (i, c) in self.contacts.iter().enumerate() { + states.push(c.evaluate(self.mesh, u, multipliers.get(i).map(Vec::as_slice).unwrap_or_default(), tangent, &mut internal, &mut blocks)?); + } + let mut k = tangent.then(|| { + if blocks.is_empty() { + self.pattern.clone() + } else { + self.graph.with_blocks(blocks.iter().map(|(nodes, _)| nodes.as_slice())).zero_matrix(layout) + } + }); + if let Some(k) = k.as_mut() { + for (nodes, data) in &blocks { + add_block(k, layout, nodes, self.mesh.translation_kinds(), data)?; + } + } let mut history = Vec::with_capacity(self.mesh.elements.len()); for round in parallel::blocks(self.mesh.elements.len(), ROUND) { let start = round.start; @@ -420,9 +466,18 @@ impl Problem<'_> { for (x, v) in internal.iter_mut().zip(&ku) { *x += v; } - // Both are on the mesh's pattern. + // On the mesh's pattern, which contact may have extended. if let Some(k) = k.as_mut() { - k.val.iter_mut().zip(&s.val).for_each(|(x, v)| *x += v); + if k.ptr == s.ptr && k.idx == s.idx { + k.val.iter_mut().zip(&s.val).for_each(|(x, v)| *x += v); + } else { + for r in 0..s.n { + let (cols, vals) = s.row(r); + for (c, v) in cols.iter().zip(vals).filter(|(_, v)| **v != 0.0) { + add_entry(k, r, *c, *v)?; + } + } + } } } let mut external: Vec = self.dead.iter().map(|f| lambda * f).collect(); @@ -440,7 +495,7 @@ impl Problem<'_> { if internal.iter().chain(&external).any(|v| !v.is_finite()) { return Err(Error::NotConverged("the forces are not finite".into())); } - Ok(Evaluation { internal, external, tangent: k, history }) + Ok(Evaluation { internal, external, tangent: k, history, contact: states }) } /// The loads at load factor 1 at displacements `u` (dead loads and follower pressures on @@ -512,18 +567,33 @@ struct Run<'a> { path_dependent: bool, /// How the undeformed model's tangent factors (Cholesky), once asked. undeformed: Option>, + /// The contacts' multipliers per secondary node (the last converged increment's, or the + /// current augmentation's), and the augmentations over the step. + multipliers: Vec>, + augmentations: usize, + /// The augmentation being solved in the current increment. + augmentation: usize, + /// The largest force of the converged states (contact's force scale, [`Run::residual`]). + converged_scale: f64, } /// The outcome of a converged increment. struct Converged { q: Vec, evaluation: Evaluation, + /// The work of the increment's first correction (the energy criterion's reference). + energy: Option, } impl Run<'_> { /// The residual over the free equations, T_fᵀ (f_ext − f_int), and the force scale. + /// + /// With contact a trial state can be far from equilibrium (a body pressed deep into another + /// before the active set settles, its contact forces huge), so the scale is not the largest + /// force of every trial but that of the converged states, the loads and the current state. fn residual(&mut self, e: &Evaluation) -> Vec { - self.force_scale = self.force_scale.max(norm_inf(&e.external)).max(norm_inf(&e.internal)); + let current = norm_inf(&e.external).max(norm_inf(&e.internal)); + self.force_scale = if self.problem.contacts.is_empty() { self.force_scale.max(current) } else { self.converged_scale.max(current) }; free_residual(self.problem.red, e) } @@ -541,7 +611,7 @@ impl Run<'_> { } fn state(&self, q: &[f64], lambda: f64, tangent: bool) -> std::result::Result { - Ok(self.problem.evaluate(&self.problem.red.expand_scaled(q, lambda), lambda, tangent, &self.committed)?) + Ok(self.problem.evaluate(&self.problem.red.expand_scaled(q, lambda), lambda, tangent, &self.committed, &self.multipliers)?) } /// How the tangent of the undeformed, unloaded model factors (Cholesky): a first tangent @@ -551,7 +621,7 @@ impl Run<'_> { if let Some(u) = &self.undeformed { return Ok(u.clone()); } - let mut e = self.problem.evaluate(&vec![0.0; self.problem.mesh.dofs()], 0.0, true, &self.committed)?; + let mut e = self.problem.evaluate(&vec![0.0; self.problem.mesh.dofs()], 0.0, true, &self.committed, &[])?; let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; let [kff] = self.problem.red.reduce([&k])?; let u = factor_raw(&kff, self.limits).map(|_| ()); @@ -600,7 +670,7 @@ impl Run<'_> { /// factored (the increment then starts from `q0`). fn tangent_predictor(&mut self, q0: &[f64], from: f64, to: f64) -> std::result::Result, f64)>, Failure> { let red = self.problem.red; - let mut e = self.problem.evaluate(&red.expand_scaled(q0, from), to, true, &self.committed)?; + let mut e = self.problem.evaluate(&red.expand_scaled(q0, from), to, true, &self.committed, &self.multipliers)?; let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; let rf: Vec = free_residual(red, &e).iter().zip(prescribed_force(red, &k, to - from)).map(|(a, b)| a - b).collect(); let mut dq = rf.clone(); @@ -642,11 +712,15 @@ impl Run<'_> { // predictor's when the increment started from one that is still in use. let mut first_energy: Option = predicted_energy.filter(|_| q.as_slice() != q0); let mut history: Vec = Vec::new(); + // Contact: the active sets of the iterations, and how often they went back to the one + // before the last (chattering). + let mut sets = (None, active_sets(&eval)); + let mut flips = 0usize; let max_iterations = usize::try_from(self.settings.max_iterations).unwrap_or(usize::MAX); for _ in 0..max_iterations { // The round-off rule: nothing left to correct. if norm_inf(&r) <= ROUND_OFF * tol.force * self.scale() && self.force_floor.is_some() { - return Ok(Converged { q, evaluation: eval }); + return Ok(Converged { q, evaluation: eval, energy: first_energy }); } if self.iterations >= self.limits.max_newton_iterations { return Err(Failure::Fatal(Error::Limit(format!( @@ -677,6 +751,16 @@ impl Run<'_> { *x += s * d; } eval = next; + if !self.problem.contacts.is_empty() { + let set = active_sets(&eval); + if sets.0 == Some(set) && set != sets.1 { + flips += 1; + if flips >= contact::CHATTER_LIMIT { + return Err(Failure::Cutback("the contact status alternates between two states (chattering)".into())); + } + } + sets = (Some(sets.1), set); + } let r_old = std::mem::replace(&mut r, self.residual(&eval)); let energy = (s * dot(&dq, &r_old)).abs(); let first = *first_energy.get_or_insert(energy); @@ -688,6 +772,7 @@ impl Run<'_> { energy: if first > 0.0 { energy / first } else { 0.0 }, line_search: s, negative_pivots: factor.inertia().negative, + augmentation: self.augmentation, }; iterations.push(it); if !it.force.is_finite() { @@ -695,7 +780,7 @@ impl Run<'_> { } let energy_met = it.energy <= tol.energy || started_converged; if (it.force <= tol.force && it.displacement <= tol.displacement && energy_met) || it.force <= ROUND_OFF * tol.force { - return Ok(Converged { q, evaluation: eval }); + return Ok(Converged { q, evaluation: eval, energy: first_energy }); } history.push(it.force); if let [.., a, b, c] = history[..] @@ -708,6 +793,53 @@ impl Run<'_> { Err(Failure::Cutback(format!("no convergence in {} iterations", self.settings.max_iterations))) } + /// The augmented Lagrangian loop of a converged increment (Uzawa: the contacts' multipliers + /// take the converged pressures, and the increment is solved again from its solution) until + /// the gaps are within [`contact::GAP_TOLERANCE`] ([`Interface::augmentation_error`]), at + /// most [`contact::MAX_AUGMENTATIONS`] times (then with a warning). Without contact it + /// returns the increment as it is. + fn augment( + &mut self, + mut c: Converged, + q0: &[f64], + lambda: f64, + size: f64, + iterations: &mut Vec, + count: &mut usize, + warnings: &mut Vec, + ) -> std::result::Result { + if self.problem.contacts.is_empty() { + return Ok(c); + } + loop { + let error = self + .problem + .contacts + .iter() + .zip(&self.multipliers) + .zip(&c.evaluation.contact) + .fold(0.0f64, |m, ((i, l), s)| m.max(i.augmentation_error(l, s))); + if error <= contact::GAP_TOLERANCE { + return Ok(c); + } + if *count >= contact::MAX_AUGMENTATIONS { + warnings.push(format!( + "contact: at load factor {lambda:.6} a gap is still {error:.1e} of its element's thickness after {} augmentations (the penalty's regularisation)", + contact::MAX_AUGMENTATIONS + )); + return Ok(c); + } + *count += 1; + self.augmentations += 1; + self.augmentation = *count; + self.multipliers = c.evaluation.contact.iter().map(|s| s.pressure.clone()).collect(); + // From the converged state, measured against the increment's start and its first + // correction like the increment itself. + let start = std::mem::take(&mut c.q); + c = self.increment(q0, Some(&start), c.energy, lambda, size, iterations)?; + } + } + /// The line search along `dq` from `q`: the accepted step and the state there (with the /// tangent). fn line_search(&mut self, q: &[f64], dq: &[f64], lambda: f64, g0: f64) -> std::result::Result<(f64, Evaluation), Failure> { @@ -785,6 +917,14 @@ fn check_supported(model: &Model, mesh: &Mesh, step: &Step, settings: &Nonlinear return Err(at(format!("centrifugal loads (load `{name}`), which depend on the deformed position"))); } plasticity(model, step, settings, warnings)?; + if let Some(name) = model.contacts().keys().next() { + if mesh.theory.is_some() { + return Err(at(format!("contact (`{name}`) in 2D models"))); + } + if settings.arc_length.is_some() { + return Err(at(format!("contact (`{name}`) with arc-length control"))); + } + } if settings.nonlinear_geometry && let Some(name) = model.couplings().keys().next() { @@ -960,6 +1100,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R } else { None }; + let contacts = Interface::all(model, &resolved, &mesh, settings.nonlinear_geometry, &mut warnings)?; let element_kinds = mesh.elements.iter().map(|e| Ok(e.kernel()?.dofs().node)).collect::>>()?; let element_dofs = mesh .elements @@ -988,6 +1129,8 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R .collect(), field, pattern: graph.zero_matrix(&mesh.layout), + graph, + contacts, }; // Elements on their linear stiffness take the thermal strain as the equivalent nodal // forces, a dead load proportional to the load factor like the temperatures. @@ -1008,7 +1151,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R .collect::>>()?; // The undeformed state first: an element that is bad before anything moves is an error in // the model, not a reason to cut increments back. - problem.evaluate(&vec![0.0; mesh.dofs()], 0.0, false, &committed)?; + problem.evaluate(&vec![0.0; mesh.dofs()], 0.0, false, &committed, &[])?; let path_dependent = problem.materials.iter().any(MaterialModel::is_path_dependent); let mut run = Run { problem: &problem, @@ -1022,6 +1165,10 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R committed, path_dependent, undeformed: None, + multipliers: problem.contacts.iter().map(|c| vec![0.0; c.nodes.len()]).collect(), + augmentations: 0, + augmentation: 0, + converged_scale: 0.0, }; let mut history = History::default(); @@ -1085,10 +1232,21 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R Err(Failure::Fatal(e)) => return Err(e), } } - match run.increment(&q, predictor.as_deref(), predicted_energy, target_lambda, change.abs(), &mut iterations) { + let multipliers = run.multipliers.clone(); + let mut augmentations = 0; + run.augmentation = 0; + let solved = run + .increment(&q, predictor.as_deref(), predicted_energy, target_lambda, change.abs(), &mut iterations) + .and_then(|c| run.augment(c, &q, target_lambda, change.abs(), &mut iterations, &mut augmentations, &mut warnings)); + if solved.is_err() { + run.multipliers = multipliers; + } + match solved { Ok(c) => { previous = Some((c.q.iter().zip(&q).map(|(a, b)| a - b).collect(), change)); - let count = iterations.len(); + // The increment's own Newton iterations size the next one (augmentations + // solve the same increment again). + let count = iterations.iter().filter(|i| i.augmentation == 0).count(); history.increments.push(Increment { number, load_factor: target_lambda, @@ -1097,11 +1255,13 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R converged: true, cutback: None, negative_pivots: 0, + augmentations, }); (t, lambda) = (target, target_lambda); q = c.q; keep_frame(&mut history, lambda, &red.expand_scaled(&q, lambda), &mut warnings); run.committed.clone_from(&c.evaluation.history); + run.converged_scale = run.converged_scale.max(norm_inf(&c.evaluation.external)).max(norm_inf(&c.evaluation.internal)); last = Some(c.evaluation); fast = if count <= FAST_ITERATIONS { fast + 1 } else { 0 }; if fast >= 2 { @@ -1121,6 +1281,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R converged: false, cutback: Some(reason.clone()), negative_pivots: 0, + augmentations: 0, }); fast = 0; size = this * CUTBACK; @@ -1139,7 +1300,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R let u = red.expand_scaled(&q, lambda); let evaluation = match last { Some(e) => e, - None => problem.evaluate(&u, lambda, false, &run.committed)?, + None => problem.evaluate(&u, lambda, false, &run.committed, &run.multipliers)?, }; if let Some((at, negative)) = history .increments @@ -1160,7 +1321,8 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R "follower pressure: the skew part of its load stiffness ({a:.1e} of it, over the free degrees of freedom) is left out of the tangent; the solution is exact, Newton may converge more slowly" )); } - finish(&mesh, step, &red, &problem, u, lambda, evaluation, &run.committed, history, run.factor_entries, warnings) + let contacts = problem.contacts.iter().zip(&evaluation.contact).map(|(c, s)| c.result(&mesh, s, run.augmentations)).collect(); + finish(&mesh, step, &red, &problem, u, lambda, evaluation, &run.committed, history, run.factor_entries, contacts, warnings) } /// The final state as a static result: displacements, reactions, Cauchy stresses and @@ -1177,6 +1339,7 @@ fn finish( committed: &[Vec], history: History, factor_entries: usize, + contacts: Vec, mut warnings: Vec, ) -> Result { let (displacement, rotation, point_motion) = split(&u, mesh); @@ -1248,7 +1411,7 @@ fn finish( } let stress: Vec<[f64; 6]> = sum.iter().zip(&count).map(|(s, c)| s.map(|v| if *c > 0 { v / f64::from(*c) } else { 0.0 })).collect(); let von_mises = stress.iter().map(von_mises).collect(); - let scale = applied_total.iter().chain(&reaction_total).fold(0.0f64, |m, v| m.max(v.abs())); + let scale = crate::system::balance_scale(&applied_total, &reaction_total, &reaction); let ground = mesh.discrete.ground_force(&u, mesh); let imbalance = (0..3).map(|c| (applied_total[c] + reaction_total[c] + ground[c]).abs()).fold(0.0, f64::max); if scale > 0.0 && imbalance > 1e-6 * scale { @@ -1273,6 +1436,7 @@ fn finish( warnings, temperature: None, nonlinear: Some(history), + contacts, }) } @@ -1403,6 +1567,7 @@ impl Run<'_> { energy: if first > 0.0 { energy / first } else { 0.0 }, line_search: 1.0, negative_pivots: factor.inertia().negative, + augmentation: 0, }); } Err(Failure::Cutback(format!("no convergence in {} iterations", self.settings.max_iterations))) @@ -1491,6 +1656,7 @@ fn arc_length( converged: true, cutback: None, negative_pivots: negative, + augmentations: 0, }); states.push((number, lambda, negative)); last = Some(evaluation); @@ -1516,6 +1682,7 @@ fn arc_length( converged: false, cutback: Some(reason.clone()), negative_pivots: 0, + augmentations: 0, }); fast = 0; size *= CUTBACK; diff --git a/crates/solver/src/sparse.rs b/crates/solver/src/sparse.rs index 7ab54d9..c6bc687 100644 --- a/crates/solver/src/sparse.rs +++ b/crates/solver/src/sparse.rs @@ -105,6 +105,30 @@ impl NodeGraph { NodeGraph { adjacency } } + /// The graph with more blocks, each coupling all its slots (contact: the slots a contact + /// element couples change with the deformation). + pub fn with_blocks<'a>(&self, blocks: impl Iterator) -> NodeGraph { + let mut adjacency = self.adjacency.clone(); + let mut touched = Vec::new(); + for b in blocks { + for &a in b { + if let Some(row) = adjacency.get_mut(a) { + row.extend_from_slice(b); + touched.push(a); + } + } + } + touched.sort_unstable(); + touched.dedup(); + for a in touched { + if let Some(row) = adjacency.get_mut(a) { + row.sort_unstable(); + row.dedup(); + } + } + NodeGraph { adjacency } + } + /// Nonzeros of the dof-level matrix: a dense block per node pair. pub fn entries(&self, layout: &DofLayout) -> usize { self.adjacency.iter().enumerate().map(|(a, row)| layout.range(a).len() * row.iter().map(|b| layout.range(*b).len()).sum::()).sum() diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index b7b2df7..1372fc1 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -426,11 +426,6 @@ impl Mesh { /// Build it from the model and its resolved view, refusing what the solver cannot represent. pub fn new(model: &Model, resolved: &Resolved, max_dofs: usize) -> Result { - if let Some(name) = model.contacts().keys().next() { - return Err(Error::NotSupported(format!( - "contact (`{name}`) in implicit analyses (it comes with M7): a model's contacts apply to all its steps, so solve its explicit dynamic steps, or remove the contacts to solve the implicit ones" - ))); - } let mut elements = Vec::new(); let mut trusses = Vec::new(); let mut structural = Vec::new(); @@ -661,6 +656,24 @@ impl Mesh { } } +/// The scale of the loads-and-reactions balance check: the largest load or reaction, nodal ones +/// too, since reactions that cancel out (a part squeezed between held faces) still set it. +pub(crate) fn balance_scale(applied: &[f64; 3], total: &[f64; 3], reaction: &[[f64; 3]]) -> f64 { + let largest = reaction.iter().flatten().fold(0.0f64, |m, v| m.max(v.abs())); + applied.iter().chain(total).fold(largest, |m, v| m.max(v.abs())) +} + +/// Refuse contact in an analysis that cannot take it (linear static, frequency, buckling): a +/// model's contacts apply to all its steps. +pub fn refuse_contact(model: &Model, analysis: &str) -> Result<()> { + match model.contacts().keys().next() { + Some(name) => Err(Error::NotSupported(format!( + "contact (`{name}`) in {analysis} steps: a model's contacts apply to all its steps, so solve its nonlinear static or explicit dynamic steps, or remove the contacts to solve this one" + ))), + None => Ok(()), + } +} + /// The prescribed value of every degree of freedom (`None` = free) under a step's boundary /// conditions; `values: false` holds them at zero (frequency steps). Rotations of solid nodes do /// not exist and are ignored, like the directions of a reference point no coupling ties diff --git a/crates/solver/tests/contact.rs b/crates/solver/tests/contact.rs new file mode 100644 index 0000000..bc52305 --- /dev/null +++ b/crates/solver/tests/contact.rs @@ -0,0 +1,405 @@ +//! Contact in nonlinear static steps against closed forms: the contact patch test on +//! non-matching meshes, Newton's convergence with contact, and what implicit contact refuses. +//! Run with `--nocapture` to see the iterations. + +use std::collections::BTreeMap; + +use fieldcraft_model::*; +use fieldcraft_solver::contact::Status; +use fieldcraft_solver::{Limits, StaticResult, solve_nonlinear}; + +/// A model built block by block. +#[derive(Default)] +struct Blocks { + model: Model, + xyz: BTreeMap, + next_node: u32, + next_element: u32, +} + +impl Blocks { + /// A structured block of `n` cells over the unit cube mapped by `map`, hex8 or tet4 (six + /// Kuhn tetrahedra per cell); its elements' labels. + fn block(&mut self, kind: ElementKind, n: [usize; 3], map: impl Fn([f64; 3]) -> [f64; 3]) -> Vec { + let first = self.next_node + 1; + let label = |i: [usize; 3]| first + (i[0] + (n[0] + 1) * (i[1] + (n[1] + 1) * i[2])) as u32; + for k in 0..=n[2] { + for j in 0..=n[1] { + for i in 0..=n[0] { + let x = map([i as f64 / n[0] as f64, j as f64 / n[1] as f64, k as f64 / n[2] as f64]); + let l = label([i, j, k]); + self.model.set_node(l, x).unwrap(); + self.xyz.insert(l, x); + self.next_node = self.next_node.max(l); + } + } + } + let mut out = Vec::new(); + for k in 0..n[2] { + for j in 0..n[1] { + for i in 0..n[0] { + let at = |d: [usize; 3]| label([i + d[0], j + d[1], k + d[2]]); + let cells: Vec> = match kind { + ElementKind::Hex8 => vec![ + [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]].iter().map(|d| at(*d)).collect(), + ], + _ => [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]] + .iter() + .map(|perm| { + let mut v = [[0usize; 3]; 4]; + for s in 0..3 { + v[s + 1] = v[s]; + v[s + 1][perm[s]] = 1; + } + let p = v.map(|c| c.map(|x| x as f64)); + let e = |i: usize| [0, 1, 2].map(|k| p[i][k] - p[0][k]); + let (a, b, c) = (e(1), e(2), e(3)); + let vol = + a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]); + if vol < 0.0 { + v.swap(1, 2); + } + v.iter().map(|d| at(*d)).collect() + }) + .collect(), + }; + for nodes in cells { + self.next_element += 1; + self.model.set_element(self.next_element, kind, nodes).unwrap(); + out.push(self.next_element); + } + } + } + } + out + } + + /// The faces of `elements` whose nodes all satisfy `f`. + fn faces(&self, elements: &[u32], f: impl Fn([f64; 3]) -> bool) -> Vec { + let mut out = Vec::new(); + for l in elements { + let e = self.model.element(*l).unwrap(); + for (i, face) in e.kind.faces().iter().enumerate() { + if face.iter().all(|k| f(self.xyz[&e.nodes[*k]])) { + out.push(ElementFace { element: *l, face: i as u8 + 1 }); + } + } + } + out + } + + fn nodes(&self, f: impl Fn([f64; 3]) -> bool) -> Vec { + self.xyz.iter().filter(|(_, x)| f(**x)).map(|(l, _)| *l).collect() + } + + fn node_set(&mut self, name: &str, f: impl Fn([f64; 3]) -> bool) { + let nodes = self.nodes(f); + self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap(); + } + + fn face_set(&mut self, name: &str, elements: &[u32], f: impl Fn([f64; 3]) -> bool) { + let faces = self.faces(elements, f); + assert!(!faces.is_empty(), "no faces for `{name}`"); + self.model.add_set(name, Set { faces, ..Set::default() }).unwrap(); + } + + fn bc(&mut self, name: &str, values: [Option; 3], f: impl Fn([f64; 3]) -> bool) { + self.node_set(name, f); + let [a, b, c] = values; + self.model + .add_boundary_condition(name, BoundaryCondition::Displacement { target: Target::Set(name.into()), values: [a, b, c, None, None, None] }) + .unwrap(); + self.model.set_boundary_condition_active("Step", name, true).unwrap(); + } + + /// One elastic material and section over every element, and the nonlinear step. + fn define(&mut self, young: f64, poisson: f64, settings: NonlinearStatic) { + let all: Vec = self.model.elements().keys().copied().collect(); + self.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap(); + self.model.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young, poisson }), ..Material::default() }).unwrap(); + self.model + .add_section( + "solid", + Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() }, + ) + .unwrap(); + self.model.add_step(Step::new("Step", Procedure::NonlinearStatic(settings)), None).unwrap(); + } + + fn pair(&mut self, name: &str, main: &str, secondary: &str) { + let c = Contact::new(ContactSurfaces::Pair { main: main.into(), secondary: secondary.into() }); + self.model.add_contact(name, c).unwrap(); + } +} + +fn near(a: f64, b: f64) -> bool { + (a - b).abs() < 1e-9 +} + +fn report(r: &StaticResult) { + let h = r.nonlinear.as_ref().unwrap(); + for inc in &h.increments { + let forces: Vec = inc.iterations.iter().map(|i| format!("{:.1e}", i.force)).collect(); + eprintln!("increment {} λ {:.3} aug {} [{}]", inc.number, inc.load_factor, inc.augmentations, forces.join(" ")); + } +} + +/// Two blocks pressed together across non-matching meshes (`lower` × `upper` cells, each of +/// `lower_kind` / `upper_kind`), between frictionless symmetry planes: a uniaxial stress. +fn patch(lower_kind: ElementKind, upper_kind: ElementKind, geometry: bool, distorted: bool) -> (Blocks, StaticResult) { + let mut b = Blocks::default(); + // A graded lower block, its cells boxes of different sizes, or distorted (its faces stay the + // planes of the cube either way). + 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 upper = [2, 2, 2]; + let up = b.block(upper_kind, upper, |p| [p[0], p[1], 1.0 + p[2]]); + b.define( + 1000.0, + 0.3, + NonlinearStatic { + nonlinear_geometry: geometry, + initial_increment: 0.5, + // Tight for small sliding (linear: exact); finite sliding converges to the default ones. + tolerances: if geometry { Tolerances::default() } else { Tolerances { force: 1e-10, displacement: 1e-8, energy: 1e-16 } }, + ..NonlinearStatic::default() + }, + ); + 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.pair("interface", "lower_top", "upper_bottom"); + b.bc("x0", [Some(0.0), None, None], |x| near(x[0], 0.0)); + b.bc("y0", [None, Some(0.0), None], |x| near(x[1], 0.0)); + b.bc("z0", [None, None, Some(0.0)], |x| near(x[2], 0.0)); + b.bc("push", [None, None, Some(-0.02)], |x| near(x[2], 2.0)); + let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap(); + report(&r); + (b, r) +} + +/// The contact patch test (Taylor & Papadopoulos 1991; Crisfield 2000): a uniform pressure +/// crosses a non-matching interface exactly, so every element of both blocks has the same +/// uniaxial stress, with nothing across. +#[test] +fn the_mortar_interface_passes_the_contact_patch_test() { + for (lower, upper, name) in [ + (ElementKind::Hex8, ElementKind::Hex8, "hex8 on hex8"), + (ElementKind::Hex8, ElementKind::Tet4, "tet4 on hex8"), + (ElementKind::Tet4, ElementKind::Hex8, "hex8 on tet4"), + ] { + for geometry in [false, true] { + // Small sliding is linear and exact; finite sliding is exact up to the gaps the + // augmentation leaves (their geometric terms, GAP_TOLERANCE) and the convergence. + let tol = if geometry { 1e-4 } else { 1e-9 }; + let (b, r) = patch(lower, upper, geometry, false); + let szz: Vec = r.elements.iter().flat_map(|e| e.nodal.iter().map(|s| s[2])).collect(); + let mean = szz.iter().sum::() / szz.len() as f64; + assert!(mean < 0.0, "{name}: {mean}"); + for (e, s) in r.elements.iter().flat_map(|e| e.nodal.iter().map(move |s| (e.label, s))) { + assert!((s[2] - mean).abs() < tol * mean.abs(), "{name} ({geometry}) element {e}: σzz {} vs {mean}", s[2]); + assert!(s[0].abs() < tol * mean.abs() && s[1].abs() < tol * mean.abs(), "{name} ({geometry}) element {e}: {s:?}"); + } + let c = &r.contacts[0]; + // The pressure is uniform, equal to the stress, at every secondary node. + assert!(c.pressure.iter().all(|p| (p + mean).abs() < tol * mean.abs()), "{name}: {:?} vs {mean}", c.pressure); + // The contact force on the upper block balances the reaction at its top. + let top: f64 = b.nodes(|x| near(x[2], 2.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()][2]).sum(); + assert!((c.force[2] + top).abs() < tol * top.abs(), "{name}: {:?} vs {top}", c.force); + assert!(c.force[0].abs() < tol * top.abs() && c.force[1].abs() < tol * top.abs(), "{name}: {:?}", c.force); + assert!((c.normal_force + top).abs() < tol * top.abs(), "{name}: {} vs {top}", c.normal_force); + // Augmented to a gap tolerance far below the penetration a penalty alone leaves. + assert!(c.penetration < 1e-5, "{name}: {}", c.penetration); + } + } +} + +/// Hertz contact of a cylinder on a flat (plane strain): a quarter model of an indenter whose +/// bottom is a circle of radius `r` on a block, symmetric about x = 0, one element thick in z +/// with z held (plane strain). Returns the model and the result. +fn hertz_cylinder(finite: bool, push: f64) -> (Blocks, StaticResult) { + let r = 10.0; + let mut b = Blocks::default(); + let grade = |s: f64| s.powf(1.6); + let foundation = b.block(ElementKind::Hex8, [24, 14, 1], |p| [5.0 * grade(p[0]), -5.0 * grade(1.0 - p[1]), 0.2 * p[2]]); + let indenter = b.block(ElementKind::Hex8, [20, 10, 1], |p| { + let x = 2.5 * grade(p[0]); + [x, r - (r * r - x * x).sqrt() + 2.5 * grade(p[1]), 0.2 * p[2]] + }); + b.define(1000.0, 0.3, NonlinearStatic { nonlinear_geometry: finite, initial_increment: 0.25, ..NonlinearStatic::default() }); + b.face_set("flat", &foundation, |x| near(x[1], 0.0)); + let bottom = |x: [f64; 3]| (x[1] - (r - (r * r - x[0] * x[0]).sqrt())).abs() < 1e-9; + b.face_set("round", &indenter, bottom); + b.pair("hertz", "flat", "round"); + b.bc("symmetry", [Some(0.0), None, None], |x| near(x[0], 0.0)); + b.bc("plane", [None, None, Some(0.0)], |_| true); + b.bc("base", [None, Some(0.0), None], |x| near(x[1], -5.0)); + let top = 2.5 + r - (r * r - 2.5 * 2.5f64).sqrt(); + b.bc("push", [None, Some(-push), None], move |x| { + x[1] > 2.5 && (x[1] - (r - (r * r - x[0] * x[0]).sqrt()) - 2.5).abs() < 1e-9 && x[1] <= top + 1e-9 + }); + let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap(); + report(&r); + (b, r) +} + +/// Hertz's line contact (H. Hertz, *Über die Berührung fester elastischer Körper*, J. reine +/// angew. Math. 92:156–171, 1882; K. L. Johnson, *Contact Mechanics*, Cambridge 1985, §4.3): +/// for a load P per length, b = √(4 P R / (π E*)) and p₀ = 2 P / (π b), E* = E / (2 (1 − ν²)) +/// for two bodies of one material. +#[test] +fn hertz_cylinder_on_a_flat() { + for finite in [false, true] { + let (b, r) = hertz_cylinder(finite, 0.03); + let c = &r.contacts[0]; + let thickness = 0.2; + // The load per length: the reaction at the base over the slab's thickness, both halves. + let base: f64 = b.nodes(|x| near(x[1], -5.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()][1]).sum(); + let p = 2.0 * base / thickness; + let e_star = 1000.0 / (2.0 * (1.0 - 0.09)); + let half_width = (4.0 * p * 10.0 / (std::f64::consts::PI * e_star)).sqrt(); + let peak = 2.0 * p / (std::f64::consts::PI * half_width); + let found = c.pressure.iter().fold(0.0f64, |m, v| m.max(*v)); + let (before, next) = edge(&b, c, 2); + eprintln!("finite {finite}: P {p:.4} b {half_width:.4} ({before:.4}..{next:.4}) p0 {peak:.4} found {found:.4}"); + assert!((found - peak).abs() < 0.03 * peak, "p0 {found} vs {peak}"); + assert!(before <= half_width && half_width <= next, "b {half_width} not in ({before}, {next})"); + if finite { + quadratic(&r); + } + } +} + +/// The contact's edge along the x axis (the nodes on the plane where coordinate `across` is 0): +/// the node before the last closed one and the first open one after it, so that the edge lies +/// within an element of the last closed node (whose status the augmentation's tolerance decides). +fn edge(b: &Blocks, c: &fieldcraft_solver::contact::ContactResult, across: usize) -> (f64, f64) { + let mut along: Vec<(f64, Status)> = + c.nodes.iter().zip(&c.status).filter(|(l, _)| near(b.xyz[l][across], 0.0)).map(|(l, s)| (b.xyz[l][0], *s)).collect(); + along.sort_by(|x, y| x.0.total_cmp(&y.0)); + let last = along.iter().rposition(|(_, s)| *s == Status::Closed).unwrap_or(0); + let before = last.checked_sub(1).and_then(|i| along.get(i)).map_or(0.0, |p| p.0); + let next = along.get(last + 1).map_or(f64::INFINITY, |p| p.0); + (before, next) +} + +/// Newton's rate with contact: within each solution of an increment (its first, and each +/// augmentation's), once the residual is small every correction squares it, up to a constant +/// (the consistent tangent of the mortar integrals, segmentation and normal included). +fn quadratic(r: &StaticResult) { + let mut checked = 0; + for inc in &r.nonlinear.as_ref().unwrap().increments { + for w in inc.iterations.windows(2) { + let (a, b) = (w[0], w[1]); + // Below 1e-10 of the force scale the sums' round-off takes over. + if a.augmentation != b.augmentation || a.force > 1e-3 || b.force < 1e-10 { + continue; + } + assert!(b.force <= 100.0 * a.force * a.force, "increment {}: {:e} after {:e}", inc.number, b.force, a.force); + checked += 1; + } + } + assert!(checked >= 3, "too few iterations near the solution to judge ({checked})"); +} + +/// Two unit blocks, the upper one starting `overlap` into the lower one, pressed together. +fn stacked(overlap: f64, penetration: InitialPenetration) -> Blocks { + let mut b = Blocks::default(); + let low = b.block(ElementKind::Hex8, [2, 2, 1], |p| p); + let up = b.block(ElementKind::Hex8, [3, 3, 1], |p| [p[0], p[1], 1.0 - overlap + p[2]]); + b.define(1000.0, 0.3, NonlinearStatic { nonlinear_geometry: true, initial_increment: 1.0, ..NonlinearStatic::default() }); + b.face_set("lower_top", &low, |x| near(x[2], 1.0)); + b.face_set("upper_bottom", &up, |x| near(x[2], 1.0 - overlap)); + let mut c = Contact::new(ContactSurfaces::Pair { main: "lower_top".into(), secondary: "upper_bottom".into() }); + c.initial_penetration = penetration; + b.model.add_contact("interface", c).unwrap(); + b.bc("x0", [Some(0.0), None, None], |x| near(x[0], 0.0)); + b.bc("y0", [None, Some(0.0), None], |x| near(x[1], 0.0)); + b.bc("z0", [None, None, Some(0.0)], |x| near(x[2], 0.0)); + b.bc("push", [None, None, Some(-0.01)], move |x| near(x[2], 2.0 - overlap)); + b +} + +#[test] +fn initial_overclosure_is_adjusted_or_refused() { + let b = stacked(1e-3, InitialPenetration::Adjust); + let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap(); + assert!(r.warnings.iter().any(|w| w.contains("start inside the main surface")), "{:?}", r.warnings); + // Adjusted, the blocks behave as if they started touching: the same stress as without the + // overlap. + let touching = solve_nonlinear(&stacked(0.0, InitialPenetration::Adjust).model, Some("Step"), &Limits::default()).unwrap(); + let (a, t) = (r.elements[0].nodal[0][2], touching.elements[0].nodal[0][2]); + assert!((a - t).abs() < 1e-4 * t.abs(), "{a} vs {t}"); + let e = solve_nonlinear(&stacked(1e-3, InitialPenetration::Refuse).model, Some("Step"), &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("starts") && e.contains("inside"), "{e}"); +} + +#[test] +fn contact_results_do_not_depend_on_the_thread_count() { + let b = stacked(0.0, InitialPenetration::Adjust); + let run = |threads| fieldcraft_parallel::with_threads(threads, || solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap()); + let (one, many) = (run(1), run(4)); + assert_eq!(one.displacement, many.displacement); + assert_eq!(one.contacts, many.contacts); +} + +#[test] +fn what_implicit_contact_refuses() { + let solve = |m: &Model| solve_nonlinear(m, Some("Step"), &Limits::default()).unwrap_err().to_string(); + // Linear analyses refuse every contact by name. + let mut b = stacked(0.0, InitialPenetration::Adjust); + b.model.add_step(Step::new("Linear", Procedure::Static {}), None).unwrap(); + let e = fieldcraft_solver::solve_static(&b.model, Some("Linear"), &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("interface") && e.contains("linear static"), "{e}"); + // General contact. + let mut b = stacked(0.0, InitialPenetration::Adjust); + b.model.add_contact("all", Contact::new(ContactSurfaces::General { surface: None })).unwrap(); + assert!(solve(&b.model).contains("general")); + // Arc-length control. + let mut b = stacked(0.0, InitialPenetration::Adjust); + let arc = NonlinearStatic { arc_length: Some(ArcLength::default()), ..NonlinearStatic::default() }; + b.model.replace_step(Step { procedure: Procedure::NonlinearStatic(arc), ..b.model.step("Step").unwrap().clone() }).unwrap(); + assert!(solve(&b.model).contains("arc-length")); +} + +/// Hertz's point contact (Johnson 1985, §4.2): a sphere of radius R on a flat of the same +/// material under a load P has the contact radius a = (3 P R / (4 E*))^(1/3) and the peak +/// pressure p₀ = 3 P / (2 π a²), E* = E / (2 (1 − ν²)). A quarter model, displacement controlled; +/// the load is the reaction at the base. +#[test] +#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"] +fn hertz_sphere_on_a_flat() { + let (r, e, nu) = (10.0, 1000.0, 0.3); + let mut b = Blocks::default(); + let grade = |s: f64| s.powf(1.5); + let foundation = b.block(ElementKind::Hex8, [16, 16, 10], |p| [5.0 * grade(p[0]), 5.0 * grade(p[1]), -5.0 * grade(1.0 - p[2])]); + let cap = move |x: f64, y: f64| r - (r * r - x * x - y * y).sqrt(); + let indenter = b.block(ElementKind::Hex8, [14, 14, 8], |p| { + let (x, y) = (2.5 * grade(p[0]), 2.5 * grade(p[1])); + let bottom = cap(x, y); + [x, y, bottom + (3.0 - bottom) * grade(p[2])] + }); + b.define(e, nu, NonlinearStatic { nonlinear_geometry: false, initial_increment: 0.25, ..NonlinearStatic::default() }); + b.face_set("flat", &foundation, |x| near(x[2], 0.0)); + b.face_set("round", &indenter, move |x| (x[2] - cap(x[0], x[1])).abs() < 1e-9 && x[2] < 1.0); + b.pair("hertz", "flat", "round"); + b.bc("x0", [Some(0.0), None, None], |x| near(x[0], 0.0)); + b.bc("y0", [None, Some(0.0), None], |x| near(x[1], 0.0)); + b.bc("base", [None, None, Some(0.0)], |x| near(x[2], -5.0)); + b.bc("push", [None, None, Some(-0.04)], |x| near(x[2], 3.0)); + let res = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap(); + report(&res); + let c = &res.contacts[0]; + let base: f64 = b.nodes(|x| near(x[2], -5.0)).iter().map(|l| res.reaction[res.nodes.binary_search(l).unwrap()][2]).sum(); + let p = 4.0 * base; + let e_star = e / (2.0 * (1.0 - nu * nu)); + let a = (3.0 * p * r / (4.0 * e_star)).cbrt(); + let peak = 3.0 * p / (2.0 * std::f64::consts::PI * a * a); + let found = c.pressure.iter().fold(0.0f64, |m, v| m.max(*v)); + let (before, next) = edge(&b, c, 1); + eprintln!("P {p:.4} a {a:.4} ({before:.4}..{next:.4}) p0 {peak:.4} found {found:.4}"); + assert!((found - peak).abs() < 0.05 * peak, "p0 {found} vs {peak}"); + assert!(before <= a && a <= next, "a {a} not in ({before}, {next})"); +}