forked from 1b-agent/fieldcraft-upstream
jyvikoandClaude Opus 5.5
85f0905db1
General contact, on a set's faces or every exterior face of the solids, is mortar contact with one face list as both sides: each face against every other one that faces it, shares no node with its element and is not reached through more than half an element, paired both ways. Tested by a nearly incompressible Mooney-Rivlin C-ring (hybrid hex8) closing on itself (no penetration beyond the mortar's weak enforcement, mirror symmetry to 1e-6) and, heavy, a rubber bellows whose folds press on each other. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo
703 lines
37 KiB
Rust
703 lines
37 KiB
Rust
//! Contact in nonlinear static steps against closed forms: the contact patch test on
|
||
//! non-matching meshes, Newton's convergence with contact, and what implicit contact refuses.
|
||
//! Run with `--nocapture` to see the iterations.
|
||
|
||
use std::collections::BTreeMap;
|
||
|
||
use fieldcraft_model::*;
|
||
use fieldcraft_solver::contact::Status;
|
||
use fieldcraft_solver::{Limits, StaticResult, solve_nonlinear};
|
||
|
||
/// A model built block by block.
|
||
#[derive(Default)]
|
||
struct Blocks {
|
||
model: Model,
|
||
xyz: BTreeMap<u32, [f64; 3]>,
|
||
next_node: u32,
|
||
next_element: u32,
|
||
}
|
||
|
||
impl Blocks {
|
||
/// A structured block of `n` cells over the unit cube mapped by `map`, hex8 or tet4 (six
|
||
/// Kuhn tetrahedra per cell); its elements' labels.
|
||
fn block(&mut self, kind: ElementKind, n: [usize; 3], map: impl Fn([f64; 3]) -> [f64; 3]) -> Vec<u32> {
|
||
let first = self.next_node + 1;
|
||
let label = |i: [usize; 3]| first + (i[0] + (n[0] + 1) * (i[1] + (n[1] + 1) * i[2])) as u32;
|
||
for k in 0..=n[2] {
|
||
for j in 0..=n[1] {
|
||
for i in 0..=n[0] {
|
||
let x = map([i as f64 / n[0] as f64, j as f64 / n[1] as f64, k as f64 / n[2] as f64]);
|
||
let l = label([i, j, k]);
|
||
self.model.set_node(l, x).unwrap();
|
||
self.xyz.insert(l, x);
|
||
self.next_node = self.next_node.max(l);
|
||
}
|
||
}
|
||
}
|
||
let mut out = Vec::new();
|
||
for k in 0..n[2] {
|
||
for j in 0..n[1] {
|
||
for i in 0..n[0] {
|
||
let at = |d: [usize; 3]| label([i + d[0], j + d[1], k + d[2]]);
|
||
let cells: Vec<Vec<u32>> = match kind {
|
||
ElementKind::Hex8 => vec![
|
||
[[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]].iter().map(|d| at(*d)).collect(),
|
||
],
|
||
_ => [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]]
|
||
.iter()
|
||
.map(|perm| {
|
||
let mut v = [[0usize; 3]; 4];
|
||
for s in 0..3 {
|
||
v[s + 1] = v[s];
|
||
v[s + 1][perm[s]] = 1;
|
||
}
|
||
let p = v.map(|c| c.map(|x| x as f64));
|
||
let e = |i: usize| [0, 1, 2].map(|k| p[i][k] - p[0][k]);
|
||
let (a, b, c) = (e(1), e(2), e(3));
|
||
let vol =
|
||
a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]);
|
||
if vol < 0.0 {
|
||
v.swap(1, 2);
|
||
}
|
||
v.iter().map(|d| at(*d)).collect()
|
||
})
|
||
.collect(),
|
||
};
|
||
for nodes in cells {
|
||
self.next_element += 1;
|
||
self.model.set_element(self.next_element, kind, nodes).unwrap();
|
||
out.push(self.next_element);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
out
|
||
}
|
||
|
||
/// The faces of `elements` whose nodes all satisfy `f`.
|
||
fn faces(&self, elements: &[u32], f: impl Fn([f64; 3]) -> bool) -> Vec<ElementFace> {
|
||
let mut out = Vec::new();
|
||
for l in elements {
|
||
let e = self.model.element(*l).unwrap();
|
||
for (i, face) in e.kind.faces().iter().enumerate() {
|
||
if face.iter().all(|k| f(self.xyz[&e.nodes[*k]])) {
|
||
out.push(ElementFace { element: *l, face: i as u8 + 1 });
|
||
}
|
||
}
|
||
}
|
||
out
|
||
}
|
||
|
||
fn nodes(&self, f: impl Fn([f64; 3]) -> bool) -> Vec<u32> {
|
||
self.xyz.iter().filter(|(_, x)| f(**x)).map(|(l, _)| *l).collect()
|
||
}
|
||
|
||
fn node_set(&mut self, name: &str, f: impl Fn([f64; 3]) -> bool) {
|
||
let nodes = self.nodes(f);
|
||
self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap();
|
||
}
|
||
|
||
fn face_set(&mut self, name: &str, elements: &[u32], f: impl Fn([f64; 3]) -> bool) {
|
||
let faces = self.faces(elements, f);
|
||
assert!(!faces.is_empty(), "no faces for `{name}`");
|
||
self.model.add_set(name, Set { faces, ..Set::default() }).unwrap();
|
||
}
|
||
|
||
fn bc(&mut self, name: &str, values: [Option<f64>; 3], f: impl Fn([f64; 3]) -> bool) {
|
||
self.node_set(name, f);
|
||
let [a, b, c] = values;
|
||
self.model
|
||
.add_boundary_condition(name, BoundaryCondition::Displacement { target: Target::Set(name.into()), values: [a, b, c, None, None, None] })
|
||
.unwrap();
|
||
self.model.set_boundary_condition_active("Step", name, true).unwrap();
|
||
}
|
||
|
||
/// One elastic material and section over every element, and the nonlinear step.
|
||
fn define(&mut self, young: f64, poisson: f64, settings: NonlinearStatic) {
|
||
let all: Vec<u32> = self.model.elements().keys().copied().collect();
|
||
self.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
|
||
self.model.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young, poisson }), ..Material::default() }).unwrap();
|
||
self.model
|
||
.add_section(
|
||
"solid",
|
||
Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() },
|
||
)
|
||
.unwrap();
|
||
self.model.add_step(Step::new("Step", Procedure::NonlinearStatic(settings)), None).unwrap();
|
||
}
|
||
|
||
fn pair(&mut self, name: &str, main: &str, secondary: &str) {
|
||
let c = Contact::new(ContactSurfaces::Pair { main: main.into(), secondary: secondary.into() });
|
||
self.model.add_contact(name, c).unwrap();
|
||
}
|
||
}
|
||
|
||
fn near(a: f64, b: f64) -> bool {
|
||
(a - b).abs() < 1e-9
|
||
}
|
||
|
||
fn report(r: &StaticResult) {
|
||
let h = r.nonlinear.as_ref().unwrap();
|
||
for inc in &h.increments {
|
||
let forces: Vec<String> = inc.iterations.iter().map(|i| format!("{:.1e}", i.force)).collect();
|
||
eprintln!("increment {} λ {:.3} aug {} [{}]", inc.number, inc.load_factor, inc.augmentations, forces.join(" "));
|
||
}
|
||
}
|
||
|
||
/// Two blocks pressed together across non-matching meshes (`lower` × `upper` cells, each of
|
||
/// `lower_kind` / `upper_kind`), between frictionless symmetry planes: a uniaxial stress.
|
||
fn patch(lower_kind: ElementKind, upper_kind: ElementKind, geometry: bool, distorted: bool) -> (Blocks, StaticResult) {
|
||
let mut b = Blocks::default();
|
||
// A graded lower block, its cells boxes of different sizes, or distorted (its faces stay the
|
||
// planes of the cube either way).
|
||
let low = if distorted {
|
||
b.block(lower_kind, [3, 3, 2], |p| [p[0] + 0.3 * p[0] * (1.0 - p[0]) * (p[1] - 0.5), p[1] + 0.3 * p[1] * (1.0 - p[1]) * (p[0] - 0.4), p[2]])
|
||
} else {
|
||
b.block(lower_kind, [3, 3, 2], |p| [p[0].powf(1.4), p[1].powf(0.8), p[2]])
|
||
};
|
||
let upper = [2, 2, 2];
|
||
let up = b.block(upper_kind, upper, |p| [p[0], p[1], 1.0 + p[2]]);
|
||
b.define(
|
||
1000.0,
|
||
0.3,
|
||
NonlinearStatic {
|
||
nonlinear_geometry: geometry,
|
||
initial_increment: 0.5,
|
||
// Tight for small sliding (linear: exact); finite sliding converges to the default ones.
|
||
tolerances: if geometry { Tolerances::default() } else { Tolerances { force: 1e-10, displacement: 1e-8, energy: 1e-16 } },
|
||
..NonlinearStatic::default()
|
||
},
|
||
);
|
||
b.face_set("lower_top", &low, |x| near(x[2], 1.0));
|
||
b.face_set("upper_bottom", &up, |x| near(x[2], 1.0));
|
||
b.pair("interface", "lower_top", "upper_bottom");
|
||
b.bc("x0", [Some(0.0), None, None], |x| near(x[0], 0.0));
|
||
b.bc("y0", [None, Some(0.0), None], |x| near(x[1], 0.0));
|
||
b.bc("z0", [None, None, Some(0.0)], |x| near(x[2], 0.0));
|
||
b.bc("push", [None, None, Some(-0.02)], |x| near(x[2], 2.0));
|
||
let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap();
|
||
report(&r);
|
||
(b, r)
|
||
}
|
||
|
||
/// The contact patch test (Taylor & Papadopoulos 1991; Crisfield 2000): a uniform pressure
|
||
/// crosses a non-matching interface exactly, so every element of both blocks has the same
|
||
/// uniaxial stress, with nothing across.
|
||
#[test]
|
||
fn the_mortar_interface_passes_the_contact_patch_test() {
|
||
for (lower, upper, name) in [
|
||
(ElementKind::Hex8, ElementKind::Hex8, "hex8 on hex8"),
|
||
(ElementKind::Hex8, ElementKind::Tet4, "tet4 on hex8"),
|
||
(ElementKind::Tet4, ElementKind::Hex8, "hex8 on tet4"),
|
||
] {
|
||
for geometry in [false, true] {
|
||
// Small sliding is linear and exact; finite sliding is exact up to the gaps the
|
||
// augmentation leaves (their geometric terms, GAP_TOLERANCE) and the convergence.
|
||
let tol = if geometry { 1e-4 } else { 1e-9 };
|
||
let (b, r) = patch(lower, upper, geometry, false);
|
||
let szz: Vec<f64> = r.elements.iter().flat_map(|e| e.nodal.iter().map(|s| s[2])).collect();
|
||
let mean = szz.iter().sum::<f64>() / szz.len() as f64;
|
||
assert!(mean < 0.0, "{name}: {mean}");
|
||
for (e, s) in r.elements.iter().flat_map(|e| e.nodal.iter().map(move |s| (e.label, s))) {
|
||
assert!((s[2] - mean).abs() < tol * mean.abs(), "{name} ({geometry}) element {e}: σzz {} vs {mean}", s[2]);
|
||
assert!(s[0].abs() < tol * mean.abs() && s[1].abs() < tol * mean.abs(), "{name} ({geometry}) element {e}: {s:?}");
|
||
}
|
||
let c = &r.contacts[0];
|
||
// The pressure is uniform, equal to the stress, at every secondary node.
|
||
assert!(c.pressure.iter().all(|p| (p + mean).abs() < tol * mean.abs()), "{name}: {:?} vs {mean}", c.pressure);
|
||
// The contact force on the upper block balances the reaction at its top.
|
||
let top: f64 = b.nodes(|x| near(x[2], 2.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()][2]).sum();
|
||
assert!((c.force[2] + top).abs() < tol * top.abs(), "{name}: {:?} vs {top}", c.force);
|
||
assert!(c.force[0].abs() < tol * top.abs() && c.force[1].abs() < tol * top.abs(), "{name}: {:?}", c.force);
|
||
assert!((c.normal_force + top).abs() < tol * top.abs(), "{name}: {} vs {top}", c.normal_force);
|
||
// Augmented to a gap tolerance far below the penetration a penalty alone leaves.
|
||
assert!(c.penetration < 1e-5, "{name}: {}", c.penetration);
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Hertz contact of a cylinder on a flat (plane strain): a quarter model of an indenter whose
|
||
/// bottom is a circle of radius `r` on a block, symmetric about x = 0, one element thick in z
|
||
/// with z held (plane strain). Returns the model and the result.
|
||
fn hertz_cylinder(finite: bool, push: f64) -> (Blocks, StaticResult) {
|
||
let r = 10.0;
|
||
let mut b = Blocks::default();
|
||
let grade = |s: f64| s.powf(1.6);
|
||
let foundation = b.block(ElementKind::Hex8, [24, 14, 1], |p| [5.0 * grade(p[0]), -5.0 * grade(1.0 - p[1]), 0.2 * p[2]]);
|
||
let indenter = b.block(ElementKind::Hex8, [20, 10, 1], |p| {
|
||
let x = 2.5 * grade(p[0]);
|
||
[x, r - (r * r - x * x).sqrt() + 2.5 * grade(p[1]), 0.2 * p[2]]
|
||
});
|
||
b.define(1000.0, 0.3, NonlinearStatic { nonlinear_geometry: finite, initial_increment: 0.25, ..NonlinearStatic::default() });
|
||
b.face_set("flat", &foundation, |x| near(x[1], 0.0));
|
||
let bottom = |x: [f64; 3]| (x[1] - (r - (r * r - x[0] * x[0]).sqrt())).abs() < 1e-9;
|
||
b.face_set("round", &indenter, bottom);
|
||
b.pair("hertz", "flat", "round");
|
||
b.bc("symmetry", [Some(0.0), None, None], |x| near(x[0], 0.0));
|
||
b.bc("plane", [None, None, Some(0.0)], |_| true);
|
||
b.bc("base", [None, Some(0.0), None], |x| near(x[1], -5.0));
|
||
let top = 2.5 + r - (r * r - 2.5 * 2.5f64).sqrt();
|
||
b.bc("push", [None, Some(-push), None], move |x| {
|
||
x[1] > 2.5 && (x[1] - (r - (r * r - x[0] * x[0]).sqrt()) - 2.5).abs() < 1e-9 && x[1] <= top + 1e-9
|
||
});
|
||
let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap();
|
||
report(&r);
|
||
(b, r)
|
||
}
|
||
|
||
/// Hertz's line contact (H. Hertz, *Über die Berührung fester elastischer Körper*, J. reine
|
||
/// angew. Math. 92:156–171, 1882; K. L. Johnson, *Contact Mechanics*, Cambridge 1985, §4.3):
|
||
/// for a load P per length, b = √(4 P R / (π E*)) and p₀ = 2 P / (π b), E* = E / (2 (1 − ν²))
|
||
/// for two bodies of one material.
|
||
#[test]
|
||
fn hertz_cylinder_on_a_flat() {
|
||
for finite in [false, true] {
|
||
let (b, r) = hertz_cylinder(finite, 0.03);
|
||
let c = &r.contacts[0];
|
||
let thickness = 0.2;
|
||
// The load per length: the reaction at the base over the slab's thickness, both halves.
|
||
let base: f64 = b.nodes(|x| near(x[1], -5.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()][1]).sum();
|
||
let p = 2.0 * base / thickness;
|
||
let e_star = 1000.0 / (2.0 * (1.0 - 0.09));
|
||
let half_width = (4.0 * p * 10.0 / (std::f64::consts::PI * e_star)).sqrt();
|
||
let peak = 2.0 * p / (std::f64::consts::PI * half_width);
|
||
let found = c.pressure.iter().fold(0.0f64, |m, v| m.max(*v));
|
||
let (before, next) = edge(&b, c, 2);
|
||
eprintln!("finite {finite}: P {p:.4} b {half_width:.4} ({before:.4}..{next:.4}) p0 {peak:.4} found {found:.4}");
|
||
assert!((found - peak).abs() < 0.03 * peak, "p0 {found} vs {peak}");
|
||
assert!(before <= half_width && half_width <= next, "b {half_width} not in ({before}, {next})");
|
||
if finite {
|
||
quadratic(&r);
|
||
}
|
||
}
|
||
}
|
||
|
||
/// The contact's edge along the x axis (the nodes on the plane where coordinate `across` is 0):
|
||
/// the node before the last closed one and the first open one after it, so that the edge lies
|
||
/// within an element of the last closed node (whose status the augmentation's tolerance decides).
|
||
fn edge(b: &Blocks, c: &fieldcraft_solver::contact::ContactResult, across: usize) -> (f64, f64) {
|
||
let mut along: Vec<(f64, Status)> =
|
||
c.nodes.iter().zip(&c.status).filter(|(l, _)| near(b.xyz[l][across], 0.0)).map(|(l, s)| (b.xyz[l][0], *s)).collect();
|
||
along.sort_by(|x, y| x.0.total_cmp(&y.0));
|
||
let last = along.iter().rposition(|(_, s)| *s == Status::Closed).unwrap_or(0);
|
||
let before = last.checked_sub(1).and_then(|i| along.get(i)).map_or(0.0, |p| p.0);
|
||
let next = along.get(last + 1).map_or(f64::INFINITY, |p| p.0);
|
||
(before, next)
|
||
}
|
||
|
||
/// Newton's rate with contact: within each solution of an increment (its first, and each
|
||
/// augmentation's), once the residual is small every correction squares it, up to a constant
|
||
/// (the consistent tangent of the mortar integrals, segmentation and normal included).
|
||
fn quadratic(r: &StaticResult) {
|
||
let mut checked = 0;
|
||
for inc in &r.nonlinear.as_ref().unwrap().increments {
|
||
for w in inc.iterations.windows(2) {
|
||
let (a, b) = (w[0], w[1]);
|
||
// Below 1e-10 of the force scale the sums' round-off takes over.
|
||
if a.augmentation != b.augmentation || a.force > 1e-3 || b.force < 1e-10 {
|
||
continue;
|
||
}
|
||
assert!(b.force <= 100.0 * a.force * a.force, "increment {}: {:e} after {:e}", inc.number, b.force, a.force);
|
||
checked += 1;
|
||
}
|
||
}
|
||
assert!(checked >= 3, "too few iterations near the solution to judge ({checked})");
|
||
}
|
||
|
||
/// Two unit blocks, the upper one starting `overlap` into the lower one, pressed together.
|
||
fn stacked(overlap: f64, penetration: InitialPenetration) -> Blocks {
|
||
let mut b = Blocks::default();
|
||
let low = b.block(ElementKind::Hex8, [2, 2, 1], |p| p);
|
||
let up = b.block(ElementKind::Hex8, [3, 3, 1], |p| [p[0], p[1], 1.0 - overlap + p[2]]);
|
||
b.define(1000.0, 0.3, NonlinearStatic { nonlinear_geometry: true, initial_increment: 1.0, ..NonlinearStatic::default() });
|
||
b.face_set("lower_top", &low, |x| near(x[2], 1.0));
|
||
b.face_set("upper_bottom", &up, |x| near(x[2], 1.0 - overlap));
|
||
let mut c = Contact::new(ContactSurfaces::Pair { main: "lower_top".into(), secondary: "upper_bottom".into() });
|
||
c.initial_penetration = penetration;
|
||
b.model.add_contact("interface", c).unwrap();
|
||
b.bc("x0", [Some(0.0), None, None], |x| near(x[0], 0.0));
|
||
b.bc("y0", [None, Some(0.0), None], |x| near(x[1], 0.0));
|
||
b.bc("z0", [None, None, Some(0.0)], |x| near(x[2], 0.0));
|
||
b.bc("push", [None, None, Some(-0.01)], move |x| near(x[2], 2.0 - overlap));
|
||
b
|
||
}
|
||
|
||
#[test]
|
||
fn initial_overclosure_is_adjusted_or_refused() {
|
||
let b = stacked(1e-3, InitialPenetration::Adjust);
|
||
let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap();
|
||
assert!(r.warnings.iter().any(|w| w.contains("start inside the main surface")), "{:?}", r.warnings);
|
||
// Adjusted, the blocks behave as if they started touching: the same stress as without the
|
||
// overlap.
|
||
let touching = solve_nonlinear(&stacked(0.0, InitialPenetration::Adjust).model, Some("Step"), &Limits::default()).unwrap();
|
||
let (a, t) = (r.elements[0].nodal[0][2], touching.elements[0].nodal[0][2]);
|
||
assert!((a - t).abs() < 1e-4 * t.abs(), "{a} vs {t}");
|
||
let e = solve_nonlinear(&stacked(1e-3, InitialPenetration::Refuse).model, Some("Step"), &Limits::default()).unwrap_err().to_string();
|
||
assert!(e.contains("starts") && e.contains("inside"), "{e}");
|
||
}
|
||
|
||
#[test]
|
||
fn contact_results_do_not_depend_on_the_thread_count() {
|
||
let b = stacked(0.0, InitialPenetration::Adjust);
|
||
let run = |threads| fieldcraft_parallel::with_threads(threads, || solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap());
|
||
let (one, many) = (run(1), run(4));
|
||
assert_eq!(one.displacement, many.displacement);
|
||
assert_eq!(one.contacts, many.contacts);
|
||
}
|
||
|
||
#[test]
|
||
fn what_implicit_contact_refuses() {
|
||
let solve = |m: &Model| solve_nonlinear(m, Some("Step"), &Limits::default()).unwrap_err().to_string();
|
||
// Linear analyses refuse every contact by name.
|
||
let mut b = stacked(0.0, InitialPenetration::Adjust);
|
||
b.model.add_step(Step::new("Linear", Procedure::Static {}), None).unwrap();
|
||
let e = fieldcraft_solver::solve_static(&b.model, Some("Linear"), &Limits::default()).unwrap_err().to_string();
|
||
assert!(e.contains("interface") && e.contains("linear static"), "{e}");
|
||
// Arc-length control.
|
||
let mut b = stacked(0.0, InitialPenetration::Adjust);
|
||
let arc = NonlinearStatic { arc_length: Some(ArcLength::default()), ..NonlinearStatic::default() };
|
||
b.model.replace_step(Step { procedure: Procedure::NonlinearStatic(arc), ..b.model.step("Step").unwrap().clone() }).unwrap();
|
||
assert!(solve(&b.model).contains("arc-length"));
|
||
}
|
||
|
||
/// Hertz's point contact (Johnson 1985, §4.2): a sphere of radius R on a flat of the same
|
||
/// material under a load P has the contact radius a = (3 P R / (4 E*))^(1/3) and the peak
|
||
/// pressure p₀ = 3 P / (2 π a²), E* = E / (2 (1 − ν²)). A quarter model, displacement controlled;
|
||
/// the load is the reaction at the base.
|
||
#[test]
|
||
#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"]
|
||
fn hertz_sphere_on_a_flat() {
|
||
let (r, e, nu) = (10.0, 1000.0, 0.3);
|
||
let mut b = Blocks::default();
|
||
let grade = |s: f64| s.powf(1.5);
|
||
let foundation = b.block(ElementKind::Hex8, [16, 16, 10], |p| [5.0 * grade(p[0]), 5.0 * grade(p[1]), -5.0 * grade(1.0 - p[2])]);
|
||
let cap = move |x: f64, y: f64| r - (r * r - x * x - y * y).sqrt();
|
||
let indenter = b.block(ElementKind::Hex8, [14, 14, 8], |p| {
|
||
let (x, y) = (2.5 * grade(p[0]), 2.5 * grade(p[1]));
|
||
let bottom = cap(x, y);
|
||
[x, y, bottom + (3.0 - bottom) * grade(p[2])]
|
||
});
|
||
b.define(e, nu, NonlinearStatic { nonlinear_geometry: false, initial_increment: 0.25, ..NonlinearStatic::default() });
|
||
b.face_set("flat", &foundation, |x| near(x[2], 0.0));
|
||
b.face_set("round", &indenter, move |x| (x[2] - cap(x[0], x[1])).abs() < 1e-9 && x[2] < 1.0);
|
||
b.pair("hertz", "flat", "round");
|
||
b.bc("x0", [Some(0.0), None, None], |x| near(x[0], 0.0));
|
||
b.bc("y0", [None, Some(0.0), None], |x| near(x[1], 0.0));
|
||
b.bc("base", [None, None, Some(0.0)], |x| near(x[2], -5.0));
|
||
b.bc("push", [None, None, Some(-0.04)], |x| near(x[2], 3.0));
|
||
let res = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap();
|
||
report(&res);
|
||
let c = &res.contacts[0];
|
||
let base: f64 = b.nodes(|x| near(x[2], -5.0)).iter().map(|l| res.reaction[res.nodes.binary_search(l).unwrap()][2]).sum();
|
||
let p = 4.0 * base;
|
||
let e_star = e / (2.0 * (1.0 - nu * nu));
|
||
let a = (3.0 * p * r / (4.0 * e_star)).cbrt();
|
||
let peak = 3.0 * p / (2.0 * std::f64::consts::PI * a * a);
|
||
let found = c.pressure.iter().fold(0.0f64, |m, v| m.max(*v));
|
||
let (before, next) = edge(&b, c, 1);
|
||
eprintln!("P {p:.4} a {a:.4} ({before:.4}..{next:.4}) p0 {peak:.4} found {found:.4}");
|
||
assert!((found - peak).abs() < 0.05 * peak, "p0 {found} vs {peak}");
|
||
assert!(before <= a && a <= next, "a {a} not in ({before}, {next})");
|
||
}
|
||
|
||
/// Two blocks, the upper one pressed down by `down` and moved sideways by `side` at its top
|
||
/// (both in proportion over the step), with Coulomb friction `friction`; the reactions at the
|
||
/// top (normal, tangential) and the result.
|
||
fn sheared(down: f64, side: f64, friction: Friction, finite: bool) -> (f64, f64, StaticResult) {
|
||
let mut b = Blocks::default();
|
||
let low = b.block(ElementKind::Hex8, [3, 3, 1], |p| p);
|
||
let up = b.block(ElementKind::Hex8, [2, 2, 1], |p| [p[0], p[1], 1.0 + p[2]]);
|
||
b.define(1000.0, 0.3, NonlinearStatic { nonlinear_geometry: finite, initial_increment: 0.25, ..NonlinearStatic::default() });
|
||
b.face_set("lower_top", &low, |x| near(x[2], 1.0));
|
||
b.face_set("upper_bottom", &up, |x| near(x[2], 1.0));
|
||
let mut c = Contact::new(ContactSurfaces::Pair { main: "lower_top".into(), secondary: "upper_bottom".into() });
|
||
c.friction = friction;
|
||
b.model.add_contact("interface", c).unwrap();
|
||
b.bc("base", [Some(0.0), Some(0.0), Some(0.0)], |x| near(x[2], 0.0));
|
||
b.bc("top", [Some(side), Some(0.0), Some(-down)], |x| near(x[2], 2.0));
|
||
let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap();
|
||
report(&r);
|
||
let sum = |k: usize| -> f64 { b.nodes(|x| near(x[2], 2.0)).iter().map(|l| r.reaction[r.nodes.binary_search(l).unwrap()][k]).sum() };
|
||
(-sum(2), sum(0), r)
|
||
}
|
||
|
||
/// Coulomb's law across a frictional interface: a tangential demand below μ N sticks (the blocks
|
||
/// shear as one, nothing slips), above it the interface slips and carries exactly μ N; with a
|
||
/// kinetic coefficient the slipping interface carries μ_k N.
|
||
#[test]
|
||
fn coulomb_friction_sticks_below_the_cone_and_slips_on_it() {
|
||
for finite in [false, true] {
|
||
let mu = 0.3;
|
||
// Stuck as one body, the shear the top displacement asks for: a frictionless trial gives
|
||
// the normal force, a stuck one with a very large μ the tangential force.
|
||
let (n, t, _) = sheared(0.01, 0.001, Friction::coulomb(5.0), finite);
|
||
let ratio = t / n;
|
||
eprintln!("finite {finite}: stuck N {n:.4} T {t:.4} ratio {ratio:.4}");
|
||
// Well below the cone (the shear concentrates at the edges of the interface, which slip
|
||
// first): sticks, the same forces as glued.
|
||
let below = 0.4 * mu / ratio;
|
||
let (n1, t1, r1) = sheared(0.01, 0.001 * below, Friction::coulomb(mu), finite);
|
||
assert!((t1 / n1 - 0.4 * mu).abs() < 1e-3 * mu, "finite {finite}: {} vs {}", t1 / n1, 0.4 * mu);
|
||
// The middle of the interface sticks (its corners, where the shear concentrates, may slip).
|
||
// (The upper block's 3 × 3 bottom nodes come first, in x then y order: the fifth.)
|
||
assert_eq!(r1.contacts[0].status[4], Status::Stick, "finite {finite}: {:?}", r1.contacts[0].status);
|
||
// Above it: slips, carrying μ N.
|
||
let (n2, t2, r2) = sheared(0.01, 0.001 * 2.0 * mu / ratio, Friction::coulomb(mu), finite);
|
||
// With finite sliding the upper block turns a little as it slides, and the contact acts
|
||
// along its face's normal, a few thousandths of a radian from the vertical.
|
||
let tol = if finite { 0.03 } else { 1e-3 };
|
||
assert!((t2 - mu * n2).abs() < tol * mu * n2, "finite {finite}: T {t2} vs μN {}", mu * n2);
|
||
assert!(r2.contacts[0].status.iter().all(|s| *s == Status::Slip), "finite {finite}: {:?}", r2.contacts[0].status);
|
||
let shear: f64 = r2.contacts[0].shear.iter().map(|s| s[0]).zip(&r2.contacts[0].pressure).map(|(s, p)| s / p).sum::<f64>()
|
||
/ r2.contacts[0].nodes.len() as f64;
|
||
// Node by node the nodal forces' ratio scatters a little about μ (the smoothing of the
|
||
// nodal pressure and shear differs at the interface's edges); the mean is μ.
|
||
assert!((shear + mu).abs() < 0.05 * mu, "finite {finite}: shear over pressure {shear}");
|
||
// Kinetic friction once slipping.
|
||
let decaying = Friction { static_coefficient: mu, kinetic: Some(0.2), decay: 1.0 };
|
||
let (n3, t3, _) = sheared(0.01, 0.001 * 2.0 * mu / ratio, decaying, finite);
|
||
assert!((t3 - 0.2 * n3).abs() < tol * 0.2 * n3, "finite {finite}: T {t3} vs μk N {}", 0.2 * n3);
|
||
}
|
||
}
|
||
|
||
/// Two blocks of non-matching meshes tied at z = 1 (`lower_kind` 3 × 3 × 2, graded or distorted,
|
||
/// below; `upper_kind` 2 × 2 × 2 above), every outer boundary node moved by the linear field
|
||
/// u = G x: the tie patch test.
|
||
fn tied(lower_kind: ElementKind, upper_kind: ElementKind, distorted: bool, procedure: Procedure) -> (Blocks, [[f64; 3]; 3]) {
|
||
let mut b = Blocks::default();
|
||
let low = if distorted {
|
||
b.block(lower_kind, [3, 3, 2], |p| [p[0] + 0.3 * p[0] * (1.0 - p[0]) * (p[1] - 0.5), p[1] + 0.3 * p[1] * (1.0 - p[1]) * (p[0] - 0.4), p[2]])
|
||
} else {
|
||
b.block(lower_kind, [3, 3, 2], |p| [p[0].powf(1.4), p[1].powf(0.8), p[2]])
|
||
};
|
||
let up = b.block(upper_kind, [2, 2, 2], |p| [p[0], p[1], 1.0 + p[2]]);
|
||
let all: Vec<u32> = b.model.elements().keys().copied().collect();
|
||
b.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
|
||
b.model.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: 1000.0, poisson: 0.3 }), ..Material::default() }).unwrap();
|
||
b.model
|
||
.add_section("solid", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() })
|
||
.unwrap();
|
||
b.model.add_step(Step::new("Step", procedure), None).unwrap();
|
||
b.face_set("lower_top", &low, |x| near(x[2], 1.0));
|
||
b.face_set("upper_bottom", &up, |x| near(x[2], 1.0));
|
||
b.model.add_tie("glue", Tie::new("lower_top", "upper_bottom")).unwrap();
|
||
let g = [[1e-3, 2e-4, -3e-4], [4e-4, -5e-4, 1e-4], [-2e-4, 3e-4, 2e-3]];
|
||
let outer = |x: [f64; 3]| near(x[0], 0.0) || near(x[0], 1.0) || near(x[1], 0.0) || near(x[1], 1.0) || near(x[2], 0.0) || near(x[2], 2.0);
|
||
for l in b.nodes(outer) {
|
||
let x = b.xyz[&l];
|
||
let u: [f64; 3] = std::array::from_fn(|i| (0..3).map(|j| g[i][j] * x[j]).sum());
|
||
let name = format!("n{l}");
|
||
b.model.add_set(&name, Set { nodes: vec![l], ..Set::default() }).unwrap();
|
||
b.model
|
||
.add_boundary_condition(
|
||
&name,
|
||
BoundaryCondition::Displacement { target: Target::Set(name.clone()), values: [Some(u[0]), Some(u[1]), Some(u[2]), None, None, None] },
|
||
)
|
||
.unwrap();
|
||
b.model.set_boundary_condition_active("Step", &name, true).unwrap();
|
||
}
|
||
(b, g)
|
||
}
|
||
|
||
/// The tie patch test (dual mortar: Wohlmuth 2000): a linear displacement field across a tie
|
||
/// of non-matching meshes is reproduced exactly, so every element has the same stress, in linear
|
||
/// static, nonlinear static and through the interface nodes the tie moves.
|
||
#[test]
|
||
fn a_tie_passes_the_patch_test() {
|
||
for (lower, upper, distorted, name) in [
|
||
(ElementKind::Hex8, ElementKind::Hex8, false, "hex8 on graded hex8"),
|
||
(ElementKind::Hex8, ElementKind::Tet4, false, "tet4 on hex8"),
|
||
(ElementKind::Tet4, ElementKind::Hex8, false, "hex8 on tet4"),
|
||
(ElementKind::Hex8, ElementKind::Hex8, true, "hex8 on distorted hex8"),
|
||
] {
|
||
let (b, g) = tied(lower, upper, distorted, Procedure::Static {});
|
||
let r = fieldcraft_solver::solve_static(&b.model, Some("Step"), &Limits::default()).unwrap();
|
||
// The interface nodes follow the field too.
|
||
let tol = if distorted { 1e-6 } else { 1e-11 };
|
||
for l in b.nodes(|x| near(x[2], 1.0)) {
|
||
let x = b.xyz[&l];
|
||
let u = r.displacement[r.nodes.binary_search(&l).unwrap()];
|
||
for i in 0..3 {
|
||
let want: f64 = (0..3).map(|j| g[i][j] * x[j]).sum();
|
||
assert!((u[i] - want).abs() < tol * 1e-2, "{name}: node {l} u {u:?}");
|
||
}
|
||
}
|
||
let first = r.elements[0].nodal[0];
|
||
for e in &r.elements {
|
||
for s in &e.nodal {
|
||
for k in 0..6 {
|
||
assert!((s[k] - first[k]).abs() < tol * 10.0, "{name}: element {} {s:?} vs {first:?}", e.label);
|
||
}
|
||
}
|
||
}
|
||
assert!(r.warnings.iter().all(|w| !w.contains("stay free")), "{name}: {:?}", r.warnings);
|
||
}
|
||
// The same in a nonlinear step with small strain.
|
||
let (b, _) = tied(
|
||
ElementKind::Hex8,
|
||
ElementKind::Hex8,
|
||
false,
|
||
Procedure::NonlinearStatic(NonlinearStatic { nonlinear_geometry: false, ..NonlinearStatic::default() }),
|
||
);
|
||
let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap();
|
||
let first = r.elements[0].nodal[0];
|
||
assert!(r.elements.iter().flat_map(|e| &e.nodal).all(|s| (0..6).all(|k| (s[k] - first[k]).abs() < 1e-8)));
|
||
}
|
||
|
||
#[test]
|
||
fn what_ties_refuse() {
|
||
let (mut b, _) = tied(ElementKind::Hex8, ElementKind::Hex8, false, Procedure::Static {});
|
||
b.model.add_step(Step::new("Heat", Procedure::HeatTransfer { transient: None }), None).unwrap();
|
||
let e = fieldcraft_solver::solve_heat(&b.model, Some("Heat"), &Limits::default()).unwrap_err().to_string();
|
||
assert!(e.contains("tie `glue`") && e.contains("heat transfer"), "{e}");
|
||
let (mut b, _) = tied(ElementKind::Hex8, ElementKind::Hex8, false, Procedure::Static {});
|
||
// (Touching surfaces are tied within any tolerance.)
|
||
b.model.replace_tie("glue", Tie { position_tolerance: Some(1e-9), ..Tie::new("lower_top", "upper_bottom") }).unwrap();
|
||
let r = fieldcraft_solver::solve_static(&b.model, Some("Step"), &Limits::default()).unwrap();
|
||
assert!(r.warnings.iter().all(|w| !w.contains("stay free")), "{:?}", r.warnings);
|
||
}
|
||
|
||
/// 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<u32> = 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<u32> = 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);
|
||
}
|