diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 99e3051..e049cee 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -13,8 +13,20 @@ //! //! In order: heading, nodes (mesh nodes, then reference points), elements grouped by element type, //! sets (a node set, an element set and an element-face surface per model set, each written when -//! the set resolves to some), materials, sections, couplings, the initial state's boundary -//! conditions, then the steps. +//! the set resolves to some), materials, sections, couplings, contact, the initial state's +//! boundary conditions, then the steps. +//! +//! # Contact +//! +//! Each contact is a `*SURFACE INTERACTION` with its friction (`*FRICTION`, or +//! `*FRICTION, EXPONENTIAL DECAY` with the static and kinetic coefficients and the decay +//! coefficient, the same law). A pair is `*CONTACT PAIR, TYPE=SURFACE TO SURFACE` with the +//! secondary surface first, and `ADJUST=0.0` when initially overclosed nodes are adjusted. General +//! contact is the model's one `*CONTACT` definition: `*CONTACT INCLUSIONS, ALL EXTERIOR` for every +//! exterior face, a surface with a blank partner for self contact, and a +//! `*CONTACT PROPERTY ASSIGNMENT` line per contact. FieldCraft's penalty scale and its refusal of +//! initial overclosure have no Abaqus counterpart and are noted in comments (Abaqus chooses its +//! own penalty and resolves overclosure in the first increment). //! //! # Steps and the activity table //! @@ -130,8 +142,8 @@ use std::collections::{BTreeMap, BTreeSet}; use std::fmt::Write as _; use fieldcraft_model::{ - BeamProfile, BoundaryCondition, CouplingKind, Elastic, Element, ElementKind, Hardening, Hyperelastic, Integration, Load, Model, OutputVariable, - Plastic, Procedure, Resolved, Section, SectionKind, Step, Symmetry, Target, + BeamProfile, BoundaryCondition, ContactSurfaces, CouplingKind, Elastic, Element, ElementKind, Hardening, Hyperelastic, InitialPenetration, + Integration, Load, Model, OutputVariable, Plastic, Procedure, Resolved, Section, SectionKind, Step, Symmetry, Target, }; use crate::{Error, Result, unsupported}; @@ -174,6 +186,7 @@ pub fn write_resolved(model: &Model, r: &Resolved) -> Result { d.sections(model)?; d.discrete_properties(&discrete)?; d.couplings(model, r); + d.contacts(model, r)?; d.steps(model, r)?; Ok(d.out) } @@ -957,6 +970,78 @@ impl Deck { } } + /// One surface interaction (the friction) per contact, then its pair, and every general + /// contact in the model's one `*CONTACT` definition. + fn contacts(&mut self, model: &Model, r: &Resolved) -> Result<()> { + if model.contacts().is_empty() { + return Ok(()); + } + let mut taken: BTreeSet = model.sets().keys().map(|n| n.to_ascii_uppercase()).collect(); + taken.extend(coupling_surfaces(model, r).into_values().map(|n| n.to_ascii_uppercase())); + self.comment("Contact: a surface interaction (Coulomb friction) per contact, contact pairs, general contact"); + let mut general = Vec::new(); + for (name, c) in model.contacts() { + let interaction = crate::unique_name(&format!("{name}-interaction"), &mut taken); + self.line(&format!("*SURFACE INTERACTION, NAME={interaction}")); + let f = c.friction; + if f.decay > 0.0 { + self.line("*FRICTION, EXPONENTIAL DECAY"); + self.numbers_line("", &[f.static_coefficient, f.kinetic_coefficient(), f.decay])?; + } else { + self.line("*FRICTION"); + self.numbers_line("", &[f.static_coefficient])?; + } + if c.penalty_scale != 1.0 { + self.comment(&format!( + "contact `{name}`: penalty scale {} of FieldCraft's explicit penalty stiffness; Abaqus chooses its own", + c.penalty_scale + )); + } + match &c.surfaces { + ContactSurfaces::Pair { main, secondary } => { + // ADJUST=0.0 moves the secondary nodes that start inside the main surface + // onto it, strain-free; Abaqus has no refusal of initial overclosure. + let adjust = match c.initial_penetration { + InitialPenetration::Adjust => ", ADJUST=0.0", + InitialPenetration::Refuse => { + self.comment(&format!( + "contact `{name}`: FieldCraft refuses initially overclosed nodes; Abaqus resolves them in the first increment" + )); + "" + } + }; + self.line(&format!("*CONTACT PAIR, INTERACTION={interaction}, TYPE=SURFACE TO SURFACE{adjust}")); + self.line(&format!("{secondary}, {main}")); + } + ContactSurfaces::General { surface } => general.push((surface.clone(), interaction)), + } + } + if general.is_empty() { + return Ok(()); + } + // General contact: every surface against itself (a blank second surface), or every + // exterior face; each with its interaction. + self.line("*CONTACT"); + if general.iter().any(|(s, _)| s.is_none()) { + self.line("*CONTACT INCLUSIONS, ALL EXTERIOR"); + } + let named: Vec<&str> = general.iter().filter_map(|(s, _)| s.as_deref()).collect(); + if !named.is_empty() { + self.line("*CONTACT INCLUSIONS"); + for s in &named { + self.line(&format!("{s},")); + } + } + self.line("*CONTACT PROPERTY ASSIGNMENT"); + for (surface, interaction) in &general { + match surface { + Some(s) => self.line(&format!("{s}, {s}, {interaction}")), + None => self.line(&format!(", , {interaction}")), + } + } + Ok(()) + } + fn boundary_lines<'n>(&mut self, model: &Model, r: &Resolved, names: impl IntoIterator) -> Result<()> { for name in names { let Some(bc) = model.boundary_condition(name) else { continue }; diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 40d26c4..6da44cf 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -22,6 +22,12 @@ //! - couplings: kinematic as `*CONSTRAINED_NODAL_RIGID_BODY` (the reference point first in its //! node set), distributing as `*CONSTRAINED_INTERPOLATION` (the reference point dependent on //! the set's translations, equal weights); +//! - contact: a pair as `*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE` between segment sets (the +//! secondary surface first), general contact as `*CONTACT_AUTOMATIC_SINGLE_SURFACE` on a segment +//! set or on every part (`surfatyp` 5); the static and kinetic coefficients and the decay +//! coefficient (`fs`, `fd`, `dc`: the same exponential decay law) and the penalty scale on both +//! sides (`sfsa`, `sfsb`, a factor on LS-DYNA's own default stiffness). LS-DYNA moves nodes +//! that start inside the other surface out, as `adjust` does; a refusal is noted in a comment; //! - per step: `*CONTROL_IMPLICIT_GENERAL`, `*CONTROL_IMPLICIT_SOLUTION` (linear) or //! `*CONTROL_IMPLICIT_EIGENVALUE`, `*CONTROL_TERMINATION` (a static step is one implicit step //! from time 0 to 1 with every load ramped by one curve, `*DEFINE_CURVE`); a nonlinear static @@ -69,8 +75,8 @@ mod elements; use std::collections::{BTreeMap, BTreeSet}; use fieldcraft_model::{ - BoundaryCondition, CouplingKind, Diagnostic, Dofs, Elastic, Element, ElementKind, Hardening, Load, Material, Model, OutputVariable, Plastic, - Procedure, Resolved, ResolvedDirection, Section, SectionKind, Step, Target, + BoundaryCondition, ContactSurfaces, CouplingKind, Diagnostic, Dofs, Elastic, Element, ElementKind, Hardening, InitialPenetration, Load, Material, + Model, OutputVariable, Plastic, Procedure, Resolved, ResolvedDirection, Section, SectionKind, Step, Target, }; use crate::{Error, Result}; @@ -177,6 +183,7 @@ fn write_with(model: &Model, r: &Resolved) -> Result { element_cards(&mut d, &plan)?; sets(&mut d, model, r, &plan)?; couplings(&mut d, model, r, &plan)?; + contacts(&mut d, model, &plan)?; discrete(&mut d, model, r, &plan)?; if let Some(parts) = &plan.gravity_parts { d.line("*SET_PART_LIST"); @@ -269,6 +276,10 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> { node_sets.insert(c.set.as_str()); } } + // Contact surfaces are segment sets. + for c in model.contacts().values() { + segment_sets.extend(c.surfaces.sets()); + } let mut gravity_sets: BTreeSet> = BTreeSet::new(); for step in model.steps() { if matches!(step.procedure, Procedure::ExplicitDynamic(_)) { @@ -689,6 +700,42 @@ fn couplings(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan) -> Result<( Ok(()) } +/// Contacts (see the module docs): a pair as an automatic surface-to-surface contact between the +/// segment sets, general contact as an automatic single-surface contact. +fn contacts(d: &mut Deck, model: &Model, plan: &Plan) -> Result<()> { + for (name, c) in model.contacts() { + let sid = |set: &str| plan.sets.get(set).copied().ok_or_else(|| broken(format!("contact `{name}`: set `{set}`"))); + d.comment(&format!("contact `{name}`")); + if c.initial_penetration == InitialPenetration::Refuse { + d.comment("FieldCraft refuses nodes that start inside the other surface; LS-DYNA moves them out (as `adjust` does)"); + } + match &c.surfaces { + ContactSurfaces::Pair { main, secondary } => { + d.line("*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE"); + d.header(&["surfa", "surfb", "surfatyp", "surfbtyp"]); + d.card([Field::from(sid(secondary)?), Field::from(sid(main)?), Field::from(0), Field::from(0)]); + } + ContactSurfaces::General { surface: Some(s) } => { + d.line("*CONTACT_AUTOMATIC_SINGLE_SURFACE"); + d.header(&["surfa", "surfb", "surfatyp"]); + d.card([Field::from(sid(s)?), Field::Blank, Field::from(0)]); + } + ContactSurfaces::General { surface: None } => { + // SURFATYP 5: every part. + d.line("*CONTACT_AUTOMATIC_SINGLE_SURFACE"); + d.header(&["surfa", "surfb", "surfatyp"]); + d.card([Field::from(0), Field::Blank, Field::from(5)]); + } + } + let f = c.friction; + d.header(&["fs", "fd", "dc"]); + d.card([Field::from(f.static_coefficient), Field::from(f.kinetic_coefficient()), Field::from(f.decay)]); + d.header(&["sfsa", "sfsb"]); + d.card([Field::from(c.penalty_scale), Field::from(c.penalty_scale)]); + } + Ok(()) +} + /// Point masses, springs and dashpots (see the module docs). fn discrete(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan) -> Result<()> { if r.masses.is_empty() && r.springs.is_empty() { diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index 28d171c..3e53e12 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -331,6 +331,9 @@ impl<'m> Deck<'m> { self.plan_constraints()?; self.plan_loads()?; self.problems.extend(crate::thermal_items(self.model)); + for name in self.model.contacts().keys() { + self.problems.insert(format!("contact `{name}` (Nastran's linear solutions SOL 101 and 103 have no contact)")); + } for sub in &mut self.subcases { if matches!(sub.step.procedure, Procedure::Frequency { .. }) { sub.method = Some(next_id(&mut self.next_sid)?); diff --git a/crates/decks/tests/contact.rs b/crates/decks/tests/contact.rs new file mode 100644 index 0000000..1021ff7 --- /dev/null +++ b/crates/decks/tests/contact.rs @@ -0,0 +1,98 @@ +//! Contact definitions in every deck format: Abaqus surface interactions, contact pairs and +//! general contact; LS-DYNA automatic contacts; Nastran refuses them by name. + +use fieldcraft_decks::{Error, abaqus, lsdyna, nastran}; +use fieldcraft_model::{ + BoundaryCondition, Contact, ContactSurfaces, Elastic, ElementFace, ElementKind, Friction, InitialPenetration, Material, Model, Procedure, + Section, SectionKind, Set, Step, Target, UnitSystem, +}; + +/// Two unit hex8 cubes, the second on top of the first (separate nodes), with face sets for the +/// touching faces and a static step. +fn model() -> Model { + let mut m = Model::default(); + m.set_units(UnitSystem::MmTonneS); + let 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.]]; + for cube in 0..2u32 { + for (i, c) in corners.iter().enumerate() { + m.set_node(8 * cube + i as u32 + 1, [c[0], c[1], c[2] + cube as f64]).unwrap(); + } + m.set_element(cube + 1, ElementKind::Hex8, (8 * cube + 1..=8 * cube + 8).collect()).unwrap(); + } + m.add_set("all", Set { elements: vec![1, 2], ..Set::default() }).unwrap(); + m.add_set("lower", Set { faces: vec![ElementFace { element: 1, face: 2 }], ..Set::default() }).unwrap(); + m.add_set("upper", Set { faces: vec![ElementFace { element: 2, face: 1 }], ..Set::default() }).unwrap(); + m.add_set("base", Set { nodes: vec![1, 2, 3, 4], ..Set::default() }).unwrap(); + m.add_material( + "Steel", + Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), ..Material::default() }, + ) + .unwrap(); + m.add_section("all", Section { set: "all".into(), material: "Steel".into(), kind: SectionKind::Solid {}, formulation: Default::default() }) + .unwrap(); + m.add_boundary_condition("fix", BoundaryCondition::Fixed { target: Target::Set("base".into()) }).unwrap(); + let mut step = Step::new("Load", Procedure::Static {}); + step.boundary_conditions.insert("fix".into()); + m.add_step(step, None).unwrap(); + let pair = Contact { + friction: Friction { static_coefficient: 0.3, kinetic: Some(0.2), decay: 5.0 }, + initial_penetration: InitialPenetration::Adjust, + ..Contact::new(ContactSurfaces::Pair { main: "lower".into(), secondary: "upper".into() }) + }; + m.add_contact("stack", pair).unwrap(); + m.add_contact("everything", Contact { friction: Friction::coulomb(0.1), ..Contact::new(ContactSurfaces::General { surface: None }) }).unwrap(); + m.add_contact("fold", Contact { penalty_scale: 0.5, ..Contact::new(ContactSurfaces::General { surface: Some("upper".into()) }) }).unwrap(); + m +} + +#[test] +fn abaqus_writes_interactions_pairs_and_general_contact() { + let deck = abaqus::write(&model()).unwrap(); + for expected in [ + "*SURFACE, TYPE=ELEMENT, NAME=lower\n1, S2\n", + "*SURFACE INTERACTION, NAME=stack-interaction\n*FRICTION, EXPONENTIAL DECAY\n0.3, 0.2, 5.0\n", + "*CONTACT PAIR, INTERACTION=stack-interaction, TYPE=SURFACE TO SURFACE, ADJUST=0.0\nupper, lower\n", + "*SURFACE INTERACTION, NAME=everything-interaction\n*FRICTION\n0.1\n", + "*SURFACE INTERACTION, NAME=fold-interaction\n*FRICTION\n0.0\n", + "*CONTACT\n*CONTACT INCLUSIONS, ALL EXTERIOR\n*CONTACT INCLUSIONS\nupper,\n*CONTACT PROPERTY ASSIGNMENT\n, , everything-interaction\nupper, upper, fold-interaction\n", + "** contact `fold`: penalty scale 0.5", + ] { + assert!(deck.contains(expected), "{expected}\n{deck}"); + } + // Model data: before the first step. + assert!(deck.find("*CONTACT PAIR") < deck.find("*STEP")); +} + +#[test] +fn lsdyna_writes_automatic_contacts_on_segment_sets() { + let deck = lsdyna::write(&model()).unwrap().text; + let after = |keyword: &str| -> Vec { + let at = deck.find(keyword).unwrap_or_else(|| panic!("{keyword}\n{deck}")); + deck[at..].lines().skip(1).filter(|l| !l.starts_with('$')).take(3).map(|l| l.split_whitespace().collect::>().join(" ")).collect() + }; + // Sets are numbered in name order: all 1, base 2, lower 3, upper 4. + assert_eq!(after("*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE"), ["4 3 0 0", "0.3 0.2 5.0", "1.0 1.0"], "{deck}"); + let single: Vec<&str> = deck.match_indices("*CONTACT_AUTOMATIC_SINGLE_SURFACE").map(|(i, _)| &deck[i..]).collect(); + assert_eq!(single.len(), 2, "{deck}"); + let cards = |text: &str| -> Vec { + text.lines().skip(1).filter(|l| !l.starts_with('$')).take(3).map(|l| l.split_whitespace().collect::>().join(" ")).collect() + }; + // `everything` (every part), then `fold` (the upper face's segment set). + assert_eq!(cards(single[0]), ["0 5", "0.1 0.1 0.0", "1.0 1.0"], "{deck}"); + assert_eq!(cards(single[1]), ["4 0", "0.0 0.0 0.0", "0.5 0.5"], "{deck}"); + for set in ["lower", "upper"] { + assert!(deck.contains(&format!("*SET_SEGMENT_TITLE\n{set}\n")), "{set}\n{deck}"); + } +} + +#[test] +fn nastran_refuses_contact_by_name() { + match nastran::write(&model()) { + Err(Error::Unsupported { what, .. }) => { + for name in ["everything", "fold", "stack"] { + assert!(what.contains(&format!("contact `{name}`")), "{what}"); + } + } + other => panic!("{other:?}"), + } +} diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index d8cbb61..3515eab 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -415,9 +415,29 @@ pub static CATALOGUE: &[Capability] = &[ &["model.coupling.edit", "model.coupling.delete"], ), soon("interaction.tie", "Tie", Interactions, "Constraints", "M7", "Tie non-matching meshes together"), - soon("interaction.contact.general", "General Contact", Interactions, "Contact", "M7", "Contact between everything that touches, with friction"), - soon("interaction.contact.pair", "Contact Pair", Interactions, "Contact", "M7", "Contact between two surfaces"), - soon("interaction.contact.find", "Find Contacts", Interactions, "Contact", "M7", "Detect contact pairs automatically"), + ready_more( + "interaction.contact.general", + "General Contact", + Interactions, + "Contact", + "Contact between everything that may touch, or self contact of a surface, with friction (explicit dynamics)", + &["interaction.contact.edit", "interaction.contact.delete"], + ), + ready_more( + "interaction.contact.pair", + "Contact Pair", + Interactions, + "Contact", + "Contact between two surfaces, with friction (explicit dynamics)", + &["interaction.contact.edit", "interaction.contact.delete"], + ), + ready( + "interaction.contact.find", + "Find Contacts", + Interactions, + "Contact", + "Find the parts of the mesh that touch, and make contact pairs of them", + ), soon("interaction.rigid_wall", "Rigid Wall", Interactions, "Contact", "M10", "Rigid walls for explicit dynamics"), soon("interaction.bolt", "Bolt Load", Interactions, "Connectors", "M7", "Bolt pretension"), ready_more( diff --git a/crates/engine/src/cmd/contact.rs b/crates/engine/src/cmd/contact.rs new file mode 100644 index 0000000..2a31216 --- /dev/null +++ b/crates/engine/src/cmd/contact.rs @@ -0,0 +1,358 @@ +//! `interaction.contact.*`: contact definitions (pairs of surfaces, general contact) and the +//! search for surfaces that touch. +//! +//! A contact's surfaces are sets of element faces ([`fieldcraft_model::contact`]). The commands +//! take the friction, penalty scale and initial-penetration choice of the model's form; `edit` +//! replaces a definition with the model's JSON form, `delete` removes it. +//! +//! `find` looks for the bodies of the mesh that touch: the solid elements are split into +//! connected parts (elements sharing a node belong to the same part), and a node of one part's +//! exterior faces within `tolerance` of an exterior face of another part makes the two parts +//! touch there. Each touching pair is reported with the faces of each side that are close to the +//! other; with `create`, each becomes two face sets and a contact pair. + +use std::collections::BTreeMap; + +use fieldcraft_fem::contact::{Grid, bounds, closest_point}; +use fieldcraft_model::{Contact, ContactSurfaces, ElementFace, Friction, InitialPenetration, Model, Set}; +use serde_json::{Value, json}; + +use super::{name, opt}; +use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; + +/// The options every contact takes, from the parameters. +fn contact(p: &Value, surfaces: ContactSurfaces) -> Result { + Ok(Contact { + surfaces, + friction: opt::(p, "friction")?.unwrap_or(Friction::NONE), + penalty_scale: opt::(p, "penalty_scale")?.unwrap_or(1.0), + initial_penetration: opt::(p, "initial_penetration")?.unwrap_or_default(), + }) +} + +/// A required name-valued parameter. +fn text<'a>(p: &'a Value, key: &str) -> Result<&'a str> { + p.get(key).and_then(Value::as_str).ok_or_else(|| EngineError::Params(format!("`{key}` must be a set name"))) +} + +pub const PAIR: CommandSpec = CommandSpec { + id: "interaction.contact.pair", + 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)"# + ), + kind: CommandKind::Model, + run: |s: &mut Session, p: &Value| { + let n = name(p)?; + let surfaces = ContactSurfaces::Pair { main: text(p, "main")?.to_string(), secondary: text(p, "secondary")?.to_string() }; + s.model.add_contact(n, contact(p, surfaces)?)?; + Ok(json!({"name": n})) + }, +}; + +pub const GENERAL: CommandSpec = CommandSpec { + id: "interaction.contact.general", + 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)"# + ), + kind: CommandKind::Model, + run: |s: &mut Session, p: &Value| { + let n = name(p)?; + let surface = match p.get("surface") { + None | Some(Value::Null) => None, + Some(_) => Some(text(p, "surface")?.to_string()), + }; + s.model.add_contact(n, contact(p, ContactSurfaces::General { surface })?)?; + Ok(json!({"name": n})) + }, +}; + +pub const EDIT: CommandSpec = CommandSpec { + id: "interaction.contact.edit", + label: "Edit Contact", + params: r#"{"name": "impact", "surfaces": {"type": "pair", "main": "plate", "secondary": "ball"} | {"type": "general", "surface": "skin" (optional)}, "friction": {…}, "penalty_scale": 1, "initial_penetration": "adjust|refuse"} — replaces the contact of that name (the model's form; see interaction.contact.pair)"#, + kind: CommandKind::Model, + run: |s: &mut Session, p: &Value| { + let n = name(p)?; + let mut obj = p.as_object().cloned().ok_or_else(|| EngineError::Params("the parameters of a contact are a JSON object".into()))?; + obj.remove("name"); + let c: Contact = serde_json::from_value(Value::Object(obj)).map_err(|e| EngineError::Params(format!("contact: {e}")))?; + s.model.replace_contact(n, c)?; + Ok(json!({"name": n})) + }, +}; + +pub const DELETE: CommandSpec = CommandSpec { + id: "interaction.contact.delete", + label: "Delete Contact", + params: r#"{"name": "impact"}"#, + kind: CommandKind::Model, + run: |s: &mut Session, p: &Value| { + s.model.remove_contact(name(p)?)?; + Ok(Value::Null) + }, +}; + +pub const FIND: CommandSpec = CommandSpec { + id: "interaction.contact.find", + label: "Find Contacts", + params: concat!( + r#"{"tolerance": 0.1 (distance within which surfaces touch; default a tenth of the mean exterior face size), "create": false (true: add two face sets and a contact pair per touching pair of parts, named find-N-main, find-N-secondary and find-N), "#, + r#""friction", "penalty_scale", "initial_penetration" (of the created pairs, as interaction.contact.pair)} — the parts of the mesh (solid elements connected through their nodes) whose exterior faces come within the tolerance of each other"# + ), + kind: CommandKind::Model, + run: find, +}; + +/// Most touching part pairs `find` reports. +pub const MAX_FOUND: usize = fieldcraft_model::contact::MAX_CONTACTS; + +/// Most point–face checks `find` makes. +pub const MAX_CHECKS: usize = 200_000_000; + +/// The connected parts of the solid elements: per element label, its part (numbered from 0 in +/// the order of their lowest element label). +fn parts(model: &Model) -> BTreeMap { + // Union–find over the node labels of the solid elements. + let solids: Vec<(u32, &[u32])> = + model.elements().iter().filter(|(_, e)| e.kind.dimension() == 3).map(|(l, e)| (*l, e.nodes.as_slice())).collect(); + let mut index: BTreeMap = BTreeMap::new(); + for (_, nodes) in &solids { + for n in *nodes { + let next = index.len(); + index.entry(*n).or_insert(next); + } + } + let mut parent: Vec = (0..index.len()).collect(); + fn root(parent: &mut [usize], mut i: usize) -> usize { + while let Some(&p) = parent.get(i) { + if p == i { + break; + } + let grand = parent.get(p).copied().unwrap_or(p); + if let Some(slot) = parent.get_mut(i) { + *slot = grand; + } + i = grand; + } + i + } + for (_, nodes) in &solids { + let mut it = nodes.iter().filter_map(|n| index.get(n).copied()); + let Some(first) = it.next() else { continue }; + for other in it { + let (a, b) = (root(&mut parent, first), root(&mut parent, other)); + if a != b { + let (lo, hi) = (a.min(b), a.max(b)); + if let Some(slot) = parent.get_mut(hi) { + *slot = lo; + } + } + } + } + // Number the parts in element order. + let mut numbers: BTreeMap = BTreeMap::new(); + let mut out = BTreeMap::new(); + for (label, nodes) in &solids { + let Some(r) = nodes.first().and_then(|n| index.get(n).copied()).map(|i| root(&mut parent, i)) else { continue }; + let next = numbers.len(); + out.insert(*label, *numbers.entry(r).or_insert(next)); + } + out +} + +fn find(s: &mut Session, p: &Value) -> Result { + let model = &s.model; + let faces = model.exterior_faces(); + let part_of = parts(model); + // Each exterior face's corner coordinates and part. + let mut face_xyz: Vec> = Vec::with_capacity(faces.len()); + for f in &faces { + let labels = model.face_nodes(*f).unwrap_or_default(); + let corners = if labels.len() == 3 || labels.len() == 6 { 3 } else { 4 }; + face_xyz.push(labels.iter().take(corners).filter_map(|l| model.node(*l)).map(|n| n.xyz).collect()); + } + let size = { + let sizes: Vec = face_xyz + .iter() + .filter_map(|x| { + let b = bounds(x, 0.0); + let d = (0..3).map(|k| b[1][k] - b[0][k]).fold(0.0, f64::max); + d.is_finite().then_some(d) + }) + .collect(); + if sizes.is_empty() { 0.0 } else { sizes.iter().sum::() / sizes.len() as f64 } + }; + let tolerance = opt::(p, "tolerance")?.unwrap_or(0.1 * size); + if !(tolerance >= 0.0 && tolerance.is_finite()) { + return Err(EngineError::Params(format!("`tolerance` must be a finite distance ≥ 0, got {tolerance}"))); + } + let create = opt::(p, "create")?.unwrap_or(false); + let options = contact(p, ContactSurfaces::General { surface: None })?; + options.validate()?; + + let boxes: Vec<_> = face_xyz.iter().map(|x| bounds(x, tolerance)).collect(); + let grid = Grid::new(&boxes, 4 * faces.len() + 64); + let part = |c: usize| faces.get(c).and_then(|f| part_of.get(&f.element).copied()); + // Each face's centroid and outward normal there. + let centroid = |x: &[[f64; 3]]| -> [f64; 3] { [0, 1, 2].map(|k| x.iter().map(|p| p[k]).sum::() / x.len().max(1) as f64) }; + let normals: Vec> = face_xyz.iter().map(|x| closest_point(x, centroid(x)).map(|q| q.normal)).collect(); + // (lower part, higher part) → the faces of each side that touch the other: a face whose + // centroid or a corner is within the tolerance of a face of the other part facing it + // (normals more than 60° apart from parallel). + let mut touching: BTreeMap<(usize, usize), [std::collections::BTreeSet; 2]> = BTreeMap::new(); + let mut checks = 0usize; + for (c, xyz) in face_xyz.iter().enumerate() { + let (Some(k), Some(normal)) = (part(c), normals.get(c).copied().flatten()) else { continue }; + for point in xyz.iter().copied().chain(std::iter::once(centroid(xyz))) { + let candidates = grid.candidates(point); + checks += candidates.len(); + if checks > MAX_CHECKS { + return Err(EngineError::Params(format!( + "the tolerance {tolerance:e} puts too many faces near each other (more than {MAX_CHECKS} checks); lower it" + ))); + } + for &d in candidates { + let d = d as usize; + let Some(other) = part(d).filter(|o| *o != k) else { continue }; + let facing = normals.get(d).copied().flatten().is_some_and(|n| (0..3).map(|i| n[i] * normal[i]).sum::() < -0.5); + let near = facing && face_xyz.get(d).and_then(|x| closest_point(x, point)).is_some_and(|q| q.distance <= tolerance); + if near { + let key = (k.min(other), k.max(other)); + touching.entry(key).or_default()[usize::from(k != key.0)].insert(c); + if touching.len() > MAX_FOUND { + return Err(EngineError::Params(format!("more than {MAX_FOUND} pairs of parts touch; lower the tolerance"))); + } + } + } + } + } + let face_list = |set: &std::collections::BTreeSet| -> Vec { set.iter().filter_map(|c| faces.get(*c).copied()).collect() }; + let mut found = Vec::new(); + let mut created = Vec::new(); + for (n, ((a, b), sides)) in touching.iter().enumerate() { + let (main, secondary) = (face_list(&sides[0]), face_list(&sides[1])); + let mut entry = json!({"parts": [a, b], "faces": [main.len(), secondary.len()]}); + if create { + let names = [format!("find-{}-main", n + 1), format!("find-{}-secondary", n + 1), format!("find-{}", n + 1)]; + created.push((names, main, secondary)); + entry["contact"] = json!(format!("find-{}", n + 1)); + } + found.push(entry); + } + let part_count = part_of.values().max().map_or(0, |m| m + 1); + for ([main_name, secondary_name, contact_name], main, secondary) in created { + s.model.add_set(&main_name, Set { faces: main, ..Set::default() })?; + s.model.add_set(&secondary_name, Set { faces: secondary, ..Set::default() })?; + let surfaces = ContactSurfaces::Pair { main: main_name, secondary: secondary_name }; + s.model.add_contact(&contact_name, Contact { surfaces, ..options.clone() })?; + } + Ok(json!({"parts": part_count, "exterior_faces": faces.len(), "tolerance": tolerance, "touching": found})) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Two unit cubes of 2×2×2 hex8 each, the second resting on the first (separate nodes), and + /// a third cube far away. + fn stacked(s: &mut Session) { + let mut next_node = 1u32; + let mut next_element = 1u32; + for (z0, x0) in [(0.0, 0.0), (1.0, 0.0), (0.0, 10.0)] { + let base = next_node; + for k in 0..3 { + for j in 0..3 { + for i in 0..3 { + let xyz = [x0 + 0.5 * i as f64, 0.5 * j as f64, z0 + 0.5 * k as f64]; + s.execute("mesh.node.add", &json!({"label": next_node, "xyz": xyz})).unwrap(); + next_node += 1; + } + } + } + let n = |i: u32, j: u32, k: u32| base + i + 3 * j + 9 * k; + for k in 0..2 { + for j in 0..2 { + for i in 0..2 { + let nodes = [ + n(i, j, k), + n(i + 1, j, k), + n(i + 1, j + 1, k), + n(i, j + 1, k), + n(i, j, k + 1), + n(i + 1, j, k + 1), + n(i + 1, j + 1, k + 1), + n(i, j + 1, k + 1), + ]; + s.execute("mesh.element.add", &json!({"label": next_element, "kind": "hex8", "nodes": nodes})).unwrap(); + next_element += 1; + } + } + } + } + } + + #[test] + fn contacts_are_created_edited_and_deleted() { + let mut s = Session::default(); + stacked(&mut s); + s.execute("model.set.create", &json!({"name": "a", "faces": [{"element": 5, "face": 2}]})).unwrap(); + s.execute("model.set.create", &json!({"name": "b", "faces": [{"element": 9, "face": 1}]})).unwrap(); + s.execute("interaction.contact.pair", &json!({"name": "ab", "main": "a", "secondary": "b", "friction": {"static": 0.2}})).unwrap(); + let c = s.model.contact("ab").unwrap(); + assert_eq!(c.friction, Friction::coulomb(0.2)); + assert!(s.execute("interaction.contact.pair", &json!({"name": "ab", "main": "a", "secondary": "b"})).is_err(), "taken"); + assert!(s.execute("interaction.contact.pair", &json!({"name": "aa", "main": "a", "secondary": "a"})).is_err(), "self pair"); + assert!(s.execute("interaction.contact.pair", &json!({"name": "x", "main": "a", "secondary": "nope"})).is_err()); + assert!(s.execute("interaction.contact.pair", &json!({"name": "x", "main": "a", "secondary": "b", "penalty_scale": 2})).is_err()); + assert!(s.execute("interaction.contact.pair", &json!({"name": "x", "main": "a", "secondary": "b", "friction": {"static": -1}})).is_err()); + s.execute("interaction.contact.general", &json!({"name": "all"})).unwrap(); + s.execute("interaction.contact.general", &json!({"name": "self", "surface": "a", "initial_penetration": "refuse"})).unwrap(); + assert_eq!(s.model.contact("self").unwrap().initial_penetration, InitialPenetration::Refuse); + s.execute("interaction.contact.edit", &json!({"name": "all", "surfaces": {"type": "general", "surface": "b"}, "penalty_scale": 0.5})) + .unwrap(); + assert_eq!(s.model.contact("all").unwrap().surfaces, ContactSurfaces::General { surface: Some("b".into()) }); + assert!(s.execute("interaction.contact.edit", &json!({"name": "nope", "surfaces": {"type": "general"}})).is_err()); + assert!(s.execute("model.set.delete", &json!({"name": "a"})).unwrap_err().to_string().contains("contact `ab`")); + s.execute("model.set.create", &json!({"name": "n", "nodes": [1, 2]})).unwrap(); + s.execute("interaction.contact.general", &json!({"name": "bad", "surface": "n"})).unwrap(); + let check = s.execute("model.check", &json!({})).unwrap(); + assert!(check.to_string().contains("set `n` has no element faces"), "{check}"); + for n in ["ab", "all", "self", "bad"] { + s.execute("interaction.contact.delete", &json!({"name": n})).unwrap(); + } + assert!(s.model.contacts().is_empty()); + } + + #[test] + fn find_reports_and_creates_the_touching_parts() { + let mut s = Session::default(); + stacked(&mut s); + let r = s.execute("interaction.contact.find", &json!({})).unwrap(); + assert_eq!(r["parts"], 3, "{r}"); + assert_eq!(r["exterior_faces"], 3 * 24, "{r}"); + // The two stacked cubes touch over four faces each; the far one touches nothing. + assert_eq!(r["touching"], json!([{"parts": [0, 1], "faces": [4, 4]}]), "{r}"); + assert!(s.model.contacts().is_empty(), "nothing created without `create`"); + let r = s.execute("interaction.contact.find", &json!({"create": true, "friction": {"static": 0.1}})).unwrap(); + assert_eq!(r["touching"][0]["contact"], "find-1"); + let c = s.model.contact("find-1").unwrap(); + assert_eq!(c.friction, Friction::coulomb(0.1)); + let main = s.model.resolve_set("find-1-main").unwrap(); + assert!(main.faces.iter().all(|f| (5..=8).contains(&f.element) && f.face == 2) && main.faces.len() == 4, "{main:?}"); + let secondary = s.model.resolve_set("find-1-secondary").unwrap(); + assert!(secondary.faces.iter().all(|f| (9..=12).contains(&f.element) && f.face == 1) && secondary.faces.len() == 4, "{secondary:?}"); + // A tight tolerance on a gap finds nothing; a second create reuses no names. + assert!(s.execute("interaction.contact.find", &json!({"create": true})).is_err(), "the names are taken"); + for bad in [json!({"tolerance": -1}), json!({"tolerance": "x"}), json!({"tolerance": f64::MAX, "create": "yes"})] { + assert!(s.execute("interaction.contact.find", &bad).is_err(), "{bad}"); + } + let empty = s.execute("interaction.contact.find", &json!({"tolerance": 0})).unwrap(); + assert_eq!(empty["touching"].as_array().map(Vec::len), Some(1), "touching faces are at distance 0: {empty}"); + } +} diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 7f370b8..d6bc36e 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -1,6 +1,7 @@ //! The command registry. Each area has its own module; [`ALL`] lists every command. mod block; +mod contact; pub(crate) mod deck; pub(crate) mod define; mod edit; @@ -67,6 +68,11 @@ pub static ALL: &[CommandSpec] = &[ define::COUPLING_CREATE, define::COUPLING_EDIT, define::COUPLING_DELETE, + contact::PAIR, + contact::GENERAL, + contact::EDIT, + contact::DELETE, + contact::FIND, define::MASS_CREATE, define::MASS_EDIT, define::MASS_DELETE, diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 0ebb31b..e02656d 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -477,6 +477,18 @@ mod tests { json!({"name": "A", "procedure": {"type": "heat_transfer", "transient": {"end_time": -1, "increment": 0}}}), json!({"name": "A", "procedure": {"type": "coupled_temperature_displacement", "extra": 1}}), json!({"step": "Heat", "frame": u64::MAX, "nodes": [1]}), + // Contact. + json!({"name": "C", "main": 7, "secondary": null}), + json!({"name": "C", "main": "end", "secondary": "root", "friction": {"static": 1e308, "kinetic": -1, "decay": -1e308}}), + json!({"name": "C", "main": "end", "secondary": "end"}), + json!({"name": "C", "main": "x".repeat(10_000), "secondary": "root", "penalty_scale": 1e308}), + json!({"name": "C", "surface": "root", "penalty_scale": 0, "initial_penetration": "ignore"}), + json!({"name": "C", "surface": 3}), + json!({"name": "C", "surfaces": {"type": "pair", "main": "end"}}), + json!({"name": "C", "surfaces": {"type": "general", "surface": "nope"}, "friction": "sticky"}), + json!({"tolerance": 1e300, "create": true, "friction": {"static": 0.1}}), + json!({"tolerance": -0.0, "create": true}), + json!({"tolerance": 1e-300, "create": 1}), json!({"geometry": {"body": "Bar"}, "elements": [1, 99], "faces": [{"element": 1, "face": 200}]}), json!({"geometry": {"body": "Bar", "vertex": "v"}, "nodes": [1], "edges": [{"element": 1, "edge": 0}]}), ]; diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index e70677c..897eb47 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -196,6 +196,12 @@ impl System { if let Some(name) = model.couplings().keys().next() { return Err(Error::NotSupported(format!("couplings (coupling `{name}`) in explicit dynamics"))); } + if let Some(name) = model.contacts().keys().next() { + return Err(Error::NotSupported(format!("contact (`{name}`) in explicit dynamics"))); + } + if let Some(w) = settings.rigid_walls.iter().find(|w| !w.friction.is_frictionless()) { + return Err(Error::NotSupported(format!("friction on rigid walls (the wall on set `{}`)", w.set))); + } // Elements by label, with their section's material. let mut materials: Vec = Vec::new(); let mut chosen: Vec<(u32, ElementKind, Vec, usize, Option)> = Vec::new(); diff --git a/crates/explicit/tests/large_deformation.rs b/crates/explicit/tests/large_deformation.rs index bde845d..5901806 100644 --- a/crates/explicit/tests/large_deformation.rs +++ b/crates/explicit/tests/large_deformation.rs @@ -318,7 +318,7 @@ fn taylor_bar_recovers_its_yield_stress() { step.output_interval = Some(end / 4.0); step.history_interval = Some(end / 400.0); step.initial_velocities = vec![InitialVelocity { set: "nodes".into(), velocity: [0.0, 0.0, -speed] }]; - step.rigid_walls = vec![RigidWall { set: "nodes".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0] }]; + step.rigid_walls = vec![RigidWall { set: "nodes".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0], friction: Friction::NONE }]; step.hourglass = Hourglass { kind: HourglassKind::Stiffness, coefficient: 0.1 }; let mut m = block([4, 4, 30], [width, width, length], ElementKind::Hex8); node_set(&mut m, "nodes", |_| true); @@ -359,7 +359,7 @@ fn taylor_with(material: Material) -> (Model, ExplicitResult) { let mut step = ExplicitDynamic::new(end); step.output_interval = Some(end / 2.0); step.initial_velocities = vec![InitialVelocity { set: "nodes".into(), velocity: [0.0, 0.0, -300_000.0] }]; - step.rigid_walls = vec![RigidWall { set: "nodes".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0] }]; + step.rigid_walls = vec![RigidWall { set: "nodes".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0], friction: Friction::NONE }]; step.hourglass = Hourglass { kind: HourglassKind::Stiffness, coefficient: 0.1 }; let mut m = block([2, 2, 12], [7.6, 7.6, 25.4], ElementKind::Hex8); node_set(&mut m, "nodes", |_| true); diff --git a/crates/explicit/tests/verification.rs b/crates/explicit/tests/verification.rs index d058cde..c0e3a7a 100644 --- a/crates/explicit/tests/verification.rs +++ b/crates/explicit/tests/verification.rs @@ -235,7 +235,7 @@ fn bar_rebounds_from_a_rigid_wall() { step.history_interval = Some(end / 3000.0); step.output_interval = Some(end); step.initial_velocities = vec![InitialVelocity { set: "all_nodes".into(), velocity: [-v0, 0.0, 0.0] }]; - step.rigid_walls = vec![RigidWall { set: "start".into(), point: [0.0; 3], normal: [1.0, 0.0, 0.0] }]; + step.rigid_walls = vec![RigidWall { set: "start".into(), point: [0.0; 3], normal: [1.0, 0.0, 0.0], friction: Friction::NONE }]; step.bulk_viscosity.linear = linear; let mut bar = Bar::new(kind, 200, ExplicitDynamic::new(end)); let nodes: Vec = bar.model.nodes().keys().copied().collect(); @@ -347,7 +347,7 @@ fn caps_and_unsupported_models_are_errors() { assert!(solve_explicit(&light.model, None, &Limits::default()).is_err()); // A node behind its wall at the start. let mut behind = ExplicitDynamic::new(end); - behind.rigid_walls = vec![RigidWall { set: "end".into(), point: [LENGTH + 1.0, 0.0, 0.0], normal: [-1.0, 0.0, 0.0] }]; + behind.rigid_walls = vec![RigidWall { set: "end".into(), point: [LENGTH + 1.0, 0.0, 0.0], normal: [-1.0, 0.0, 0.0], friction: Friction::NONE }]; let mut wall = Bar::new(ElementKind::Hex8, 4, ExplicitDynamic::new(end)); wall.model.replace_step(Step::new("Dynamic", Procedure::ExplicitDynamic(behind.clone()))).unwrap(); assert!(solve_explicit(&wall.model, None, &Limits::default()).is_ok(), "the wall is beyond the end node, which is on its allowed side"); diff --git a/crates/fem/src/contact.rs b/crates/fem/src/contact.rs new file mode 100644 index 0000000..714cafa --- /dev/null +++ b/crates/fem/src/contact.rs @@ -0,0 +1,398 @@ +//! Contact geometry shared by the solvers: the closest point of an element face to a point (the +//! local search) and a bucket grid of face boxes (the global search). +//! +//! **Closest point.** A face is a flat triangle or a bilinear quadrilateral given by its corner +//! nodes, numbered so that their right-hand normal points out of the element. On a triangle the +//! closest point is found exactly by the region tests of C. Ericson, *Real-Time Collision +//! Detection*, Morgan Kaufmann 2005, §5.1.5. On a bilinear quadrilateral x(ξ, η) it is the +//! minimum of ½|x(ξ, η) − p|² over [−1, 1]², found by Newton's method on the stationarity +//! conditions x,ξ · (x − p) = x,η · (x − p) = 0 with the exact Hessian (the bilinear patch's +//! only second derivative is x,ξη), projected onto the square after each step: the +//! closest-point projection of node-to-segment contact (J. O. Hallquist, G. L. Goudreau, D. J. +//! Benson, *Sliding interfaces with contact-impact in large-scale Lagrangian computations*, +//! CMAME 51:107–137, 1985; P. Wriggers, *Computational Contact Mechanics*, 2nd ed., Springer +//! 2006, §4.1). The signed gap is measured along the face's unit normal at the closest point; the +//! lateral offset (how far the point lies beside the face rather than above or below it) tells a +//! point over the face from one past its edge. +//! +//! **Global search.** Each face's bounding box, grown by a margin, goes into every cell of a +//! uniform grid it overlaps (a bucket sort: D. J. Benson, J. O. Hallquist, *A single surface +//! contact algorithm for the post-buckling analysis of shell structures*, CMAME 78:141–163, +//! 1990); a point's candidates are the faces in its cell. The cell is twice the mean box size, +//! grown until the grid and its entries fit their caps, so the search is linear in the number of +//! faces and its memory bounded whatever the input. Entries are listed in face order, so the +//! candidates (and everything decided from them) do not depend on how the work is split. + +pub use crate::frame::Vec3; +use crate::frame::{cross, dot, norm, sub}; + +/// Newton iterations of the bilinear closest-point search. +const MAX_NEWTON: usize = 30; + +/// The closest point of a face to a point. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Projection { + /// Face coordinates of the closest point: (ξ, η) ∈ [−1, 1]² on a quadrilateral, the + /// barycentric coordinates of the second and third nodes on a triangle. + pub xi: [f64; 2], + /// The face's shape functions at the closest point, per node (0 past the face's nodes). + pub weights: [f64; 4], + pub point: Vec3, + /// The unit outward normal at the closest point. + pub normal: Vec3, + /// (p − point) · normal: negative when the point is behind the face. + pub gap: f64, + /// |p − point|. + pub distance: f64, + /// The part of p − point across the normal: 0 when the point is over or under the face, its + /// distance from the face's edge when it is beside it. + pub lateral: f64, +} + +/// The closest point of a face (three or four corner nodes, in its node order) to `p`; `None` +/// for another number of nodes, non-finite coordinates or a face collapsed at that point. +pub fn closest_point(face: &[Vec3], p: Vec3) -> Option { + if face.iter().chain(std::iter::once(&p)).flatten().any(|c| !c.is_finite()) { + return None; + } + let (xi, weights, normal) = match face { + [a, b, c] => triangle(*a, *b, *c, p)?, + [a, b, c, d] => quadrilateral([*a, *b, *c, *d], p)?, + _ => return None, + }; + let mut point = [0.0; 3]; + for (x, w) in face.iter().zip(weights) { + for k in 0..3 { + point[k] += w * x[k]; + } + } + let d = sub(p, point); + let gap = dot(d, normal); + let distance = norm(d); + let lateral = (distance * distance - gap * gap).max(0.0).sqrt(); + Some(Projection { xi, weights, point, normal, gap, distance, lateral }) +} + +/// The unit vector along `v`, unless it is too short against `scale` (a squared length). +fn unit(v: Vec3, scale: f64) -> Option { + let n = norm(v); + (n > 1e-12 * scale && n.is_finite()).then(|| v.map(|c| c / n)) +} + +/// Ericson's closest point on triangle abc: (barycentric v, w), weights, unit normal. +fn triangle(a: Vec3, b: Vec3, c: Vec3, p: Vec3) -> Option<([f64; 2], [f64; 4], Vec3)> { + let (ab, ac, ap) = (sub(b, a), sub(c, a), sub(p, a)); + let normal = unit(cross(ab, ac), dot(ab, ab).max(dot(ac, ac)))?; + let (d1, d2) = (dot(ab, ap), dot(ac, ap)); + let bp = sub(p, b); + let (d3, d4) = (dot(ab, bp), dot(ac, bp)); + let cp = sub(p, c); + let (d5, d6) = (dot(ab, cp), dot(ac, cp)); + let (va, vb, vc) = (d3 * d6 - d5 * d4, d5 * d2 - d1 * d6, d1 * d4 - d3 * d2); + let (v, w) = if d1 <= 0.0 && d2 <= 0.0 { + (0.0, 0.0) + } else if d3 >= 0.0 && d4 <= d3 { + (1.0, 0.0) + } else if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 { + (d1 / (d1 - d3), 0.0) + } else if d6 >= 0.0 && d5 <= d6 { + (0.0, 1.0) + } else if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 { + (0.0, d2 / (d2 - d6)) + } else if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 { + let t = (d4 - d3) / ((d4 - d3) + (d5 - d6)); + (1.0 - t, t) + } else { + let s = va + vb + vc; + (vb / s, vc / s) + }; + if !(v.is_finite() && w.is_finite()) { + return None; + } + Some(([v, w], [1.0 - v - w, v, w, 0.0], normal)) +} + +/// The bilinear shape functions and their derivatives at (ξ, η), nodes at (−1, −1), (1, −1), +/// (1, 1), (−1, 1). +fn bilinear(s: f64, t: f64) -> ([f64; 4], [f64; 4], [f64; 4]) { + let n = [(1.0 - s) * (1.0 - t), (1.0 + s) * (1.0 - t), (1.0 + s) * (1.0 + t), (1.0 - s) * (1.0 + t)].map(|v| 0.25 * v); + let ds = [-(1.0 - t), 1.0 - t, 1.0 + t, -(1.0 + t)].map(|v| 0.25 * v); + let dt = [-(1.0 - s), -(1.0 + s), 1.0 + s, 1.0 - s].map(|v| 0.25 * v); + (n, ds, dt) +} + +fn combine(x: &[Vec3; 4], w: &[f64; 4]) -> Vec3 { + let mut out = [0.0; 3]; + for (xi, wi) in x.iter().zip(w) { + for k in 0..3 { + out[k] += wi * xi[k]; + } + } + out +} + +/// The closest point of a bilinear quadrilateral by projected Newton iterations. +fn quadrilateral(x: [Vec3; 4], p: Vec3) -> Option<([f64; 2], [f64; 4], Vec3)> { + let size = dot(sub(x[2], x[0]), sub(x[2], x[0])).max(dot(sub(x[3], x[1]), sub(x[3], x[1]))); + if !(size > 0.0) { + return None; + } + // x,ξη: constant over the patch. + let twist = [0, 1, 2].map(|k| 0.25 * (x[0][k] - x[1][k] + x[2][k] - x[3][k])); + let (mut s, mut t) = (0.0f64, 0.0f64); + for _ in 0..MAX_NEWTON { + let (n, ds, dt) = bilinear(s, t); + let r = sub(combine(&x, &n), p); + let (a, b) = (combine(&x, &ds), combine(&x, &dt)); + let g = [dot(a, r), dot(b, r)]; + let (haa, hbb) = (dot(a, a), dot(b, b)); + let mut hab = dot(a, b) + dot(r, twist); + // Far from the patch the exact Hessian may be indefinite: Gauss–Newton then. + if haa * hbb - hab * hab <= 1e-14 * haa * hbb { + hab = dot(a, b); + } + let det = haa * hbb - hab * hab; + if !(det > 0.0) { + return None; + } + let step = [-(hbb * g[0] - hab * g[1]) / det, -(haa * g[1] - hab * g[0]) / det]; + let (s1, t1) = ((s + step[0]).clamp(-1.0, 1.0), (t + step[1]).clamp(-1.0, 1.0)); + let moved = (s1 - s).abs().max((t1 - t).abs()); + (s, t) = (s1, t1); + if !moved.is_finite() { + return None; + } + if moved < 1e-13 { + break; + } + } + let (n, ds, dt) = bilinear(s, t); + let normal = unit(cross(combine(&x, &ds), combine(&x, &dt)), size)?; + Some(([s, t], n, normal)) +} + +/// A box: its lowest and highest corners. +pub type Aabb = [Vec3; 2]; + +/// The bounding box of a face's nodes, grown by `margin` on every side. +pub fn bounds(face: &[Vec3], margin: f64) -> Aabb { + let mut b = [[f64::INFINITY; 3], [f64::NEG_INFINITY; 3]]; + for x in face { + for k in 0..3 { + b[0][k] = b[0][k].min(x[k] - margin); + b[1][k] = b[1][k].max(x[k] + margin); + } + } + b +} + +/// A uniform grid of boxes (a bucket sort), for the global contact search. +#[derive(Debug, Clone, Default)] +pub struct Grid { + lowest: Vec3, + cell: f64, + dims: [usize; 3], + /// Per cell, the range of its entries in `items` (compressed rows). + start: Vec, + items: Vec, +} + +/// Most cells of a grid (whatever `max_cells` says). +pub const MAX_GRID_CELLS: usize = 1 << 24; + +/// Most entries of a grid (box × cell overlaps). +pub const MAX_GRID_ITEMS: usize = 1 << 27; + +impl Grid { + /// The grid of `boxes` (their indices are the entries; boxes that are not finite are left + /// out), with at most `max_cells` cells (and [`MAX_GRID_CELLS`]). + pub fn new(boxes: &[Aabb], max_cells: usize) -> Grid { + let finite = |b: &Aabb| b.iter().flatten().all(|c| c.is_finite()) && (0..3).all(|k| b[0][k] <= b[1][k]); + let mut lowest = [f64::INFINITY; 3]; + let mut highest = [f64::NEG_INFINITY; 3]; + let (mut count, mut extent) = (0usize, 0.0); + for b in boxes.iter().filter(|b| finite(b)) { + for k in 0..3 { + lowest[k] = lowest[k].min(b[0][k]); + highest[k] = highest[k].max(b[1][k]); + } + extent += (0..3).map(|k| b[1][k] - b[0][k]).fold(0.0, f64::max); + count += 1; + } + if count == 0 { + return Grid::default(); + } + let span = (0..3).map(|k| highest[k] - lowest[k]).fold(0.0, f64::max); + let max_cells = max_cells.clamp(1, MAX_GRID_CELLS); + // Twice the mean box size, but not below a millionth of the span (point-like boxes). + let mut cell = (2.0 * extent / count as f64).max(1e-6 * span).max(f64::MIN_POSITIVE); + for _ in 0..200 { + let dims = [0, 1, 2].map(|k| ((highest[k] - lowest[k]) / cell).floor() + 1.0); + let cells = dims[0] * dims[1] * dims[2]; + if !(cells <= max_cells as f64) { + cell *= (cells / max_cells as f64).cbrt().max(1.01) * 1.01; + continue; + } + let dims = dims.map(|d| d as usize); + let mut grid = Grid { lowest, cell, dims, start: vec![0; dims[0] * dims[1] * dims[2] + 1], items: Vec::new() }; + // Count the entries per cell, then fill (compressed rows, boxes in order). + let mut total = 0usize; + for b in boxes.iter().filter(|b| finite(b)) { + let (lo, hi) = (grid.coordinates(b[0]), grid.coordinates(b[1])); + total = total.saturating_add((0..3).map(|k| hi[k] - lo[k] + 1).product::()); + } + if total > MAX_GRID_ITEMS { + cell *= 2.0; + continue; + } + for b in boxes.iter().filter(|b| finite(b)) { + grid.each_cell(b, |c, start| start[c + 1] += 1); + } + for c in 0..grid.start.len() - 1 { + grid.start[c + 1] += grid.start[c]; + } + let mut next = grid.start.clone(); + let mut items = vec![0u32; total]; + for (i, b) in boxes.iter().enumerate().filter(|(_, b)| finite(b)) { + grid.each_cell(b, |c, _| { + if let (Some(slot), Some(item)) = (next.get_mut(c), u32::try_from(i).ok()) { + if let Some(e) = items.get_mut(*slot) { + *e = item; + } + *slot += 1; + } + }); + } + grid.items = items; + return grid; + } + Grid::default() + } + + /// The cell coordinates of a point inside the grid's span (clamped to it). + fn coordinates(&self, x: Vec3) -> [usize; 3] { + [0, 1, 2].map(|k| { + let c = ((x[k] - self.lowest[k]) / self.cell).floor(); + if c > 0.0 { (c as usize).min(self.dims[k].saturating_sub(1)) } else { 0 } + }) + } + + fn index(&self, c: [usize; 3]) -> usize { + c[0] + self.dims[0] * (c[1] + self.dims[1] * c[2]) + } + + /// Call `f(cell, start)` for every cell a box overlaps. + fn each_cell(&mut self, b: &Aabb, mut f: impl FnMut(usize, &mut Vec)) { + let (lo, hi) = (self.coordinates(b[0]), self.coordinates(b[1])); + for k in lo[2]..=hi[2] { + for j in lo[1]..=hi[1] { + for i in lo[0]..=hi[0] { + let c = self.index([i, j, k]); + f(c, &mut self.start); + } + } + } + } + + /// The boxes that may hold `p`: those listed in its cell, in box order (empty outside the + /// grid). + pub fn candidates(&self, p: Vec3) -> &[u32] { + if self.items.is_empty() || p.iter().any(|c| !c.is_finite()) { + return &[]; + } + let outside = (0..3).any(|k| p[k] < self.lowest[k] || p[k] > self.lowest[k] + self.cell * self.dims[k] as f64); + if outside { + return &[]; + } + let c = self.index(self.coordinates(p)); + match (self.start.get(c), self.start.get(c + 1)) { + (Some(a), Some(b)) => self.items.get(*a..*b).unwrap_or_default(), + _ => &[], + } + } + + /// The number of entries (box × cell overlaps). + pub fn entries(&self) -> usize { + self.items.len() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn close(a: f64, b: f64) -> bool { + (a - b).abs() < 1e-12 + } + + #[test] + fn triangle_projection_covers_the_face_its_edges_and_corners() { + let t = [[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 2.0, 0.0]]; + // Above the inside: the foot of the perpendicular, positive gap along +z. + let p = closest_point(&t, [0.5, 0.25, 0.3]).unwrap(); + assert!(close(p.gap, 0.3) && close(p.lateral, 0.0) && close(p.point[0], 0.5) && close(p.point[1], 0.25)); + assert!(close(p.weights.iter().sum::(), 1.0) && p.normal == [0.0, 0.0, 1.0]); + // Below: a negative gap. + assert!(close(closest_point(&t, [0.5, 0.5, -0.1]).unwrap().gap, -0.1)); + // Past the hypotenuse: the closest point is on it, the offset lateral. + let p = closest_point(&t, [2.0, 2.0, 0.0]).unwrap(); + assert!(close(p.point[0], 1.0) && close(p.point[1], 1.0) && close(p.lateral, 2f64.sqrt()) && close(p.gap, 0.0)); + // Past a corner. + let p = closest_point(&t, [-1.0, -1.0, 1.0]).unwrap(); + assert_eq!(p.weights, [1.0, 0.0, 0.0, 0.0]); + assert!(close(p.distance, 3f64.sqrt())); + // Degenerate and hostile faces give nothing. + assert!(closest_point(&[[0.0; 3], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]], [0.0, 1.0, 0.0]).is_none()); + assert!(closest_point(&t, [f64::NAN, 0.0, 0.0]).is_none()); + assert!(closest_point(&t[..2], [0.0; 3]).is_none()); + } + + #[test] + fn quadrilateral_projection_finds_the_foot_on_a_warped_patch() { + let flat = [[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [2.0, 1.0, 0.0], [0.0, 1.0, 0.0]]; + let p = closest_point(&flat, [1.5, 0.25, -0.2]).unwrap(); + assert!(close(p.xi[0], 0.5) && close(p.xi[1], -0.5) && close(p.gap, -0.2) && close(p.lateral, 0.0)); + assert!(close(p.weights.iter().sum::(), 1.0)); + // A warped (hyperbolic paraboloid) patch: the residual is normal to the surface. + let warped = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.2], [1.0, 1.0, 0.0], [0.0, 1.0, 0.2]]; + let q = [0.3, 0.6, 0.5]; + let p = closest_point(&warped, q).unwrap(); + let (_, ds, dt) = bilinear(p.xi[0], p.xi[1]); + let r = sub(q, p.point); + assert!(dot(r, combine(&warped, &ds)).abs() < 1e-12 && dot(r, combine(&warped, &dt)).abs() < 1e-12, "{p:?}"); + assert!(p.lateral < 1e-9 && p.gap > 0.0); + // Beyond an edge: clamped, with the offset lateral. + let p = closest_point(&flat, [3.0, 0.5, 0.0]).unwrap(); + assert!(close(p.xi[0], 1.0) && close(p.lateral, 1.0)); + // A quadrilateral collapsed to a segment has no normal. + assert!(closest_point(&[[0.0; 3], [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0; 3]], [0.5, 0.5, 0.0]).is_none()); + } + + #[test] + fn the_grid_finds_every_box_holding_a_point_and_stays_bounded() { + let boxes: Vec = (0..50).map(|i| [[i as f64, 0.0, 0.0], [i as f64 + 1.5, 1.0, 1.0]]).collect(); + let grid = Grid::new(&boxes, 1000); + for k in 0..500 { + let p = [k as f64 * 0.1 + 0.01, 0.5, 0.5]; + let found = grid.candidates(p); + for (i, b) in boxes.iter().enumerate() { + let inside = (0..3).all(|d| p[d] >= b[0][d] && p[d] <= b[1][d]); + assert!(!inside || found.contains(&(i as u32)), "point {p:?} misses box {i}"); + } + assert!(found.windows(2).all(|w| w[0] < w[1]), "candidates in box order"); + } + assert!(grid.candidates([-5.0, 0.0, 0.0]).is_empty() && grid.candidates([f64::NAN; 3]).is_empty()); + // One huge box among small ones, a cap of one cell, boxes that are not finite: still + // bounded and correct. + let mut mixed = boxes.clone(); + mixed.push([[-1e6, -1e6, -1e6], [1e6, 1e6, 1e6]]); + mixed.push([[f64::NAN, 0.0, 0.0], [1.0, 1.0, 1.0]]); + for cells in [1, 10, 1 << 30] { + let g = Grid::new(&mixed, cells); + assert!(g.candidates([10.2, 0.5, 0.5]).contains(&50) && g.candidates([10.2, 0.5, 0.5]).contains(&9)); + assert!(!g.candidates([10.2, 0.5, 0.5]).contains(&51)); + assert!(g.entries() <= MAX_GRID_ITEMS); + } + assert_eq!(Grid::new(&[], 10).candidates([0.0; 3]), &[] as &[u32]); + } +} diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index d561a1d..43b4e4d 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -28,6 +28,8 @@ //! - [`beam`]: two-node Timoshenko beams with profile section properties; //! - [`frame`]: local axes of shells and beams and the rotation of their matrices; //! - [`face`]: consistent nodal forces of pressure on element faces; +//! - [`contact`]: contact geometry: the closest point of an element face to a point and the +//! bucket grid of the global search; //! - [`heat`]: heat conduction (conductivity, capacity, sources, surface flux and convection, //! heat flux); thermal strain loads and stresses are in [`element`]; //! - [`uniform`]: uniform-strain (one-point) hexahedron and tetrahedron kinematics, hourglass @@ -54,6 +56,7 @@ pub mod basis; pub mod beam; pub mod constitutive; +pub mod contact; pub mod dense; pub mod dof; pub mod elastoplastic; diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs index e279e33..15ad54b 100644 --- a/crates/model/src/check.rs +++ b/crates/model/src/check.rs @@ -387,6 +387,22 @@ impl Model { } } + // Contacts: their surfaces are element faces of solids. + for (name, c) in &self.contacts { + for set in c.surfaces.sets() { + let faces = &set_of(set).faces; + let solid = |f: &crate::ElementFace| self.elements.get(&f.element).is_some_and(|e| e.kind.dimension() == 3); + if faces.is_empty() { + d.error(format!("contact `{name}`"), format!("set `{set}` has no element faces (a contact surface is a set of faces)")); + } else if !faces.iter().all(solid) { + d.error( + format!("contact `{name}`"), + format!("set `{set}` has faces of elements that are not solids (contact acts on the faces of solid elements so far)"), + ); + } + } + } + // Nodes with rotational degrees of freedom: those of shells and beams, and reference points // (worked out only when a boundary condition or load needs it). let rotational = std::cell::OnceCell::new(); diff --git a/crates/model/src/contact.rs b/crates/model/src/contact.rs new file mode 100644 index 0000000..1cf8e09 --- /dev/null +++ b/crates/model/src/contact.rs @@ -0,0 +1,361 @@ +//! 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 +//! 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); +//! - **general** contact: every face against every other face of one surface (single-surface or +//! self contact, for folding and buckling), or of every exterior face of the model's solid +//! elements when no surface is named. +//! +//! **Friction** is Coulomb's law with a coefficient that falls from its static to its kinetic +//! value with the slip rate γ̇: μ = μ_k + (μ_s − μ_k) e^{−d γ̇} (the exponential decay of +//! E. Rabinowicz's observations, as LS-DYNA's and Abaqus's exponential decay friction write it). +//! The same friction applies to the planar rigid walls of an explicit step +//! ([`crate::RigidWall`]). +//! +//! **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. +//! +//! **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 +//! the analysis is refused naming the first node, its depth and the face it is behind +//! (`refuse`). + +use std::collections::BTreeMap; + +use serde::{Deserialize, Serialize}; + +use crate::{Error, Model, Result, context, insert_new, remove_existing, replace_existing, require}; + +/// Most contact definitions in a model. +pub const MAX_CONTACTS: usize = 1000; + +/// Largest friction coefficient (a coefficient of 10 already locks every practical interface). +pub const MAX_FRICTION: f64 = 10.0; + +/// Largest exponential decay coefficient (time per length: 1/d is the slip rate at which the +/// coefficient has gone 63% of the way from static to kinetic). +pub const MAX_DECAY: f64 = 1e15; + +/// Penalty scale factors outside this range are refused (above 1 the contact springs would +/// shorten the stable time increment; see the module docs). +pub const PENALTY_SCALE_RANGE: (f64, f64) = (1e-6, 1.0); + +/// Coulomb friction with an exponential decay from the static to the kinetic coefficient. +#[derive(Debug, Clone, Copy, PartialEq, Default, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Friction { + /// The static coefficient μ_s (0: frictionless). + #[serde(rename = "static", default)] + pub static_coefficient: f64, + /// The kinetic coefficient μ_k (by default the static one: no decay). + #[serde(default, skip_serializing_if = "Option::is_none")] + pub kinetic: Option, + /// The decay coefficient d of μ = μ_k + (μ_s − μ_k) e^{−d γ̇}, per slip rate (time per length). + #[serde(default, skip_serializing_if = "is_zero")] + pub decay: f64, +} + +fn is_zero(v: &f64) -> bool { + *v == 0.0 +} + +impl Friction { + /// No friction. + pub const NONE: Friction = Friction { static_coefficient: 0.0, kinetic: None, decay: 0.0 }; + + /// A constant coefficient. + pub fn coulomb(mu: f64) -> Friction { + Friction { static_coefficient: mu, ..Friction::NONE } + } + + pub fn is_frictionless(&self) -> bool { + self.static_coefficient == 0.0 && self.kinetic.is_none_or(|k| k == 0.0) + } + + /// The kinetic coefficient (the static one when none is given). + pub fn kinetic_coefficient(&self) -> f64 { + self.kinetic.unwrap_or(self.static_coefficient) + } + + /// The coefficient at slip rate `rate` (≥ 0). + pub fn coefficient(&self, rate: f64) -> f64 { + let (s, k) = (self.static_coefficient, self.kinetic_coefficient()); + if self.decay == 0.0 { s } else { k + (s - k) * (-self.decay * rate.abs()).exp() } + } + + pub fn validate(&self) -> Result<()> { + for (what, v) in [("static friction coefficient", self.static_coefficient), ("kinetic friction coefficient", self.kinetic_coefficient())] { + if !(0.0..=MAX_FRICTION).contains(&v) { + return Err(Error::Invalid(format!("the {what} is in 0..={MAX_FRICTION}, got {v}"))); + } + } + if !(0.0..=MAX_DECAY).contains(&self.decay) { + return Err(Error::Invalid(format!("the friction decay coefficient is in 0..={MAX_DECAY:e}, got {}", self.decay))); + } + if self.decay == 0.0 && self.kinetic_coefficient() != self.static_coefficient { + return Err(Error::Invalid("a kinetic friction coefficient other than the static one needs a decay coefficient".into())); + } + Ok(()) + } +} + +/// Which faces a contact puts against which. +#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum ContactSurfaces { + /// The faces of two sets, each surface's nodes against the other's faces. + Pair { main: String, secondary: String }, + /// Every face against every other face of a set's faces (single-surface or self contact), or + /// of every exterior face of the model's solid elements when no set is named. + General { + #[serde(default, skip_serializing_if = "Option::is_none")] + surface: Option, + }, +} + +impl ContactSurfaces { + /// The sets it names. + pub fn sets(&self) -> Vec<&str> { + match self { + ContactSurfaces::Pair { main, secondary } => vec![main.as_str(), secondary.as_str()], + ContactSurfaces::General { surface } => surface.iter().map(String::as_str).collect(), + } + } +} + +/// What happens to nodes that start inside the surface they are in contact with. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum InitialPenetration { + /// Move them onto the surface before the analysis (strain-free), and report them. + #[default] + Adjust, + /// Refuse the analysis, naming the first one. + Refuse, +} + +/// A contact definition. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Contact { + pub surfaces: ContactSurfaces, + #[serde(default, skip_serializing_if = "Friction::is_frictionless")] + pub friction: Friction, + /// A factor on the solver's penalty stiffness. + #[serde(default = "one")] + pub penalty_scale: f64, + #[serde(default)] + pub initial_penetration: InitialPenetration, +} + +fn one() -> f64 { + 1.0 +} + +impl Contact { + /// A frictionless contact with the default penalty and initial adjustment. + pub fn new(surfaces: ContactSurfaces) -> Contact { + Contact { surfaces, friction: Friction::NONE, penalty_scale: 1.0, initial_penetration: InitialPenetration::Adjust } + } + + pub fn validate(&self) -> Result<()> { + if let ContactSurfaces::Pair { main, secondary } = &self.surfaces + && main == secondary + { + return Err(Error::Invalid(format!( + "a contact pair needs two different surfaces, got `{main}` twice (self contact is general contact on that surface)" + ))); + } + self.friction.validate()?; + let (lo, hi) = PENALTY_SCALE_RANGE; + if !(lo..=hi).contains(&self.penalty_scale) { + return Err(Error::Invalid(format!("the penalty scale is in {lo:e}..={hi}, got {}", self.penalty_scale))); + } + Ok(()) + } +} + +impl Model { + pub fn contacts(&self) -> &BTreeMap { + &self.contacts + } + + pub fn contact(&self, name: &str) -> Option<&Contact> { + self.contacts.get(name) + } + + fn check_contact(&self, name: &str, c: &Contact) -> Result<()> { + c.validate().map_err(|e| context("contact", name, e))?; + for set in c.surfaces.sets() { + require(&self.sets, "set", set).map_err(|e| context("contact", name, e))?; + } + Ok(()) + } + + pub fn add_contact(&mut self, name: &str, contact: Contact) -> Result<()> { + self.check_contact(name, &contact)?; + if self.contacts.len() >= MAX_CONTACTS { + return Err(Error::Limit(format!("a model holds at most {MAX_CONTACTS} contacts"))); + } + insert_new(&mut self.contacts, "contact", name, contact) + } + + pub fn replace_contact(&mut self, name: &str, contact: Contact) -> Result<()> { + self.check_contact(name, &contact)?; + replace_existing(&mut self.contacts, "contact", name, contact) + } + + pub fn remove_contact(&mut self, name: &str) -> Result<()> { + remove_existing(&mut self.contacts, "contact", name).map(|_| ()) + } + + /// The exterior faces of the solid elements (the faces no other solid element shares), in + /// element and face order: the surface of general contact without a named set. + pub fn exterior_faces(&self) -> Vec { + let mut keyed: Vec<([u32; 4], crate::ElementFace)> = Vec::new(); + for (label, e) in self.elements.iter().filter(|(_, e)| e.kind.dimension() == 3) { + for (i, local) in e.kind.faces().iter().enumerate() { + // A face is known by its sorted corner nodes (three or four). + let corners = if local.len() == 6 || local.len() == 3 { 3 } else { 4 }; + let mut key = [u32::MAX; 4]; + for (k, n) in key.iter_mut().zip(local.iter().take(corners)) { + *k = e.nodes.get(*n).copied().unwrap_or(u32::MAX); + } + key.sort_unstable(); + keyed.push((key, crate::ElementFace { element: *label, face: i as u8 + 1 })); + } + } + keyed.sort_unstable(); + let mut out: Vec = keyed + .chunk_by(|a, b| a.0 == b.0) + .filter_map(|run| match run { + [(_, face)] => Some(*face), + _ => None, + }) + .collect(); + out.sort_unstable(); + out + } + + /// The contacts that name a set, as ``contact `name` `` strings. + pub(crate) fn contact_users<'a>(&'a self, set: &'a str) -> impl Iterator + 'a { + self.contacts.iter().filter(move |(_, c)| c.surfaces.sets().contains(&set)).map(|(n, _)| format!("contact `{n}`")) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::Set; + + fn model() -> Model { + let mut m = Model::default(); + for s in ["top", "bottom"] { + m.add_set(s, Set { nodes: vec![1], ..Set::default() }).unwrap(); + } + m + } + + #[test] + fn contacts_validate_refer_to_sets_and_round_trip() { + let mut m = model(); + let json = serde_json::json!({ + "surfaces": {"type": "pair", "main": "top", "secondary": "bottom"}, + "friction": {"static": 0.3, "kinetic": 0.2, "decay": 5.0}, "penalty_scale": 0.5 + }); + let c: Contact = serde_json::from_value(json).unwrap(); + assert_eq!((c.friction.static_coefficient, c.friction.kinetic_coefficient(), c.initial_penetration), (0.3, 0.2, InitialPenetration::Adjust)); + assert_eq!(serde_json::from_value::(serde_json::to_value(&c).unwrap()).unwrap(), c); + m.add_contact("impact", c.clone()).unwrap(); + assert!(matches!(m.remove_set("top"), Err(Error::InUse(_))), "the contact refers to the set"); + assert!(m.add_contact("impact", c.clone()).is_err(), "names stay unique"); + let general: Contact = serde_json::from_value(serde_json::json!({"surfaces": {"type": "general"}})).unwrap(); + assert_eq!(general.surfaces, ContactSurfaces::General { surface: None }); + assert!(general.friction.is_frictionless() && general.penalty_scale == 1.0); + m.add_contact("everything", general).unwrap(); + m.remove_contact("impact").unwrap(); + m.remove_set("top").unwrap(); + assert!(m.remove_contact("impact").is_err()); + assert!(serde_json::from_value::(serde_json::json!({"surfaces": {"type": "general"}, "mu": 1})).is_err(), "typos are errors"); + + let base = Contact::new(ContactSurfaces::Pair { main: "bottom".into(), secondary: "top".into() }); + let bad = [ + Contact { friction: Friction::coulomb(-0.1), ..base.clone() }, + Contact { friction: Friction::coulomb(f64::NAN), ..base.clone() }, + Contact { friction: Friction { kinetic: Some(f64::INFINITY), ..Friction::coulomb(0.1) }, ..base.clone() }, + Contact { friction: Friction { decay: -1.0, ..Friction::coulomb(0.1) }, ..base.clone() }, + Contact { penalty_scale: 0.0, ..base.clone() }, + Contact { penalty_scale: f64::NAN, ..base.clone() }, + Contact { penalty_scale: 1.5, ..base.clone() }, + Contact::new(ContactSurfaces::Pair { main: "bottom".into(), secondary: "bottom".into() }), + ]; + for b in bad { + assert!(b.validate().is_err(), "{b:?}"); + } + let missing = Contact::new(ContactSurfaces::General { surface: Some("nope".into()) }); + assert!(matches!(m.add_contact("x", missing), Err(Error::NotFound(_)))); + } + + #[test] + fn friction_decays_from_static_to_kinetic() { + let f = Friction { static_coefficient: 0.5, kinetic: Some(0.2), decay: 2.0 }; + assert_eq!(f.coefficient(0.0), 0.5); + assert!((f.coefficient(1.0) - (0.2 + 0.3 * (-2.0f64).exp())).abs() < 1e-15); + assert!((f.coefficient(1e6) - 0.2).abs() < 1e-15); + assert_eq!(Friction::coulomb(0.3).coefficient(10.0), 0.3); + // A kinetic coefficient acts through the decay. + assert!(Friction { kinetic: Some(0.1), ..Friction::coulomb(0.3) }.validate().is_err()); + assert!(Friction { kinetic: Some(0.3), ..Friction::coulomb(0.3) }.validate().is_ok()); + assert!(Friction::NONE.is_frictionless() && !Friction::coulomb(0.1).is_frictionless()); + } + + #[test] + fn the_check_wants_solid_faces_on_contact_surfaces() { + let mut m = Model::default(); + let 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.]]; + for (i, x) in corners.iter().enumerate() { + m.set_node(i as u32 + 1, *x).unwrap(); + } + m.set_element(1, crate::ElementKind::Hex8, (1..=8).collect()).unwrap(); + m.add_set("nodes", Set { nodes: vec![1, 2], ..Set::default() }).unwrap(); + m.add_set("top", Set { faces: vec![crate::ElementFace { element: 1, face: 2 }], ..Set::default() }).unwrap(); + m.add_contact("c", Contact::new(ContactSurfaces::Pair { main: "top".into(), secondary: "nodes".into() })).unwrap(); + let errors: Vec = m.check().iter().filter(|d| d.severity == crate::Severity::Error).map(ToString::to_string).collect(); + assert!(errors.iter().any(|e| e.contains("contact `c`") && e.contains("set `nodes` has no element faces")), "{errors:?}"); + assert!(!errors.iter().any(|e| e.contains("set `top`")), "{errors:?}"); + } + + #[test] + fn exterior_faces_are_those_of_one_solid() { + // Two hexahedra sharing a face, and a tetrahedron on the top of the second. + let mut m = Model::default(); + for k in 0..3 { + for (i, xy) in [[0., 0.], [1., 0.], [1., 1.], [0., 1.]].iter().enumerate() { + m.set_node(4 * k + i as u32 + 1, [xy[0], xy[1], k as f64]).unwrap(); + } + } + m.set_node(13, [0.5, 0.5, 3.0]).unwrap(); + m.set_element(1, crate::ElementKind::Hex8, (1..=8).collect()).unwrap(); + m.set_element(2, crate::ElementKind::Hex8, (5..=12).collect()).unwrap(); + m.set_element(3, crate::ElementKind::Tet4, vec![9, 10, 11, 13]).unwrap(); + let faces = m.exterior_faces(); + let of = |e: u32| faces.iter().filter(|f| f.element == e).map(|f| f.face).collect::>(); + assert_eq!(of(1), [1, 3, 4, 5, 6], "all but the shared top"); + assert_eq!(of(2), [2, 3, 4, 5, 6], "all but the shared bottom (its top meets the tet over half its area)"); + assert_eq!(of(3), [1, 2, 3, 4], "a triangle never matches a quad"); + assert!(Model::default().exterior_faces().is_empty()); + } + + #[test] + fn a_model_holds_a_bounded_number_of_contacts() { + let mut m = model(); + for i in 0..MAX_CONTACTS { + m.add_contact(&format!("c{i}"), Contact::new(ContactSurfaces::General { surface: None })).unwrap(); + } + assert!(matches!(m.add_contact("one-more", Contact::new(ContactSurfaces::General { surface: None })), Err(Error::Limit(_)))); + } +} diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index 958b516..cbee8d5 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -34,6 +34,7 @@ pub mod bcs; pub mod check; pub mod constraints; +pub mod contact; pub mod discrete; pub mod fluids; pub mod geometry; @@ -53,6 +54,7 @@ use serde::{Deserialize, Serialize}; pub use bcs::{BoundaryCondition, Dofs, Symmetry, TEMPERATURE_DOF, Target}; pub use check::{Diagnostic, Resolved, ResolvedSet, Severity}; pub use constraints::{Coupling, CouplingKind, ReferencePoint}; +pub use contact::{Contact, ContactSurfaces, Friction, InitialPenetration}; pub use discrete::{Direction, Mass, PointMass, ResolvedDirection, ResolvedMass, ResolvedSpring, Spring}; pub use fluids::{BoxFace, BoxFaces, FaceCondition, FluidBoundary, FluidDomain, Profile, SurfaceBoundary, SurfaceCondition}; pub use geometry::{Classification, Entity, GeometryRef}; @@ -123,6 +125,8 @@ pub struct Model { reference_points: BTreeMap, couplings: BTreeMap, #[serde(skip_serializing_if = "BTreeMap::is_empty")] + contacts: BTreeMap, + #[serde(skip_serializing_if = "BTreeMap::is_empty")] masses: BTreeMap, #[serde(skip_serializing_if = "BTreeMap::is_empty")] springs: BTreeMap, @@ -157,6 +161,8 @@ struct ModelData { #[serde(default)] couplings: BTreeMap, #[serde(default)] + contacts: BTreeMap, + #[serde(default)] masses: BTreeMap, #[serde(default)] springs: BTreeMap, @@ -215,6 +221,9 @@ impl TryFrom for Model { for (n, v) in d.couplings { m.add_coupling(&n, v)?; } + for (n, v) in d.contacts { + m.add_contact(&n, v)?; + } for (n, v) in d.masses { m.add_mass(&n, v)?; } @@ -330,6 +339,7 @@ impl Model { sections: self.sections.clone(), reference_points: self.reference_points.clone(), couplings: self.couplings.clone(), + contacts: self.contacts.clone(), masses: self.masses.clone(), springs: self.springs.clone(), boundary_conditions: self.boundary_conditions.clone(), @@ -340,13 +350,14 @@ impl Model { } /// Counts of each kind of definition, for summaries. - pub fn definition_counts(&self) -> [(&'static str, usize); 11] { + pub fn definition_counts(&self) -> [(&'static str, usize); 12] { [ ("sets", self.sets.len()), ("materials", self.materials.len()), ("sections", self.sections.len()), ("reference_points", self.reference_points.len()), ("couplings", self.couplings.len()), + ("contacts", self.contacts.len()), ("masses", self.masses.len()), ("springs", self.springs.len()), ("boundary_conditions", self.boundary_conditions.len()), diff --git a/crates/model/src/sets.rs b/crates/model/src/sets.rs index 87e2cfb..84964c7 100644 --- a/crates/model/src/sets.rs +++ b/crates/model/src/sets.rs @@ -122,6 +122,7 @@ impl Model { let mut users = Vec::new(); users.extend(self.sections.iter().filter(|(_, s)| s.set == set).map(|(n, _)| format!("section `{n}`"))); users.extend(self.couplings.iter().filter(|(_, c)| c.set == set).map(|(n, _)| format!("coupling `{n}`"))); + users.extend(self.contact_users(set)); users.extend(self.discrete_users(|t| t.set() == Some(set))); users .extend(self.boundary_conditions.iter().filter(|(_, b)| b.target().set() == Some(set)).map(|(n, _)| format!("boundary condition `{n}`"))); diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index e0a6f90..67e808e 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -10,7 +10,7 @@ use std::collections::BTreeSet; use serde::{Deserialize, Serialize}; -use crate::{Error, Model, Result, check_name, context, finite, positive, require}; +use crate::{Error, Friction, Model, Result, check_name, context, finite, positive, require}; /// Most steps in a model. pub const MAX_STEPS: usize = 1000; @@ -317,7 +317,7 @@ pub struct ExplicitDynamic { /// the last one listed). #[serde(default, skip_serializing_if = "Vec::is_empty")] pub initial_velocities: Vec, - /// Fixed, frictionless rigid planes the nodes of sets cannot pass through. + /// Fixed rigid planes the nodes of sets cannot pass through, with Coulomb friction. #[serde(default, skip_serializing_if = "Vec::is_empty")] pub rigid_walls: Vec, /// Hourglass control of one-point hexahedra. @@ -431,6 +431,7 @@ impl ExplicitDynamic { if !(w.normal.iter().map(|c| c * c).sum::() > 0.0) { return Err(Error::Invalid("a rigid wall's normal must not be zero".into())); } + w.friction.validate().map_err(|e| Error::Invalid(format!("rigid wall on set `{}`: {e}", w.set)))?; } self.hourglass.validate()?; self.bulk_viscosity.validate() @@ -445,14 +446,16 @@ pub struct InitialVelocity { pub velocity: [f64; 3], } -/// A fixed, frictionless rigid plane through `point`. Its `normal` points to the side the nodes -/// of `set` stay on. +/// A fixed rigid plane through `point`. Its `normal` points to the side the nodes of `set` stay +/// on; nodes sliding on it feel `friction` (frictionless by default). #[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] #[serde(deny_unknown_fields)] pub struct RigidWall { pub set: String, pub point: [f64; 3], pub normal: [f64; 3], + #[serde(default, skip_serializing_if = "Friction::is_frictionless")] + pub friction: Friction, } /// How one-point hexahedra resist their hourglass (zero-energy) modes. @@ -870,7 +873,7 @@ mod tests { let json = serde_json::json!({ "type": "explicit_dynamic", "end_time": 1e-3, "history_nodes": "front", "initial_velocities": [{"set": "front", "velocity": [-10, 0, 0]}], - "rigid_walls": [{"set": "front", "point": [0, 0, 0], "normal": [1, 0, 0]}], + "rigid_walls": [{"set": "front", "point": [0, 0, 0], "normal": [1, 0, 0], "friction": {"static": 0.2}}], "hourglass": {"type": "stiffness"} }); let p: Procedure = serde_json::from_value(json).unwrap(); @@ -895,7 +898,14 @@ mod tests { ExplicitDynamic { max_increments: 0, ..base.clone() }, ExplicitDynamic { output_interval: Some(1e-9), ..base.clone() }, ExplicitDynamic { history_interval: Some(-1.0), ..base.clone() }, - ExplicitDynamic { rigid_walls: vec![RigidWall { set: "front".into(), point: [0.0; 3], normal: [0.0; 3] }], ..base.clone() }, + ExplicitDynamic { + rigid_walls: vec![RigidWall { set: "front".into(), point: [0.0; 3], normal: [0.0; 3], friction: Friction::NONE }], + ..base.clone() + }, + ExplicitDynamic { + rigid_walls: vec![RigidWall { set: "front".into(), point: [0.0; 3], normal: [1.0, 0.0, 0.0], friction: Friction::coulomb(f64::NAN) }], + ..base.clone() + }, ExplicitDynamic { initial_velocities: vec![InitialVelocity { set: "front".into(), velocity: [f64::INFINITY, 0.0, 0.0] }], ..base.clone() diff --git a/crates/render/src/glyphs.rs b/crates/render/src/glyphs.rs index 85cbebb..1992ddf 100644 --- a/crates/render/src/glyphs.rs +++ b/crates/render/src/glyphs.rs @@ -874,12 +874,12 @@ mod tests { #[test] fn rigid_walls_back_their_set_and_initial_velocities_leave_it() { - use fieldcraft_model::{ExplicitDynamic, InitialVelocity, Procedure, RigidWall, Step}; + use fieldcraft_model::{ExplicitDynamic, Friction, InitialVelocity, Procedure, RigidWall, Step}; let mut m = bar(4); m.add_set("all", Set { nodes: m.nodes().keys().copied().collect(), ..Set::default() }).unwrap(); let e = ExplicitDynamic { initial_velocities: vec![InitialVelocity { set: "all".into(), velocity: [-10.0, 0.0, 0.0] }], - rigid_walls: vec![RigidWall { set: "root".into(), point: [0.0, 0.0, 0.0], normal: [1.0, 0.0, 0.0] }], + rigid_walls: vec![RigidWall { set: "root".into(), point: [0.0, 0.0, 0.0], normal: [1.0, 0.0, 0.0], friction: Friction::NONE }], ..ExplicitDynamic::new(1.0) }; m.add_step(Step::new("Impact", Procedure::ExplicitDynamic(e)), None).unwrap(); diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 7657768..9ca9391 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -418,6 +418,11 @@ 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(); diff --git a/crates/solver/tests/couplings.rs b/crates/solver/tests/couplings.rs index 0309f28..0f3b394 100644 --- a/crates/solver/tests/couplings.rs +++ b/crates/solver/tests/couplings.rs @@ -390,3 +390,14 @@ fn oversized_reductions_are_refused_before_allocating() { assert!(e.contains("reduced matrices need more than"), "{e}"); assert!(solve_static(&g.model, Some("Step"), &Limits::default()).is_ok()); } + +/// Contact belongs to explicit dynamics so far: the implicit solves refuse it by name. +#[test] +fn implicit_solves_refuse_contact_by_name() { + let mut g = grid([1, 1, 1], ElementKind::Hex8, |x| x); + g.define(iso(1000.0, 0.3), Some(1.0), full(), Procedure::Static {}); + g.face_set("top", |x| near(x[2], 1.0)); + g.model.add_contact("touch", fieldcraft_model::Contact::new(fieldcraft_model::ContactSurfaces::General { surface: Some("top".into()) })).unwrap(); + let e = solve_static(&g.model, Some("Step"), &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("contact (`touch`)") && e.contains("not supported"), "{e}"); +} diff --git a/crates/ui-egui/src/browser.rs b/crates/ui-egui/src/browser.rs index 3aa88f0..b680339 100644 --- a/crates/ui-egui/src/browser.rs +++ b/crates/ui-egui/src/browser.rs @@ -33,6 +33,7 @@ pub enum Def { Load, ReferencePoint, Coupling, + Contact, FluidDomain, } @@ -48,6 +49,7 @@ impl Def { Def::Load => "load", Def::ReferencePoint => "reference_point", Def::Coupling => "coupling", + Def::Contact => "contact", Def::FluidDomain => "fluid_domain", } } @@ -62,6 +64,7 @@ impl Def { Def::Load => "Load", Def::ReferencePoint => "Reference Point", Def::Coupling => "Coupling", + Def::Contact => "Contact", Def::FluidDomain => "Fluid Domain", } } @@ -77,13 +80,19 @@ impl Def { Def::Load => "model.load.create", Def::ReferencePoint => "model.reference_point.create", Def::Coupling => "model.coupling.create", + Def::Contact => "interaction.contact.pair", Def::FluidDomain => "model.fluid_domain.create", } } - /// The command of this kind with `verb` (`edit`, `delete`). + /// The command of this kind with `verb` (`create`, `edit`, `delete`); contacts are + /// interactions, created as pairs or general contact. pub fn command(self, verb: &str) -> String { - format!("model.{}.{verb}", self.stem()) + match (self, verb) { + (Def::Contact, "create") => "interaction.contact.pair".into(), + (Def::Contact, _) => format!("interaction.contact.{verb}"), + _ => format!("model.{}.{verb}", self.stem()), + } } /// How many definitions of this kind the model has. @@ -97,6 +106,7 @@ impl Def { Def::Load => m.loads().len(), Def::ReferencePoint => m.reference_points().len(), Def::Coupling => m.couplings().len(), + Def::Contact => m.contacts().len(), Def::FluidDomain => m.fluid_domains().len(), } } @@ -112,6 +122,7 @@ impl Def { Def::Load => m.loads().contains_key(name), Def::ReferencePoint => m.reference_points().contains_key(name), Def::Coupling => m.couplings().contains_key(name), + Def::Contact => m.contacts().contains_key(name), Def::FluidDomain => m.fluid_domain(name).is_some(), } } @@ -127,6 +138,7 @@ impl Def { Def::Load => m.loads().keys().cloned().collect(), Def::ReferencePoint => m.reference_points().keys().cloned().collect(), Def::Coupling => m.couplings().keys().cloned().collect(), + Def::Contact => m.contacts().keys().cloned().collect(), Def::FluidDomain => m.fluid_domains().keys().cloned().collect(), } } @@ -143,6 +155,7 @@ pub fn definition(m: &Model, kind: Def, name: &str) -> Option { Def::Load => serde_json::to_value(m.loads().get(name)?), Def::ReferencePoint => serde_json::to_value(m.reference_points().get(name)?), Def::Coupling => serde_json::to_value(m.couplings().get(name)?), + Def::Contact => serde_json::to_value(m.contacts().get(name)?), Def::FluidDomain => serde_json::to_value(m.fluid_domain(name)?), } .ok()?; @@ -442,7 +455,7 @@ fn rows(app: &mut ModelApp, ui: &mut egui::Ui, acts: &mut Vec) { ("sections", "Sections", &[Def::Section]), ("steps", "Steps", &[Def::Step]), ("loads", "Loads & BCs", &[Def::BoundaryCondition, Def::Load]), - ("interactions", "Interactions", &[Def::ReferencePoint, Def::Coupling]), + ("interactions", "Interactions", &[Def::ReferencePoint, Def::Coupling, Def::Contact]), ("fluids", "Fluid Domains", &[Def::FluidDomain]), ]; for (key, label, kinds) in groups { @@ -550,6 +563,7 @@ mod tests { Def::Load, Def::ReferencePoint, Def::Coupling, + Def::Contact, Def::FluidDomain, ] { for verb in ["create", "edit", "delete"] { diff --git a/docs/features.md b/docs/features.md index 8470c1a..071152d 100644 --- a/docs/features.md +++ b/docs/features.md @@ -6,7 +6,7 @@ Every entry is a toolbar button. **Ready** entries run their command; **coming s shown greyed with the roadmap milestone that brings them (ROADMAP.md). Scripts, the command line, the control channel and AI agents over MCP (docs/mcp.md) run the same commands. -56 ready, 48 coming soon. +59 ready, 45 coming soon. ## Geometry @@ -103,9 +103,9 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands. | Constraints | Reference Point — A reference point for couplings, loads and boundary conditions | ready | `model.reference_point.create`, `model.reference_point.edit`, `model.reference_point.delete` | | 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 touches, with friction | coming soon (M7) | | -| Contact | Contact Pair — Contact between two surfaces | coming soon (M7) | | -| Contact | Find Contacts — Detect contact pairs automatically | 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 | 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) | | | Connectors | Spring — Spring or dashpot to the ground or between two points, along a direction or their line | ready | `model.spring.create`, `model.spring.edit`, `model.spring.delete` |