diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index e37cc8a..a2ce808 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -428,7 +428,7 @@ pub static CATALOGUE: &[Capability] = &[ "Contact Pair", Interactions, "Contact", - "Contact between two surfaces, with friction (explicit dynamics)", + "Contact between two surfaces: with friction in explicit dynamics, frictionless mortar contact in nonlinear static steps", &["interaction.contact.edit", "interaction.contact.delete"], ), ready( diff --git a/crates/engine/src/cmd/contact.rs b/crates/engine/src/cmd/contact.rs index 911c777..ad731ec 100644 --- a/crates/engine/src/cmd/contact.rs +++ b/crates/engine/src/cmd/contact.rs @@ -40,8 +40,8 @@ pub const PAIR: CommandSpec = CommandSpec { label: "Contact Pair", params: concat!( r#"{"name": "impact", "main": "plate" (a set of element faces), "secondary": "ball" (another), "#, - r#""friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb, μ = μk + (μs − μk) e^(−decay × slip rate); default frictionless), "penalty_scale": 1 (0 < scale ≤ 1: a factor on the stiffest penalty that leaves the time increment alone), "initial_penetration": "adjust" (move nodes that start inside onto the surface, strain-free, and report them) | "refuse" (refuse the analysis, naming the first)"#, - r#"} — contact between two surfaces, each one's nodes kept out of the other's faces (explicit dynamics; fails if the name is taken)"# + r#""friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb, μ = μk + (μs − μk) e^(−decay × slip rate); default frictionless), "penalty_scale": 1 (0 < scale ≤ 1: a factor on the penalty: explicit, the stiffest that leaves the time increment alone; nonlinear static, the augmented Lagrangian's), "initial_penetration": "adjust" (nodes that start inside: explicit moves them onto the surface, strain-free; nonlinear static takes the overlap out of their gap; both report them) | "refuse" (refuse the analysis, naming the first)"#, + r#"} — contact between two surfaces: in explicit dynamics each one's nodes are kept out of the other's faces; in nonlinear static steps `secondary` is the mortar (non-mortar) side whose pressure field is enforced against `main` by an augmented Lagrangian, frictionless so far (fails if the name is taken)"# ), form: Some(super::forms::CONTACT_PAIR), kind: CommandKind::Model, @@ -58,8 +58,8 @@ pub const GENERAL: CommandSpec = CommandSpec { label: "General Contact", params: concat!( r#"{"name": "all", "surface": "skin" (optional: a set of element faces, each in contact with every other: self contact; default every exterior face of the solids), "#, - r#""friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb, μ = μk + (μs − μk) e^(−decay × slip rate); default frictionless), "penalty_scale": 1 (0 < scale ≤ 1: a factor on the stiffest penalty that leaves the time increment alone), "initial_penetration": "adjust" (move nodes that start inside onto the surface, strain-free, and report them) | "refuse" (refuse the analysis, naming the first)"#, - r#"} — single-surface contact of everything that may touch (explicit dynamics; fails if the name is taken)"# + r#""friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb, μ = μk + (μs − μk) e^(−decay × slip rate); default frictionless), "penalty_scale": 1 (0 < scale ≤ 1: a factor on the penalty: explicit, the stiffest that leaves the time increment alone; nonlinear static, the augmented Lagrangian's), "initial_penetration": "adjust" (nodes that start inside: explicit moves them onto the surface, strain-free; nonlinear static takes the overlap out of their gap; both report them) | "refuse" (refuse the analysis, naming the first)"#, + r#"} — single-surface contact of everything that may touch (explicit dynamics; nonlinear static steps refuse it: use pairs; fails if the name is taken)"# ), form: Some(super::forms::CONTACT_GENERAL), kind: CommandKind::Model, diff --git a/crates/engine/src/cmd/forms.rs b/crates/engine/src/cmd/forms.rs index dfeee85..1bee87b 100644 --- a/crates/engine/src/cmd/forms.rs +++ b/crates/engine/src/cmd/forms.rs @@ -880,6 +880,9 @@ pub const CONTOUR: &[Field] = &[ ("stress_xy", "Stress XY"), ("stress_xz", "Stress XZ"), ("stress_yz", "Stress YZ"), + ("contact_pressure", "Contact pressure"), + ("contact_gap", "Contact gap"), + ("contact_status", "Contact status"), ("none", "None"), ], ) diff --git a/crates/engine/src/cmd/results.rs b/crates/engine/src/cmd/results.rs index cac209f..e9db9ec 100644 --- a/crates/engine/src/cmd/results.rs +++ b/crates/engine/src/cmd/results.rs @@ -20,7 +20,7 @@ pub const DEFORMED: CommandSpec = CommandSpec { pub const CONTOUR: CommandSpec = CommandSpec { id: "results.contour", label: "Contours", - params: r#"{"step": "Load", "field": "von_mises", "bands": 12} — contours of a nodal field: displacement, displacement_x/y/z, von_mises, max/mid/min_principal, stress_xx/yy/zz/xy/xz/yz (stresses for static steps); "none" turns them off; 2 to 24 bands"#, + params: r#"{"step": "Load", "field": "von_mises", "bands": 12} — contours of a nodal field: displacement, displacement_x/y/z, von_mises, max/mid/min_principal, stress_xx/yy/zz/xy/xz/yz (stresses for static steps), contact_pressure, contact_gap, contact_status (0 open, 1 closed; nonlinear static steps with contact, at the contacts' secondary surfaces); "none" turns them off; 2 to 24 bands"#, form: Some(super::forms::CONTOUR), kind: CommandKind::Session, run: contour, diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 737f534..a367510 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -36,7 +36,7 @@ pub const STATIC: CommandSpec = CommandSpec { pub const NONLINEAR: CommandSpec = CommandSpec { id: "solve.nonlinear", label: "Solve Nonlinear Static Step", - params: r#"{"step": "Inflate", "max_dofs": 6000000, "max_factor_entries": 400000000, "max_newton_iterations": 100000, "threads": 8} — solves a nonlinear static step (the only one when `step` is left out) in load increments with Newton–Raphson iterations and keeps the result in the session: the final displacements, reactions, Cauchy stresses and logarithmic strains (small strains without nonlinear geometry), the equivalent plastic strain of plastic materials (the largest in the summary), every increment's load factor, iterations and residuals (force, displacement and energy, relative to their tolerances' scales), and a displacement frame per converged increment; the caps may only be lowered (`max_newton_iterations` over the whole step); `threads` as for solve.static"#, + params: r#"{"step": "Inflate", "max_dofs": 6000000, "max_factor_entries": 400000000, "max_newton_iterations": 100000, "threads": 8} — solves a nonlinear static step (the only one when `step` is left out) in load increments with Newton–Raphson iterations and keeps the result in the session: the final displacements, reactions, Cauchy stresses and logarithmic strains (small strains without nonlinear geometry), the equivalent plastic strain of plastic materials (the largest in the summary), every increment's load factor, iterations and residuals (force, displacement and energy, relative to their tolerances' scales), and a displacement frame per converged increment, and with contact pairs (interaction.contact.pair: mortar contact by augmented Lagrangian, frictionless) each contact's normal force, largest pressure, closed nodes, deepest penetration and augmentations (solve.probe gives each secondary node's pressure, gap and status); the caps may only be lowered (`max_newton_iterations` over the whole step); `threads` as for solve.static"#, form: Some(super::forms::SOLVE_NONLINEAR), kind: CommandKind::Model, run: run_nonlinear, @@ -239,6 +239,15 @@ fn run_nonlinear(s: &mut Session, p: &Value) -> Result { "frames": history.map(|h| h.load_factors()), "load_factor": history.and_then(|h| h.frames.last()).map(|f| f.load_factor), "critical_points": history.map(|h| &h.critical_points), + "contacts": r.contacts.iter().map(|c| json!({ + "name": c.name, + "normal_force": c.normal_force, + "force": c.force, + "max_pressure": largest(c.nodes.iter().copied().zip(c.pressure.iter().copied())), + "closed_nodes": c.closed(), + "penetration": c.penetration, + "augmentations": c.augmentations, + })).collect::>(), "max_displacement": largest(r.nodes.iter().zip(&r.displacement).map(|(l, u)| (*l, crate::results::magnitude(*u)))), "max_von_mises": largest(r.nodes.iter().zip(&r.von_mises).map(|(l, v)| (*l, *v))), "max_plastic_strain": largest(r.elements.iter().flat_map(|e| e.plastic_strain.iter().map(|v| (e.label, *v)))), @@ -609,16 +618,31 @@ fn probe_implicit(result: &ImplicitResult, step: &str, nodes: &[u32], p: &Value) let p = r.reference_points.get(k); json!({"node": l, "reference_point": p.map(|p| &p.name), "displacement": p.map(|p| p.displacement), "reaction": p.map(|p| p.reaction)}) } - Location::Node(i) => json!({ - "node": l, - "displacement": r.displacement.get(i), - "rotation": r.rotation.get(i), - "reaction": r.reaction.get(i), - "reaction_moment": r.reaction_moment.get(i), - "stress": r.stress.get(i), - "von_mises": r.von_mises.get(i), - "temperature": r.temperature.as_ref().and_then(|t| t.get(i)), - }), + Location::Node(i) => { + let mut v = json!({ + "node": l, + "displacement": r.displacement.get(i), + "rotation": r.rotation.get(i), + "reaction": r.reaction.get(i), + "reaction_moment": r.reaction_moment.get(i), + "stress": r.stress.get(i), + "von_mises": r.von_mises.get(i), + "temperature": r.temperature.as_ref().and_then(|t| t.get(i)), + }); + // The contacts the node is a secondary node of: pressure, gap, status. + let contact: Vec = r + .contacts + .iter() + .filter_map(|c| { + let k = c.nodes.binary_search(&l).ok()?; + Some(json!({"contact": c.name, "pressure": c.pressure.get(k), "gap": c.gap.get(k), "status": c.status.get(k)})) + }) + .collect(); + if !contact.is_empty() { + v["contact"] = Value::Array(contact); + } + v + } }) }) .collect::>()?, @@ -1509,4 +1533,59 @@ pub(crate) mod tests { let all = s.execute("solve.run", &json!({})).unwrap(); assert_eq!(all["steps"].as_array().map(Vec::len), Some(3), "{all}"); } + + /// Two blocks pressed together across a contact pair, as commands. + fn pressed_blocks(s: &mut Session) { + let script = json!([ + {"command": "mesh.block", "params": {"body": "Lower", "kind": "hex8", "corners": [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]], "divisions": [2, 2, 1], "faces": {"1": "bottom", "2": "top"}}}, + {"command": "mesh.block", "params": {"body": "Upper", "kind": "hex8", "corners": [[0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1], [0, 0, 2], [1, 0, 2], [1, 1, 2], [0, 1, 2]], "divisions": [3, 3, 1], "faces": {"1": "bottom", "2": "top"}, "merge": 0}}, + {"command": "model.set.create", "params": {"name": "all", "geometry": [{"body": "Lower"}, {"body": "Upper"}]}}, + {"command": "model.set.create", "params": {"name": "lower_top", "geometry": [{"body": "Lower", "face": "top"}]}}, + {"command": "model.set.create", "params": {"name": "upper_bottom", "geometry": [{"body": "Upper", "face": "bottom"}]}}, + {"command": "model.set.create", "params": {"name": "base", "geometry": [{"body": "Lower", "face": "bottom"}]}}, + {"command": "model.set.create", "params": {"name": "top", "geometry": [{"body": "Upper", "face": "top"}]}}, + {"command": "model.material.create", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 1000.0, "poisson": 0.0}}}, + {"command": "model.section.create", "params": {"name": "all", "set": "all", "material": "Mat", "kind": {"type": "solid"}}}, + {"command": "interaction.contact.pair", "params": {"name": "interface", "main": "lower_top", "secondary": "upper_bottom"}}, + {"command": "model.boundary_condition.create", "params": {"name": "base", "type": "fixed", "target": {"set": "base"}}}, + {"command": "model.boundary_condition.create", "params": {"name": "push", "type": "displacement", "target": {"set": "top"}, "values": [0, 0, -0.02, null, null, null]}}, + {"command": "model.step.create", "params": {"name": "Squeeze", "procedure": {"type": "nonlinear_static", "nonlinear_geometry": false}, "boundary_conditions": ["base", "push"], "loads": []}} + ]); + s.run_script(&script).unwrap(); + } + + #[test] + fn solve_and_probe_a_contact() { + let mut s = Session::default(); + pressed_blocks(&mut s); + let r = s.execute("solve.nonlinear", &json!({"step": "Squeeze"})).unwrap(); + let c = &r["contacts"][0]; + assert_eq!(c["name"], "interface"); + // ν = 0 and held ends: uniaxial strain 1% of E = 1000, over the unit area. + let force = c["normal_force"].as_f64().unwrap(); + assert!((force - 10.0).abs() < 2e-3, "{r}"); + assert_eq!(c["closed_nodes"], 16, "{r}"); + // A secondary node gives its contact pressure, gap and status; a main one gives none. + let third = 1.0 / 3.0; + let probe = s.execute("solve.probe", &json!({"step": "Squeeze", "at": [[third, third, 1.0], [0.5, 0.5, 1.0]]})).unwrap(); + let contact = &probe["nodes"][0]["contact"][0]; + assert_eq!(contact["status"], "closed", "{probe}"); + assert!((contact["pressure"].as_f64().unwrap() - 10.0).abs() < 2e-3, "{probe}"); + assert!(contact["gap"].as_f64().unwrap() < 0.0, "{probe}"); + assert!(probe["nodes"][1].get("contact").is_none(), "{probe}"); + // The contact fields on the results display. + let v = s.execute("results.contour", &json!({"step": "Squeeze", "field": "contact_pressure"})).unwrap(); + let range = &v["contour"]["range"]; + assert!((range["max"]["value"].as_f64().unwrap() - 10.0).abs() < 2e-3 && range["min"]["value"] == 0.0, "{v}"); + let status = s.execute("results.contour", &json!({"step": "Squeeze", "field": "contact_status"})).unwrap(); + assert_eq!(status["contour"]["range"]["max"]["value"], 1.0, "{status}"); + // The linear analyses refuse the contact by name. + s.execute( + "model.step.create", + &json!({"name": "Linear", "procedure": {"type": "static"}, "boundary_conditions": ["base", "push"], "loads": []}), + ) + .unwrap(); + let e = s.execute("solve.static", &json!({"step": "Linear"})).unwrap_err().to_string(); + assert!(e.contains("interface"), "{e}"); + } } diff --git a/crates/engine/src/results.rs b/crates/engine/src/results.rs index 22baa0b..7cd9866 100644 --- a/crates/engine/src/results.rs +++ b/crates/engine/src/results.rs @@ -87,10 +87,17 @@ pub enum Field { Temperature, /// Heat flux magnitude (heat transfer steps, at the end of the step). HeatFlux, + /// Contact pressure (nonlinear static steps with contact: 0 off the contacts' secondary + /// surfaces). + ContactPressure, + /// The normalised contact gap of the secondary surfaces' nodes (negative: penetration). + ContactGap, + /// Contact status of the secondary surfaces' nodes: 0 open, 1 closed. + ContactStatus, } impl Field { - pub const ALL: [Field; 17] = [ + pub const ALL: [Field; 20] = [ Field::Displacement, Field::DisplacementX, Field::DisplacementY, @@ -108,6 +115,9 @@ impl Field { Field::StressYz, Field::Temperature, Field::HeatFlux, + Field::ContactPressure, + Field::ContactGap, + Field::ContactStatus, ]; /// The id commands take (`von_mises`…). @@ -130,6 +140,9 @@ impl Field { Field::StressYz => "stress_yz", Field::Temperature => "temperature", Field::HeatFlux => "heat_flux", + Field::ContactPressure => "contact_pressure", + Field::ContactGap => "contact_gap", + Field::ContactStatus => "contact_status", } } @@ -152,6 +165,9 @@ impl Field { Field::StressYz => "Stress YZ", Field::Temperature => "Temperature", Field::HeatFlux => "Heat flux magnitude", + Field::ContactPressure => "Contact pressure", + Field::ContactGap => "Contact gap", + Field::ContactStatus => "Contact status (0 open, 1 closed)", } } @@ -170,7 +186,12 @@ impl Field { | Field::Velocity | Field::Temperature | Field::HeatFlux - ) + ) && !self.is_contact() + } + + /// Whether it is a contact quantity. + pub fn is_contact(self) -> bool { + matches!(self, Field::ContactPressure | Field::ContactGap | Field::ContactStatus) } /// Whether it is a heat transfer quantity. @@ -181,6 +202,7 @@ impl Field { /// The fields a result has: displacements everywhere, stresses for static steps and explicit /// frames that kept them (mode shapes have no stress scale), velocity for explicit frames. pub fn available(r: &StepResult) -> impl Iterator + '_ { + let contact = matches!(r.kind(), Kind::Static(s) if !s.contacts.is_empty()); let (stress, velocity, temperature, flux) = match r.kind() { Kind::Static(s) => (true, false, s.temperature.is_some(), false), Kind::Frequency(_) => (false, false, false, false), @@ -197,6 +219,7 @@ impl Field { Field::Velocity => velocity, Field::Temperature => temperature, Field::HeatFlux => flux, + f if f.is_contact() => contact, f => stress || !f.is_stress(), }) .collect::>() @@ -333,6 +356,10 @@ pub fn field_values(model: &Model, r: &StepResult, field: Field, pick: Pick) -> Kind::Heat(h) => Ok(h.flux.iter().map(|v| magnitude(*v)).collect()), _ => Err(unavailable()), }, + Field::ContactPressure | Field::ContactGap | Field::ContactStatus => match r.kind() { + Kind::Static(s) if !s.contacts.is_empty() => Ok(contact_values(s, field)), + _ => Err(unavailable()), + }, Field::VonMises => each(&fieldcraft_fem::von_mises), Field::MaxPrincipal => principal(0), Field::MidPrincipal => principal(1), @@ -346,6 +373,32 @@ pub fn field_values(model: &Model, r: &StepResult, field: Field, pick: Pick) -> } } +/// A contact field at every node of a static result: the largest pressure, the smallest gap or +/// the closed status over the contacts a node is a secondary node of; off them 0 pressure and no +/// gap or status (NaN). +fn contact_values(s: &StaticResult, field: Field) -> Vec { + let mut out = vec![if field == Field::ContactPressure { 0.0 } else { f64::NAN }; s.nodes.len()]; + for c in &s.contacts { + for (k, label) in c.nodes.iter().enumerate() { + let Some(v) = s.nodes.binary_search(label).ok().and_then(|i| out.get_mut(i)) else { continue }; + match field { + Field::ContactPressure => *v = v.max(c.pressure.get(k).copied().unwrap_or(0.0)), + Field::ContactGap => { + if let Some(Some(g)) = c.gap.get(k) { + *v = if v.is_nan() { *g } else { v.min(*g) }; + } + } + Field::ContactStatus => { + let closed = c.status.get(k).is_some_and(|s| *s == fieldcraft_solver::contact::Status::Closed); + *v = if closed || (!v.is_nan() && *v > 0.0) { 1.0 } else { 0.0 }; + } + _ => {} + } + } + } + out +} + /// Element values (in `elements` order) averaged at the nodes (in `nodes` order) of the /// elements around them; a node no element reaches gets NaN. fn nodal_average(model: &Model, nodes: &[u32], elements: &[u32], values: &[[f64; 6]]) -> Vec<[f64; 6]> { @@ -472,7 +525,7 @@ mod tests { } let (modes, _) = s.latest_result("Modes").unwrap(); assert_eq!(Field::available(modes).count(), 4); - assert_eq!(Field::available(bend).count(), Field::ALL.len() - 3, "all but velocity and the thermal fields"); + assert_eq!(Field::available(bend).count(), Field::ALL.len() - 6, "all but velocity, the thermal and the contact fields"); assert!(field_values(&s.model, modes, Field::VonMises, Pick { mode: 1, frame: None }).is_err()); assert!(field_values(&s.model, modes, Field::Displacement, Pick { mode: 9, frame: None }).is_err()); assert_eq!(field_values(&s.model, modes, Field::Displacement, Pick { mode: 2, frame: None }).unwrap().len(), modes.nodes().len()); diff --git a/crates/model/src/contact.rs b/crates/model/src/contact.rs index 1cf8e09..0f992c1 100644 --- a/crates/model/src/contact.rs +++ b/crates/model/src/contact.rs @@ -1,7 +1,8 @@ //! Contact: which surfaces may touch, and the friction between them. //! -//! A contact is a model definition, active in every step that can take it (explicit dynamics; -//! the implicit solver refuses contact by name for now). Its surfaces are sets with element +//! A contact is a model definition, active in every step that can take it (explicit dynamics, +//! and nonlinear static steps for pairs; the linear analyses refuse contact by name). Its +//! surfaces are sets with element //! faces, so they survive remeshing like a pressure load does: //! - a **pair**: the faces of `main` against those of `secondary`, both ways (each surface's //! nodes are kept out of the other's faces, so neither surface is privileged); @@ -17,10 +18,13 @@ //! //! **Penalty.** Explicit contact pushes penetrating nodes back with a penalty spring whose //! stiffness the solver chooses from the masses and the stable time increment, the stiffest -//! that leaves the increment alone; `penalty_scale` (at most 1) softens it. +//! that leaves the increment alone; `penalty_scale` (at most 1) softens it. Nonlinear static +//! steps enforce a pair's mortar pressure by an augmented Lagrangian whose penalty the solver +//! chooses from the elements under the secondary surface; `penalty_scale` scales it too. //! //! **Initial penetration.** Nodes of a surface that start inside the other surface are either -//! moved onto it before the analysis, strain-free (`adjust`, the default, each one reported), or +//! moved onto it before the analysis, strain-free (`adjust`, the default, each one reported; +//! nonlinear static steps take the overlap out of their gap instead), or //! the analysis is refused naming the first node, its depth and the face it is behind //! (`refuse`). diff --git a/docs/features.md b/docs/features.md index 8d7bb7f..de29df0 100644 --- a/docs/features.md +++ b/docs/features.md @@ -104,7 +104,7 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands. | Constraints | Coupling — Kinematic or distributing coupling of a set to a reference point | ready | `model.coupling.create`, `model.coupling.edit`, `model.coupling.delete` | | Constraints | Tie — Tie non-matching meshes together | coming soon (M7) | | | Contact | General Contact — Contact between everything that may touch, or self contact of a surface, with friction (explicit dynamics) | ready | `interaction.contact.general`, `interaction.contact.edit`, `interaction.contact.delete` | -| Contact | Contact Pair — Contact between two surfaces, with friction (explicit dynamics) | ready | `interaction.contact.pair`, `interaction.contact.edit`, `interaction.contact.delete` | +| Contact | Contact Pair — Contact between two surfaces: with friction in explicit dynamics, frictionless mortar contact in nonlinear static steps | ready | `interaction.contact.pair`, `interaction.contact.edit`, `interaction.contact.delete` | | Contact | Find Contacts — Find the parts of the mesh that touch, and make contact pairs of them | ready | `interaction.contact.find` | | Contact | Rigid Wall — Rigid walls for explicit dynamics | coming soon (M10) | | | Connectors | Bolt Load — Bolt pretension | coming soon (M7) | |