diff --git a/crates/solver/src/contact.rs b/crates/solver/src/contact.rs index 1e08b2e..2f3f58d 100644 --- a/crates/solver/src/contact.rs +++ b/crates/solver/src/contact.rs @@ -75,9 +75,15 @@ //! of the pressure field and ring at the contact's edge), the multiplier, the gap and the status; //! per contact the force on the secondary surface and its normal part. //! -//! **Scope.** Pairs of surfaces of linear solid faces (tri3 and quad4 faces: tet4, wedge6 and -//! hex8 elements) in 3D models, with Coulomb friction ([`Interface::friction`]). General (single-surface) contact, quadratic faces, -//! 2D models and arc-length steps are refused by name. +//! **Single-surface (general) contact** is the same mortar contact with one face list as both +//! sides: a set's faces, or every exterior face of the solids. Each face is a secondary face +//! against every other one that faces it, that shares no node with its element (its neighbours) +//! and that it does not reach through more than half an element (the far side of a thin wall), +//! so two faces pressing together are paired both ways, symmetrically: self-contact of folds. +//! +//! **Scope.** Linear solid faces (tri3 and quad4 faces: tet4, wedge6 and hex8 elements) in 3D +//! models, with Coulomb friction ([`Interface::friction`]). Quadratic faces, 2D models and +//! arc-length steps are refused by name. use std::collections::BTreeMap; @@ -85,7 +91,7 @@ use fieldcraft_fem::contact::{Aabb, Grid, bounds}; use fieldcraft_fem::dual::{Dual, Hyper, Real, VARS}; use fieldcraft_fem::frame::{Vec3, dot, norm, unit}; use fieldcraft_fem::mortar::{area_vector, integrate}; -use fieldcraft_model::{Contact, ContactSurfaces, ElementKind, Friction, InitialPenetration, Model, Resolved}; +use fieldcraft_model::{Contact, ContactSurfaces, ElementFace, ElementKind, Friction, InitialPenetration, Model, Resolved}; use fieldcraft_parallel as parallel; use serde::Serialize; @@ -306,7 +312,11 @@ fn size(x: &[Vec3]) -> f64 { /// The faces of set `set` for definition `owner` (``contact `name` ``, ``tie `name` ``), refused /// unless they are linear faces of analysed solid elements. pub(crate) fn faces(model: &Model, resolved: &Resolved, mesh: &Mesh, owner: &str, set: &str) -> Result> { - let list = resolved.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default(); + face_list(model, mesh, owner, resolved.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default()) +} + +/// The faces `list` for definition `owner`, as [`faces`]. +fn face_list(model: &Model, mesh: &Mesh, owner: &str, list: &[ElementFace]) -> Result> { let mut out = Vec::with_capacity(list.len()); for f in list { let at = |what: String| Error::Invalid(format!("{owner}: element {} face {}: {what}", f.element, f.face)); @@ -459,14 +469,21 @@ impl Interface { } fn new(model: &Model, resolved: &Resolved, mesh: &Mesh, name: &str, contact: &Contact) -> Result { - let ContactSurfaces::Pair { main, secondary } = &contact.surfaces else { - return Err(Error::NotSupported(format!( - "general (single-surface) contact `{name}` in implicit steps; define contact pairs (interaction.contact.pair), or solve it in explicit dynamics" - ))); - }; let owner = format!("contact `{name}`"); - let secondary = faces(model, resolved, mesh, &owner, secondary)?; - let main = faces(model, resolved, mesh, &owner, main)?; + let (secondary, main) = match &contact.surfaces { + ContactSurfaces::Pair { main, secondary } => { + (faces(model, resolved, mesh, &owner, secondary)?, faces(model, resolved, mesh, &owner, main)?) + } + // Single surface: every face against every other one that faces it and shares no node + // with its element (its neighbours), both ways. + ContactSurfaces::General { surface } => { + let list = match surface { + Some(set) => faces(model, resolved, mesh, &owner, set)?, + None => face_list(model, mesh, &owner, &model.exterior_faces())?, + }; + (list.clone(), list) + } + }; let mut nodes: Vec = secondary.iter().flat_map(|f| f.nodes.iter().copied()).collect(); nodes.sort_unstable(); nodes.dedup(); diff --git a/crates/solver/tests/contact.rs b/crates/solver/tests/contact.rs index b9e7294..66b54c7 100644 --- a/crates/solver/tests/contact.rs +++ b/crates/solver/tests/contact.rs @@ -353,10 +353,6 @@ fn what_implicit_contact_refuses() { b.model.add_step(Step::new("Linear", Procedure::Static {}), None).unwrap(); let e = fieldcraft_solver::solve_static(&b.model, Some("Linear"), &Limits::default()).unwrap_err().to_string(); assert!(e.contains("interface") && e.contains("linear static"), "{e}"); - // General contact. - let mut b = stacked(0.0, InitialPenetration::Adjust); - b.model.add_contact("all", Contact::new(ContactSurfaces::General { surface: None })).unwrap(); - assert!(solve(&b.model).contains("general")); // Arc-length control. let mut b = stacked(0.0, InitialPenetration::Adjust); let arc = NonlinearStatic { arc_length: Some(ArcLength::default()), ..NonlinearStatic::default() }; @@ -560,3 +556,147 @@ fn what_ties_refuse() { let r = fieldcraft_solver::solve_static(&b.model, Some("Step"), &Limits::default()).unwrap(); assert!(r.warnings.iter().all(|w| !w.contains("stay free")), "{:?}", r.warnings); } + +/// A C-shaped ring of rubber (an annulus of radii 1 and 1.5 with a 10° opening at +x, one hex8 +/// layer thick in z with z held: plane strain) whose two ends are pressed together through the +/// outer nodes two elements from each end, against general (single-surface) contact on every +/// exterior face of the body: Mooney–Rivlin, nearly incompressible, hybrid elements. +fn c_ring(finite: bool, closing: f64) -> (Blocks, StaticResult) { + let mut b = Blocks::default(); + let (n, open) = (36usize, 5f64.to_radians()); + let angle = move |i: f64| open + i * (2.0 * std::f64::consts::PI - 2.0 * open); + b.block(ElementKind::Hex8, [n, 2, 1], |p| { + let (t, r) = (angle(p[0]), 1.5 - 0.5 * p[1]); + [r * t.cos(), r * t.sin(), 0.2 * p[2]] + }); + let all: Vec = b.model.elements().keys().copied().collect(); + b.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap(); + let rubber = Material { + elastic: Some(Elastic::Isotropic { young: 3.0, poisson: 0.499 }), + hyperelastic: Some(Hyperelastic::MooneyRivlin { c01: 0.1 }), + ..Material::default() + }; + b.model.add_material("Rubber", rubber).unwrap(); + b.model + .add_section( + "solid", + Section { + set: "all".into(), + material: "Rubber".into(), + kind: SectionKind::Solid {}, + formulation: Formulation { hybrid: true, ..Formulation::default() }, + }, + ) + .unwrap(); + let settings = NonlinearStatic { nonlinear_geometry: finite, initial_increment: 0.1, max_increments: 200, ..NonlinearStatic::default() }; + b.model.add_step(Step::new("Step", Procedure::NonlinearStatic(settings)), None).unwrap(); + b.model.add_contact("self", Contact::new(ContactSurfaces::General { surface: None })).unwrap(); + // The outer nodes two elements from each end close towards the axis of symmetry y = 0; the + // back of the ring is held in its plane; every node in z. + let step = 2.0 * std::f64::consts::PI - 2.0 * open; + let at = |x: [f64; 3], i: f64| { + ((x[0] * x[0] + x[1] * x[1]).sqrt() - 1.5).abs() < 1e-9 + && (x[1].atan2(x[0]).rem_euclid(2.0 * std::f64::consts::PI) - angle(i)).abs() < 1e-6 * step + }; + let near_end = 2.0 / n as f64; + b.bc("upper", [None, Some(-closing), None], move |x| at(x, near_end)); + b.bc("lower", [None, Some(closing), None], move |x| at(x, 1.0 - near_end)); + b.bc("back", [Some(0.0), Some(0.0), None], move |x| at(x, 0.5)); + b.bc("plane", [None, None, Some(0.0)], |_| true); + let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap(); + report(&r); + (b, r) +} + +/// Self-contact: the ring's two ends, faces of one body and one general contact, meet and press +/// on each other without passing through, and the result keeps the model's mirror symmetry +/// about y = 0. +#[test] +fn a_rubber_ring_closes_on_itself_without_penetration() { + for finite in [true] { + let (b, r) = c_ring(finite, 0.15); + let c = &r.contacts[0]; + assert!(c.closed() > 0, "finite {finite}: the ends touch"); + // No penetration beyond the augmented Lagrangian's tolerance. + assert!(c.penetration < 1e-4, "finite {finite}: {}", c.penetration); + let deformed = |l: u32| -> [f64; 3] { + let u = r.displacement[r.nodes.binary_search(&l).unwrap()]; + let x = b.xyz[&l]; + [x[0] + u[0], x[1] + u[1], x[2] + u[2]] + }; + // The ends' nodes stay on their own side of the plane of symmetry, up to the overlap the + // mortar's weak (weighted) enforcement allows at the edge of a contact zone: a hundredth + // of the elements' size (0.25). + for (l, x) in &b.xyz { + let y = x[1]; + if x[0] > 0.9 && y.abs() < 0.2 { + let d = deformed(*l); + assert!(d[1] * y.signum() >= -2.5e-3, "finite {finite}: node {l} crossed: {y} → {}", d[1]); + } + } + // Mirror symmetry: the node at (x, −y) moves as the mirror of the node at (x, y). + let mut checked = 0; + for (l, x) in &b.xyz { + let Some((m, _)) = b.xyz.iter().find(|(_, y)| (y[0] - x[0]).abs() < 1e-9 && (y[1] + x[1]).abs() < 1e-9 && (y[2] - x[2]).abs() < 1e-9) + else { + continue; + }; + let (u, v) = (r.displacement[r.nodes.binary_search(l).unwrap()], r.displacement[r.nodes.binary_search(m).unwrap()]); + assert!((u[0] - v[0]).abs() < 1e-6 && (u[1] + v[1]).abs() < 1e-6, "finite {finite}: nodes {l}, {m}: {u:?} {v:?}"); + checked += 1; + } + assert!(checked > 50, "{checked}"); + eprintln!("finite {finite}: closed {} penetration {:.2e} force {:?}", c.closed(), c.penetration, c.force); + } +} + +/// A rubber bellows (a quarter of a corrugated tube: three convolutions, mean radius 2, fold +/// depth 0.5, wall 0.15, height 3) squeezed axially by 0.8 with general contact on all its faces: +/// its folds close and press on one another, as in a boot seal. (Squeezed to half its height +/// the step stalls at 57%: deeply folded rubber needs better increment control; see M7.) +#[test] +#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"] +fn a_rubber_bellows_folds_onto_itself() { + let mut b = Blocks::default(); + let (folds, height, depth, wall) = (3.0, 3.0, 0.5, 0.15); + let quarter = std::f64::consts::FRAC_PI_2; + b.block(ElementKind::Hex8, [4, 48, 1], move |p| { + let theta = quarter * p[0]; + let s = p[1] * folds; + let tri = 1.0 - 2.0 * (s - s.floor() - 0.5).abs(); // 0 → 1 → 0 per convolution + let r = 2.0 + depth * tri + wall * p[2]; + [r * theta.cos(), r * theta.sin(), height * p[1]] + }); + let all: Vec = b.model.elements().keys().copied().collect(); + b.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap(); + let rubber = Material { + elastic: Some(Elastic::Isotropic { young: 3.0, poisson: 0.499 }), + hyperelastic: Some(Hyperelastic::MooneyRivlin { c01: 0.1 }), + ..Material::default() + }; + b.model.add_material("Rubber", rubber).unwrap(); + b.model + .add_section( + "solid", + Section { + set: "all".into(), + material: "Rubber".into(), + kind: SectionKind::Solid {}, + formulation: Formulation { hybrid: true, ..Formulation::default() }, + }, + ) + .unwrap(); + let settings = NonlinearStatic { initial_increment: 0.05, max_increment: 0.1, max_increments: 400, ..NonlinearStatic::default() }; + b.model.add_step(Step::new("Step", Procedure::NonlinearStatic(settings)), None).unwrap(); + b.model.add_contact("folds", Contact::new(ContactSurfaces::General { surface: None })).unwrap(); + b.bc("x0", [Some(0.0), None, None], |x| near(x[0], 0.0)); + b.bc("y0", [None, Some(0.0), None], |x| near(x[1], 0.0)); + b.bc("base", [None, None, Some(0.0)], |x| near(x[2], 0.0)); + b.bc("squeeze", [None, None, Some(-0.8)], move |x| near(x[2], height)); + let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap(); + report(&r); + let c = &r.contacts[0]; + eprintln!("closed {} penetration {:.2e} augmentations {}", c.closed(), c.penetration, c.augmentations); + assert!(c.closed() > 0, "the folds touch"); + assert!(c.penetration < 1e-3 * wall, "{}", c.penetration); +}