From 98eaede2dd1885f1751233c31d3870606af8cf68 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 17:54:43 -0400 Subject: [PATCH] WIP (parked): folding robustness, stabilization, mesh.revolve, boot-seal profile Parked work for #182/#185 (not for merge as is): - contact: held active set with primal-dual updates replaces the inconsistent ENGAGED_BAND tangent after chattering (moves the nl-friction-sliding-blocks fixture by 4.5e-8; not re-blessed) - nonlinear: opt-in viscous stabilization (step 'stabilization', calibrated on the first increment), Abaqus STABILIZE, LS-DYNA refusal - ordering: graph (level-structure) nested dissection, chosen when it fills less: half the fill on the revolved boot, factor time 62 s -> 25 s - engine: mesh.revolve command (2D profile -> hex8/wedge6, tags kept) - engine: sample/boot.rs boot-seal profile, axisymmetric and 3D scripts Scratch still in the code: FC_TRACE/TIMES timing prints, NO_BAND/STAB/ SQUEEZE env hooks in tests, boot_axisymmetric_scratch test. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/decks/src/abaqus.rs | 11 +- crates/decks/src/lsdyna/mod.rs | 3 + crates/decks/tests/abaqus.rs | 6 + crates/decks/tests/lsdyna.rs | 3 + crates/engine/src/cmd/define.rs | 2 +- crates/engine/src/cmd/forms.rs | 7 + crates/engine/src/cmd/mod.rs | 2 + crates/engine/src/cmd/revolve.rs | 216 +++++++++++++++++++++++++++ crates/engine/src/sample.rs | 43 ++++++ crates/engine/src/sample/boot.rs | 248 +++++++++++++++++++++++++++++++ crates/model/src/steps.rs | 19 +++ crates/solver/src/analysis.rs | 14 +- crates/solver/src/cholesky.rs | 5 + crates/solver/src/contact.rs | 154 ++++++++++++++++--- crates/solver/src/ldlt.rs | 2 + crates/solver/src/nonlinear.rs | 238 +++++++++++++++++++++++++---- crates/solver/src/ordering.rs | 176 ++++++++++++++++++++++ crates/solver/tests/contact.rs | 19 ++- 18 files changed, 1111 insertions(+), 57 deletions(-) create mode 100644 crates/engine/src/cmd/revolve.rs create mode 100644 crates/engine/src/sample/boot.rs diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 010be25..b9bc6f4 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -1240,7 +1240,16 @@ impl Deck { self.out.push('\n'); } } else { - self.line("*STATIC"); + match n.stabilization { + // Abaqus's automatic stabilization takes the dissipated energy + // fraction, which FieldCraft's damping fraction is. + Some(a) => { + self.out.push_str("*STATIC, STABILIZE="); + push_number(&mut self.out, a)?; + self.out.push('\n'); + } + None => self.line("*STATIC"), + } self.numbers_line("", &[n.initial_increment, 1.0, n.min_increment, n.max_increment])?; } if n.line_search { diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 929da99..05b42ae 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -290,6 +290,9 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> { refused.insert(format!("step `{}`: explicit dynamic steps (this writer writes implicit analyses so far)", step.name)); } if let Procedure::NonlinearStatic(n) = &step.procedure { + if n.stabilization.is_some() { + refused.insert(format!("step `{}`: stabilization (artificial damping of a static step; not written yet)", step.name)); + } if !n.nonlinear_geometry { refused.insert(format!( "step `{}`: a nonlinear step without nonlinear geometry (LS-DYNA's nonlinear implicit analysis always follows the deformed shape)", diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs index 74ebf37..d7d0f11 100644 --- a/crates/decks/tests/abaqus.rs +++ b/crates/decks/tests/abaqus.rs @@ -776,6 +776,12 @@ fn nonlinear_steps_and_hyperelastic_materials_are_written() { assert_eq!(bend.get(at + 1..at + 4), Some(&["0.05, 1.0, 0.0001, 0.5", "*CONTROLS, PARAMETERS=LINE SEARCH", "5,"][..]), "{bend:?}"); assert!(deck.contains("*ELASTIC, TYPE=ISOTROPIC\n210000.0, 0.3\n"), "St Venant–Kirchhoff is the elastic material"); + // Stabilization: Abaqus's dissipated energy fraction. + let stabilized = NonlinearStatic { stabilization: Some(2e-4), ..controls.clone() }; + m.replace_step(step("Bend", Procedure::NonlinearStatic(stabilized), &["clamp"], &["end"])).unwrap(); + let deck = abaqus::write(&m).unwrap(); + assert!(block(&deck, "Bend").contains(&"*STATIC, STABILIZE=0.0002"), "{deck}"); + let small = NonlinearStatic { nonlinear_geometry: false, line_search: false, ..controls }; m.replace_step(step("Bend", Procedure::NonlinearStatic(small), &["clamp"], &["end"])).unwrap(); m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::NeoHookean {}), ..steel() }).unwrap(); diff --git a/crates/decks/tests/lsdyna.rs b/crates/decks/tests/lsdyna.rs index 2a86c63..f6a0c32 100644 --- a/crates/decks/tests/lsdyna.rs +++ b/crates/decks/tests/lsdyna.rs @@ -774,6 +774,9 @@ fn nonlinear_static_steps() { .collect(); assert_eq!(curve, [vec![0.0, 0.0], vec![0.5, 1.0], vec![1.0, -1.0]], "{deck}"); + m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic { stabilization: Some(2e-4), ..controls.clone() }), &["clamp"], &["tip_load"])) + .unwrap(); + unsupported(&m, "stabilization"); m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic { nonlinear_geometry: false, ..controls }), &["clamp"], &["tip_load"])) .unwrap(); unsupported(&m, "a nonlinear step without nonlinear geometry"); diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index 59e10cb..e320ad7 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -292,7 +292,7 @@ fn fluid_boundary(s: &mut Session, p: &Value) -> Result { pub const STEP_CREATE: CommandSpec = CommandSpec { id: "model.step.create", label: "New Step", - params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "nonlinear_static", "nonlinear_geometry": true, "initial_increment": 0.1, "min_increment": 1e-5, "max_increment": 1.0 (fractions of the step), "max_increments": 100, "max_iterations": 16 (per increment), "tolerances": {"force": 1e-6, "displacement": 1e-5, "energy": 1e-10}, "line_search": false, "load_path": [[0, 0], [0.5, 1], [1, 0]] (optional: the load factor against the step fraction, piecewise linear from [0, 0] to fraction 1; increments end at its points)} (loads and prescribed displacements are the load factor times their values; it ramps from 0 to 1 without a load path) | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "fluid", "domain": "Flow", "end_time": 1.0 | "steps": 1000, "time_step": (default: chosen for Mach 0.087 and τ ≤ 1), "frames": 10, "magic": 0.1875 (the two-relaxation-time Λ)} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal", "friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb; default frictionless)}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}, "mass_scaling": {"target_increment": 1e-6, "set": "…" (optional, default every element)} (opt-in: elements below the target increment get the density that reaches it; the added mass is reported)}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, + params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "nonlinear_static", "nonlinear_geometry": true, "initial_increment": 0.1, "min_increment": 1e-5, "max_increment": 1.0 (fractions of the step), "max_increments": 100, "max_iterations": 16 (per increment), "tolerances": {"force": 1e-6, "displacement": 1e-5, "energy": 1e-10}, "line_search": false, "load_path": [[0, 0], [0.5, 1], [1, 0]] (optional: the load factor against the step fraction, piecewise linear from [0, 0] to fraction 1; increments end at its points), "stabilization": 2e-4 (optional: viscous damping of local instabilities, folding and snap-through, as the fraction of the step's work it may dissipate; calibrated on the first increment)} (loads and prescribed displacements are the load factor times their values; it ramps from 0 to 1 without a load path) | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "fluid", "domain": "Flow", "end_time": 1.0 | "steps": 1000, "time_step": (default: chosen for Mach 0.087 and τ ≤ 1), "frames": 10, "magic": 0.1875 (the two-relaxation-time Λ)} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal", "friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb; default frictionless)}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}, "mass_scaling": {"target_increment": 1e-6, "set": "…" (optional, default every element)} (opt-in: elements below the target increment get the density that reaches it; the added mass is reported)}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, form: Some(super::forms::STEP_CREATE), kind: CommandKind::Model, run: step_create, diff --git a/crates/engine/src/cmd/forms.rs b/crates/engine/src/cmd/forms.rs index 5af8fb1..b83fe63 100644 --- a/crates/engine/src/cmd/forms.rs +++ b/crates/engine/src/cmd/forms.rs @@ -440,6 +440,7 @@ const PROCEDURE: Field = variants( ), flag("line_search", "Line search").def("false"), json("load_path", "Load path").hint("[[fraction, factor], …]; a ramp when empty"), + positive("stabilization", "Stabilization", D::None).hint("damping fraction (e.g. 2e-4) for folding and snap-through; none when empty"), group( "imperfection", "Imperfection", @@ -755,6 +756,12 @@ pub const PLANAR: &[Field] = &[ int("max_elements", "Max elements", 1, 4_000_000), ]; +pub const REVOLVE: &[Field] = &[ + num("angle", "Angle (°)", D::None).def("360"), + int("sectors", "Sectors", 1, 3600).def("36"), + num("start", "Start (°)", D::None).def("0"), +]; + pub const QUALITY: &[Field] = &[json("elements", "Elements").hint("[labels]; all when empty"), int("bins", "Histogram bins", 1, 1000).def("10")]; pub const SURFACE_CHECK: &[Field] = diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 94a136a..b93d824 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -15,6 +15,7 @@ pub use mesh::analytic_surfaces; mod model; pub(crate) mod octree; mod results; +mod revolve; pub(crate) mod solve; mod tags; @@ -40,6 +41,7 @@ pub static ALL: &[CommandSpec] = &[ mesh::DELAUNAY_POINTS, mesh::SURFACE_CHECK, mesh::PLANAR, + revolve::REVOLVE, mesh::TET_GENERATE, octree::HEX_OCTREE, block::BLOCK, diff --git a/crates/engine/src/cmd/revolve.rs b/crates/engine/src/cmd/revolve.rs new file mode 100644 index 0000000..d514412 --- /dev/null +++ b/crates/engine/src/cmd/revolve.rs @@ -0,0 +1,216 @@ +//! `mesh.revolve`: a 2D mesh of a profile in the XY plane (x the radius, Y the axis, as +//! axisymmetric sections take it) swept about the axis into a solid mesh. +//! +//! The profile's point (r, y) at angle θ becomes (r cos θ, r sin θ, y): the axis becomes Z (up) +//! and the profile's half-plane at θ = `start` the first ring of nodes. Each quad4 becomes a hex8 +//! per sector and each tri3 a wedge6, the profile's node order reversed so the swept elements are +//! positively oriented. The 2D elements and their nodes are replaced (new labels after the +//! largest ones). Geometry tags carry over: an element filling face F of body B fills body B, an +//! element edge on edge E of body B gives the swept faces F named E of body B, and when the +//! sweep is not a whole turn its first and last sides are faces `start` and `end` of the body, +//! so sets named by geometry resolve on the solid. A whole turn joins its last ring to the first. + +use std::collections::BTreeMap; + +use fieldcraft_model::{ElementKind, Entity, GeometryRef}; +use serde_json::{Value, json}; + +use super::mesh::optional_count; +use super::optional_number; +use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; + +pub const REVOLVE: CommandSpec = CommandSpec { + id: "mesh.revolve", + label: "Revolve 2D Mesh", + params: r#"{"angle": 360 (degrees, up to a whole turn), "sectors": 36 (elements around, 1..=3600), "start": 0 (degrees: where the profile's half-plane lies, from +X towards +Y)} — sweeps the 2D mesh (quad4 and tri3 in the XY plane, x ≥ 0 the radius, Y the axis) about the axis into hex8 and wedge6 with the axis along Z: (r, y) at angle θ becomes (r cos θ, r sin θ, y). Replaces the 2D elements and their nodes; geometry tags carry over (a face's elements fill its body, an edge's element edges sweep to faces of that name, and a partial turn's first and last sides are faces `start` and `end`), so sets named by geometry resolve on the solid"#, + form: Some(super::forms::REVOLVE), + kind: CommandKind::Model, + run: revolve, +}; + +/// Most sectors of a sweep. +const MAX_SECTORS: u64 = 3600; + +/// The hex8 face (1-based) swept by each quad4 edge, and the wedge6 face swept by each tri3 +/// edge, with the profile's nodes reversed (see [`swept`]); faces 1 and 2 are the first and +/// last sides of a sector. +const QUAD_EDGE_FACE: [u8; 4] = [6, 5, 4, 3]; +const TRI_EDGE_FACE: [u8; 3] = [5, 4, 3]; + +/// A swept element's nodes: the profile's nodes reversed on ring `k`, then on ring `k + 1`. +fn swept(profile: &[u32], ring: impl Fn(u32, usize) -> u32, k: usize) -> Vec { + let reversed: Vec = profile.first().copied().into_iter().chain(profile.iter().skip(1).rev().copied()).collect(); + reversed.iter().map(|n| ring(*n, k)).chain(reversed.iter().map(|n| ring(*n, k + 1))).collect() +} + +fn revolve(s: &mut Session, p: &Value) -> Result { + let angle = optional_number(p, "angle", 360.0)?; + if !(angle > 0.0 && angle <= 360.0) { + return Err(EngineError::Params(format!("`angle` must be in (0, 360] degrees, got {angle}"))); + } + let start = optional_number(p, "start", 0.0)?; + if !start.is_finite() { + return Err(EngineError::Params("`start` must be a number of degrees".into())); + } + let sectors = optional_count(p, "sectors", 36, 1..=MAX_SECTORS)?; + let whole = angle >= 360.0; + if whole && sectors < 3 { + return Err(EngineError::Params("a whole turn takes at least 3 sectors".into())); + } + // The profile: every element a linear 2D element in the XY plane, off the axis. + let profile: Vec<(u32, ElementKind, Vec)> = s.model.elements().iter().map(|(l, e)| (*l, e.kind, e.nodes.clone())).collect(); + if profile.is_empty() { + return Err(EngineError::Params("there is no 2D mesh to revolve".into())); + } + if let Some((l, k, _)) = profile.iter().find(|(_, k, _)| !matches!(k, ElementKind::Quad4 | ElementKind::Tri3)) { + return Err(EngineError::Params(format!( + "not supported yet: revolving {} elements (element {l}): mesh.revolve takes a mesh of quad4 and tri3 elements", + k.name() + ))); + } + let used: Vec = { + let mut v: Vec = profile.iter().flat_map(|(_, _, n)| n.iter().copied()).collect(); + v.sort_unstable(); + v.dedup(); + v + }; + let mut points = Vec::with_capacity(used.len()); + for n in &used { + let xyz = s.model.node(*n).map(|n| n.xyz).ok_or_else(|| EngineError::Params(format!("node {n} is missing")))?; + if xyz[2] != 0.0 { + return Err(EngineError::Params(format!("node {n} is off the XY plane (z = {})", xyz[2]))); + } + if !(xyz[0] > 0.0) { + return Err(EngineError::Params(format!( + "not supported yet: revolving a profile that reaches the axis (node {n} at x = {}): its elements would collapse; keep the profile at x > 0", + xyz[0] + ))); + } + points.push(xyz); + } + let rings = if whole { sectors } else { sectors + 1 }; + let count = |a: usize, b: usize| a.checked_mul(b).ok_or_else(|| EngineError::Params("the revolved mesh is too large".into())); + let (node_count, element_count) = (count(used.len(), rings)?, count(profile.len(), sectors)?); + if element_count > fieldcraft_model::MAX_ELEMENTS { + return Err(EngineError::Params(format!("the revolved mesh would have {element_count} elements, more than {}", fieldcraft_model::MAX_ELEMENTS))); + } + // The tags to carry over, read before the profile is removed. + let tags = s.model.classification().clone(); + let entity = |tag: u32| tags.entity(tag).cloned(); + let (first_node, first_element) = super::fresh_labels(&s.model, node_count, element_count)?; + let index: BTreeMap = used.iter().enumerate().map(|(i, n)| (*n, i)).collect(); + let ring = |n: u32, k: usize| -> u32 { + let i = index.get(&n).copied().unwrap_or(0); + first_node + ((k % rings) * used.len() + i) as u32 + }; + let (start_rad, step) = (start.to_radians(), angle.to_radians() / sectors as f64); + for k in 0..rings { + let (sin, cos) = (start_rad + step * k as f64).sin_cos(); + for (i, x) in points.iter().enumerate() { + s.model.set_node(first_node + (k * used.len() + i) as u32, [x[0] * cos, x[0] * sin, x[1]])?; + } + } + let mut label = first_element; + let mut new_tags: Vec<(u32, u8, GeometryRef)> = Vec::new(); + let mut regions: Vec<(u32, GeometryRef)> = Vec::new(); + for k in 0..sectors { + for (old, kind, nodes) in &profile { + let (solid, edge_face) = match kind { + ElementKind::Quad4 => (ElementKind::Hex8, &QUAD_EDGE_FACE[..]), + _ => (ElementKind::Wedge6, &TRI_EDGE_FACE[..]), + }; + s.model.set_element(label, solid, swept(nodes, ring, k))?; + if let Some(t) = tags.elements().get(old) { + if let Some(region) = t.region.and_then(entity) { + let body = GeometryRef::body(region.body.clone())?; + regions.push((label, body.clone())); + if !whole && k == 0 { + new_tags.push((label, 1, GeometryRef::face(body.body.clone(), "start".to_string())?)); + } + if !whole && k + 1 == sectors { + new_tags.push((label, 2, GeometryRef::face(body.body.clone(), "end".to_string())?)); + } + } + for (edge, tag) in &t.edges { + let (Some(e), Some(face)) = (entity(*tag), edge_face.get(usize::from(*edge).saturating_sub(1))) else { continue }; + if let Entity::Edge(name) = &e.entity { + new_tags.push((label, *face, GeometryRef::face(e.body.clone(), name.clone())?)); + } + } + } + label += 1; + } + } + // The profile goes; its tags with it. + for (old, _, _) in &profile { + s.model.remove_element(*old)?; + } + for n in &used { + s.model.remove_node(*n)?; + } + for (element, body) in regions { + let tag = s.model.geometry_tag(&body)?; + s.model.tag_element(element, tag)?; + } + for (element, face, reference) in new_tags { + let tag = s.model.geometry_tag(&reference)?; + s.model.tag_element_face(element, face, tag)?; + } + Ok(json!({ + "nodes": node_count, + "elements": element_count, + "first_node": first_node, + "first_element": first_element, + "sectors": sectors, + "angle": angle, + })) +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A ring of rectangular section (1 ≤ r ≤ 2, 0 ≤ y ≤ 1) as two quads and two triangles, + /// revolved: its volume is π (2² − 1²) · 1 = 3π, every element valid, the tags carried. + #[test] + fn a_revolved_profile_is_a_valid_solid_of_the_right_volume() { + for (angle, sectors) in [(360.0f64, 48usize), (90.0, 12)] { + let mut s = Session::default(); + s.execute("mesh.planar", &json!({"size": 0.5, "polygon": [[1, 0], [2, 0], [2, 1], [1, 1]], "name": "Ring"})).unwrap(); + let profile = s.model.elements().len(); + let r = s.execute("mesh.revolve", &json!({"angle": angle, "sectors": sectors})).unwrap(); + assert_eq!(r["elements"].as_u64(), Some((profile * sectors) as u64)); + let q = s.execute("mesh.quality", &json!({})).unwrap(); + assert_eq!(q["invalid"]["count"].as_u64(), Some(0), "{q}"); + let volume: f64 = q["volume"].as_f64().unwrap(); + // Chords: a polygon of `sectors` sides over the arc. + let chord = (angle.to_radians() / sectors as f64).sin() * sectors as f64 / angle.to_radians(); + let exact = 3.0 * std::f64::consts::PI * angle / 360.0; + assert!((volume - exact * chord).abs() < 1e-9 * exact, "{volume} {}", exact * chord); + // The bottom edge (e1, y = 0) sweeps to face e1: a set on it holds the z = 0 faces. + s.execute("model.set.create", &json!({"name": "bottom", "geometry": [{"body": "Ring", "face": "e1"}]})).unwrap(); + let set = s.execute("model.set.resolve", &json!({"name": "bottom"})).unwrap(); + assert_eq!(set["faces"]["count"].as_u64(), Some(2 * sectors as u64), "two element edges along e1: {set}"); + if angle < 360.0 { + s.execute("model.set.create", &json!({"name": "cut", "geometry": [{"body": "Ring", "face": "start"}]})).unwrap(); + let set = s.execute("model.set.resolve", &json!({"name": "cut"})).unwrap(); + assert_eq!(set["faces"]["count"].as_u64(), Some(profile as u64), "{set}"); + } + } + } + + #[test] + fn what_mesh_revolve_refuses() { + let mut s = Session::default(); + assert!(s.execute("mesh.revolve", &json!({})).is_err(), "no mesh"); + s.execute("mesh.planar", &json!({"size": 0.5, "polygon": [[0, 0], [2, 0], [2, 1], [0, 1]]})).unwrap(); + let e = s.execute("mesh.revolve", &json!({})).unwrap_err().to_string(); + assert!(e.contains("axis"), "{e}"); + let mut s = Session::default(); + s.execute("mesh.planar", &json!({"size": 0.5, "kind": "quad8", "polygon": [[1, 0], [2, 0], [2, 1], [1, 1]]})).unwrap(); + assert!(s.execute("mesh.revolve", &json!({})).unwrap_err().to_string().contains("quad8")); + for bad in [json!({"angle": 0}), json!({"angle": 400}), json!({"sectors": 0}), json!({"angle": 360, "sectors": 2}), json!({"start": "x"})] { + assert!(s.execute("mesh.revolve", &bad).is_err(), "{bad}"); + } + } +} diff --git a/crates/engine/src/sample.rs b/crates/engine/src/sample.rs index e49913c..e8da7ef 100644 --- a/crates/engine/src/sample.rs +++ b/crates/engine/src/sample.rs @@ -15,6 +15,8 @@ //! geometry is a SolidCraft script run when the sample opens ([`Sample::design`]): a corpus //! recipe (`tests/corpus`) or [`UPRIGHT_BRACKET`]. +pub mod boot; + use fieldcraft_geometry::Shape; use serde_json::{Value, json}; @@ -792,4 +794,45 @@ mod tests { let kept = pz(last) / pz(0); assert!(kept > 0.4 && kept < 0.6 && kinetic(last) < kinetic(0), "{} {}", pz(last), kinetic(last)); } + + #[test] + #[ignore = "scratch"] + fn boot_axisymmetric_scratch() { + let travel: f64 = std::env::var("TRAVEL").ok().and_then(|v| v.parse().ok()).unwrap_or(30.0); + let mut s = Session::default(); + let script = match std::env::var("SECTORS").ok().and_then(|v| v.parse::().ok()) { + Some(n) => super::boot::solid(travel, n, std::env::var("FULL").is_err()), + None => super::boot::axisymmetric(travel), + }; + let t0 = std::time::Instant::now(); + for c in script.as_array().unwrap() { + s.execute(c["command"].as_str().unwrap(), &c["params"]).unwrap(); + } + eprintln!("built {} elements in {:?}, height {}", s.model.elements().len(), t0.elapsed(), super::boot::height()); + if let Ok(path) = std::env::var("DUMP_MODEL") { + std::fs::write(path, serde_json::to_string(&s.model).unwrap()).unwrap(); + } + if std::env::var("NO_CONTACT").is_ok() { + s.execute("interaction.contact.delete", &json!({"name": "Folds"})).unwrap(); + } + let check = s.execute("model.check", &json!({})).unwrap(); + eprintln!("{check}"); + let q = s.execute("mesh.quality", &json!({})).unwrap(); + eprintln!("quality {}", q.to_string().chars().take(300).collect::()); + let t0 = std::time::Instant::now(); + let r = s.execute("solve.nonlinear", &json!({})); + eprintln!("solve {:?}: {}", t0.elapsed(), match &r { Ok(v) => v.to_string().chars().take(3000).collect::(), Err(e) => e.to_string() }); + if r.is_ok() { + let labels: Vec = s.model.nodes().keys().copied().collect(); + let mut out = Vec::new(); + for chunk in labels.chunks(5000) { + let p = s.execute("solve.probe", &json!({"step": "Squeeze", "nodes": chunk})).unwrap(); + for (l, n) in chunk.iter().zip(p["nodes"].as_array().unwrap()) { + out.push(json!({"label": l, "xyz": s.model.nodes()[l].xyz, "u": n["displacement"], "p": n["contact"][0]["pressure"]})); + } + } + let elements: Vec = s.model.elements().iter().map(|(l, e)| json!({"label": l, "nodes": e.nodes})).collect(); + std::fs::write("/private/tmp/claude-501/fc/boot_state.json", json!({"nodes": out, "elements": elements}).to_string()).unwrap(); + } + } } diff --git a/crates/engine/src/sample/boot.rs b/crates/engine/src/sample/boot.rs new file mode 100644 index 0000000..e1f0d98 --- /dev/null +++ b/crates/engine/src/sample/boot.rs @@ -0,0 +1,248 @@ +//! The boot seal: a conical rubber bellows of our own design, its six rounded convolutions +//! squeezed axially by a rigid collar until they fold onto each other. +//! +//! **Profile.** The wall's mid-surface in the (r, y) half-plane, y the axis: a base cylinder +//! (r = [`BASE_RADIUS`]), then six V-shaped convolutions zig-zagging out to a crest and back in to +//! a root, the crests and roots tapering linearly towards the top, and the top collar, a cylinder +//! at the last root's radius. Every corner is rounded with radius [`BEND`]. The wall is the +//! mid-surface offset by ± [`WALL`]/2 along its normal: a structured quad mesh, [`LAYERS`] +//! elements through the wall and about [`STEP`] mm along it. + +use serde_json::{Value, json}; + +/// Mid-surface radius of the base cylinder (mm). +pub const BASE_RADIUS: f64 = 27.0; +/// Wall thickness (mm). +pub const WALL: f64 = 1.8; +/// Radius of the rounded corners of the mid-surface (mm). +pub const BEND: f64 = 1.6; +/// Convolutions, and their pitch along the axis (mm). +pub const FOLDS: usize = 6; +const PITCH: f64 = 9.5; +/// Height of the straight base, and the boot's height (mm). +const BASE_LENGTH: f64 = 4.0; +pub const HEIGHT: f64 = 70.0; +/// Crest and root radii of the first and last convolutions (mm). +const CREST: [f64; 2] = [38.0, 26.0]; +const ROOT: [f64; 2] = [28.0, 17.0]; +/// Elements through the wall and the target element length along it (mm). +pub const LAYERS: usize = 3; +pub const STEP: f64 = 0.6; + +/// A piece of the mid-surface: a straight line, or an arc about `center` of `radius` from angle +/// `start` through `sweep` (radians, counter-clockwise positive). +#[derive(Debug, Clone, Copy)] +enum Piece { + Line { from: [f64; 2], to: [f64; 2] }, + Arc { center: [f64; 2], radius: f64, start: f64, sweep: f64 }, +} + +impl Piece { + fn length(&self) -> f64 { + match *self { + Piece::Line { from, to } => (to[0] - from[0]).hypot(to[1] - from[1]), + Piece::Arc { radius, sweep, .. } => radius * sweep.abs(), + } + } + + /// The point and unit tangent at fraction `s` of the piece. + fn at(&self, s: f64) -> ([f64; 2], [f64; 2]) { + match *self { + Piece::Line { from, to } => { + let l = self.length(); + ([from[0] + s * (to[0] - from[0]), from[1] + s * (to[1] - from[1])], [(to[0] - from[0]) / l, (to[1] - from[1]) / l]) + } + Piece::Arc { center, radius, start, sweep } => { + let a = start + s * sweep; + let (sin, cos) = a.sin_cos(); + let turn = sweep.signum(); + ([center[0] + radius * cos, center[1] + radius * sin], [-turn * sin, turn * cos]) + } + } + } + +} + +/// The crest or root radius of convolution `k` (linear from the first to the last). +fn taper(ends: [f64; 2], k: usize) -> f64 { + let f = k as f64 / (FOLDS - 1) as f64; + ends[0] + f * (ends[1] - ends[0]) +} + +/// The corners of the mid-surface's polyline, base to collar top. +fn corners() -> Vec<[f64; 2]> { + let mut out = vec![[BASE_RADIUS, 0.0], [BASE_RADIUS, BASE_LENGTH]]; + for k in 0..FOLDS { + let y = BASE_LENGTH + k as f64 * PITCH; + out.push([taper(CREST, k), y + 0.5 * PITCH]); + out.push([taper(ROOT, k), y + PITCH]); + } + let collar = out.last().map_or(0.0, |c| c[0]); + out.push([collar, HEIGHT]); + out +} + +/// The mid-surface, base to collar: the polyline of [`corners`] with each inner corner rounded +/// by an arc of radius [`BEND`] tangent to both sides (`None` when a side is too short for its +/// arcs). +fn mid_surface() -> Option> { + let c = corners(); + let unit = |a: [f64; 2], b: [f64; 2]| { + let l = (b[0] - a[0]).hypot(b[1] - a[1]); + [(b[0] - a[0]) / l, (b[1] - a[1]) / l] + }; + let mut out = Vec::new(); + let mut from = *c.first()?; + for w in c.windows(3) { + let (a, v, b) = (w[0], w[1], w[2]); + let (d1, d2) = (unit(a, v), unit(v, b)); + let cross = d1[0] * d2[1] - d1[1] * d2[0]; + let turn = cross.atan2(d1[0] * d2[0] + d1[1] * d2[1]); + let tangent = BEND * (0.5 * turn.abs()).tan(); + let enter = [v[0] - tangent * d1[0], v[1] - tangent * d1[1]]; + // The side must hold both of its arcs' tangent lengths. + if tangent > (v[0] - from[0]).hypot(v[1] - from[1]) { + return None; + } + out.push(Piece::Line { from, to: enter }); + // The centre on the inside of the turn: left of d1 for a left (counter-clockwise) turn. + let side = turn.signum(); + let center = [enter[0] - side * BEND * d1[1], enter[1] + side * BEND * d1[0]]; + let start = (enter[1] - center[1]).atan2(enter[0] - center[0]); + out.push(Piece::Arc { center, radius: BEND, start, sweep: turn }); + from = [v[0] + tangent * d2[0], v[1] + tangent * d2[1]]; + } + out.push(Piece::Line { from, to: *c.last()? }); + Some(out) +} + +/// The mesh of the wall in the (r, y) half-plane: nodes row by row along the mid-surface +/// (`LAYERS + 1` per row, from the outer face in), and the number of rows. +pub fn wall_nodes() -> (Vec<[f64; 2]>, usize) { + let mut rows: Vec<([f64; 2], [f64; 2])> = Vec::new(); + for piece in mid_surface().unwrap_or_default() { + let n = (piece.length() / STEP).round().max(1.0) as usize; + let first = if rows.is_empty() { 0 } else { 1 }; + for i in first..=n { + rows.push(piece.at(i as f64 / n as f64)); + } + } + let mut nodes = Vec::with_capacity(rows.len() * (LAYERS + 1)); + for (p, t) in &rows { + // The left normal of the tangent: through the wall from the outer face in, so the quads + // run counter-clockwise. + let n = [-t[1], t[0]]; + for j in 0..=LAYERS { + let d = WALL * (j as f64 / LAYERS as f64 - 0.5); + nodes.push([p[0] + d * n[0], p[1] + d * n[1]]); + } + } + (nodes, rows.len()) +} + +/// The top of the collar (mm): the boot's height. +pub fn height() -> f64 { + HEIGHT +} + +/// The wall's quad4 mesh in the (r, y) = (x, y) plane, tagged as face `wall` of body `Boot` +/// with its edges `base`, `collar` (the ends), `outer` and `inner` (the faces of the wall). +fn profile(steps: &mut Vec) { + let (nodes, rows) = wall_nodes(); + let per = LAYERS + 1; + let label = |i: usize, j: usize| (i * per + j + 1) as u64; + for (k, p) in nodes.iter().enumerate() { + steps.push(json!({"command": "mesh.node.add", "params": {"label": k + 1, "xyz": [p[0], p[1], 0.0]}})); + } + let mut element = 0u64; + let mut edges: [Vec; 4] = Default::default(); + for i in 0..rows.saturating_sub(1) { + for j in 0..LAYERS { + element += 1; + let quad = [label(i, j), label(i + 1, j), label(i + 1, j + 1), label(i, j + 1)]; + steps.push(json!({"command": "mesh.element.add", "params": {"label": element, "kind": "quad4", "nodes": quad}})); + // Edge 1 runs along the outer face (j = 0), 3 along the inner one, 4 and 2 across + // the wall at the base and the collar. + let mut on = |k: usize, edge: u8| edges[k].push(json!({"element": element, "edge": edge})); + if i == 0 { + on(0, 4); + } + if i + 2 == rows { + on(1, 2); + } + if j == 0 { + on(2, 1); + } + if j + 1 == LAYERS { + on(3, 3); + } + } + } + steps.push(json!({"command": "mesh.tag", "params": {"geometry": {"body": "Boot", "face": "wall"}, "elements": (1..=element).collect::>()}})); + for (name, list) in ["base", "collar", "outer", "inner"].iter().zip(edges) { + steps.push(json!({"command": "mesh.tag", "params": {"geometry": {"body": "Boot", "edge": name}, "edges": list}})); + } +} + +/// The rubber, its section, general contact on every exterior face or edge, the clamped base +/// and the collar pushed down by `travel` along the axis (`axis`: 1 for Y in 2D, 2 for Z in 3D) +/// and held across it (a rigid collar). +fn analysis(steps: &mut Vec, section: Value, entity: &str, axis: usize, travel: f64, extra: &[&str]) { + // The elements: the profile's face, or the solid body. + let region = if axis == 1 { json!({"body": "Boot", "face": "wall"}) } else { json!({"body": "Boot"}) }; + let mut push = [json!(0), json!(0), json!(0), Value::Null, Value::Null, Value::Null]; + if axis == 1 { + push[2] = Value::Null; + } + push[axis] = json!(-travel); + let mut held = vec!["Base", "Push"]; + held.extend_from_slice(extra); + steps.extend([ + json!({"command": "model.set.create", "params": {"name": "boot", "geometry": [region]}}), + json!({"command": "model.set.create", "params": {"name": "base", "geometry": [{"body": "Boot", entity: "base"}]}}), + json!({"command": "model.set.create", "params": {"name": "collar", "geometry": [{"body": "Boot", entity: "collar"}]}}), + json!({"command": "model.material.create", "params": {"name": "Rubber", "elastic": {"type": "isotropic", "young": RUBBER[0], "poisson": RUBBER[1]}, "density": 1.1e-9, + "hyperelastic": {"type": "mooney_rivlin", "c01": RUBBER[2]}}}), + json!({"command": "model.section.create", "params": {"name": "Boot", "set": "boot", "material": "Rubber", "kind": section, "formulation": {"hybrid": true}}}), + json!({"command": "interaction.contact.general", "params": {"name": "Folds"}}), + json!({"command": "model.boundary_condition.create", "params": {"name": "Base", "type": "fixed", "target": {"set": "base"}}}), + json!({"command": "model.boundary_condition.create", "params": {"name": "Push", "type": "displacement", "target": {"set": "collar"}, "values": push}}), + json!({"command": "model.step.create", "params": {"name": "Squeeze", "procedure": {"type": "nonlinear_static", "initial_increment": 0.02, "max_increment": 0.05, "max_increments": 1000, "line_search": true, "stabilization": STABILIZATION}, "boundary_conditions": held}}), + ]); +} + +/// The rubber: Young's modulus (MPa), Poisson's ratio and Mooney–Rivlin C01 (MPa), so +/// μ = E / (2 (1 + ν)) = 1.0 MPa, C10 = μ/2 − C01 = 0.4 MPa: a soft generic rubber, nearly +/// incompressible (κ/μ ≈ 830). +pub const RUBBER: [f64; 3] = [3.0, 0.499, 0.1]; + +/// The squeeze's stabilization (the fraction of its work the damping may dissipate). +pub const STABILIZATION: f64 = 2e-3; + +/// The axisymmetric boot (hybrid quad4) squeezed by `travel` mm. +pub fn axisymmetric(travel: f64) -> Value { + let mut steps = vec![json!({"command": "model.units", "params": {"system": "mm-t-s-N-MPa"}})]; + profile(&mut steps); + analysis(&mut steps, json!({"type": "axisymmetric"}), "edge", 1, travel, &[]); + Value::Array(steps) +} + +/// The boot in 3D: the axisymmetric wall mesh revolved about the axis (`mesh.revolve`, Z up) +/// into hybrid hex8, `sectors` around the circumference, or around the half x ≥ 0 (from −90° to +/// 90°) with the cut held on the symmetry plane x = 0. +pub fn solid(travel: f64, sectors: usize, half: bool) -> Value { + let mut steps = vec![json!({"command": "model.units", "params": {"system": "mm-t-s-N-MPa"}})]; + profile(&mut steps); + let (angle, start) = if half { (180, -90) } else { (360, 0) }; + steps.push(json!({"command": "mesh.revolve", "params": {"angle": angle, "sectors": sectors, "start": start}})); + let mut extra = Vec::new(); + if half { + steps.extend([ + json!({"command": "model.set.create", "params": {"name": "cut", "geometry": [{"body": "Boot", "face": "start"}, {"body": "Boot", "face": "end"}]}}), + json!({"command": "model.boundary_condition.create", "params": {"name": "Symmetry", "type": "symmetry", "target": {"set": "cut"}, "plane": "x"}}), + ]); + extra.push("Symmetry"); + } + analysis(&mut steps, json!({"type": "solid"}), "face", 2, travel, &extra); + Value::Array(steps) +} diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index ae456d3..69de3fb 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -155,6 +155,13 @@ pub struct NonlinearStatic { /// lengths in proportion. `None`: load control. #[serde(skip_serializing_if = "Option::is_none")] pub arc_length: Option, + /// Viscous stabilization of local instabilities (wrinkling, folding, snap-through of a + /// part under displacement control): damping forces α D (u − uₙ)/Δt against each + /// increment's motion, with D the diagonal of the step's first tangent and Δt the increment + /// as a fraction of the step, so α is about the fraction of the forces the damping takes + /// where the motion is smooth. `None`: none. + #[serde(skip_serializing_if = "Option::is_none")] + pub stabilization: Option, } /// The constraint of an arc-length step. @@ -242,6 +249,7 @@ impl Default for NonlinearStatic { load_path: None, imperfection: None, arc_length: None, + stabilization: None, } } } @@ -273,6 +281,14 @@ impl NonlinearStatic { } finite("imperfection scale", imp.scale)?; } + if let Some(a) = self.stabilization { + if !(a > 0.0 && a <= MAX_STABILIZATION) { + return Err(Error::Invalid(format!("the stabilization is a damping fraction in (0, {MAX_STABILIZATION}], got {a}"))); + } + if self.arc_length.is_some() { + return Err(Error::Invalid("an arc-length step follows instabilities along its path; it takes no stabilization".into())); + } + } let initial = fraction("initial increment", self.initial_increment)?; let min = fraction("minimum increment", self.min_increment)?; let max = fraction("maximum increment", self.max_increment)?; @@ -331,6 +347,9 @@ impl NonlinearStatic { } } +/// The largest damping fraction of a nonlinear step's stabilization. +pub const MAX_STABILIZATION: f64 = 1.0; + /// Most time increments of a transient heat transfer step. pub const MAX_HEAT_INCREMENTS: f64 = 100_000.0; diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 42804ce..120b8a3 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -15,7 +15,7 @@ use crate::cholesky::{Factor, FactorError, UpperColumns}; use crate::constraints::{Budget, Reduction, coupling_rows}; use crate::eigen::{Operator, Wanted, lanczos}; use crate::ldlt::{Kind, Symmetric}; -use crate::ordering::nested_dissection; +use crate::ordering::{fill, graph_dissection, nested_dissection}; use crate::sparse::{Csr, NodeGraph, add_block, add_diagonal, add_scalar}; use crate::system::{ Mesh, Structural, StructuralKernel, add_thermal_loads, inelastic_warnings, loads, prescribed, structural_body_forces, temperature_loads, @@ -311,12 +311,12 @@ pub(crate) fn pattern(mesh: &Mesh) -> (NodeGraph, Vec) { let graph = NodeGraph::new(mesh.block_count(), mesh.blocks()); // The order covers the mesh nodes; reference point equations follow them. let n = mesh.labels.len(); - let order = if graph.adjacency.len() == n { - nested_dissection(&mesh.xyz, &graph.adjacency) - } else { - let mesh_part: Vec> = graph.adjacency.iter().take(n).map(|row| row.iter().copied().filter(|b| *b < n).collect()).collect(); - nested_dissection(&mesh.xyz, &mesh_part) - }; + let mesh_part: Vec> = graph.adjacency.iter().take(n).map(|row| row.iter().copied().filter(|b| *b < n).collect()).collect(); + // The dissection that fills the factor less: by coordinate planes, or by the graph's own + // level structures (better on thin curved walls). + let by_planes = nested_dissection(&mesh.xyz, &mesh_part); + let by_graph = graph_dissection(&mesh_part); + let order = if fill(&mesh_part, &by_graph) < fill(&mesh_part, &by_planes) { by_graph } else { by_planes }; (graph, order) } diff --git a/crates/solver/src/cholesky.rs b/crates/solver/src/cholesky.rs index 6b2bc73..0b675d9 100644 --- a/crates/solver/src/cholesky.rs +++ b/crates/solver/src/cholesky.rs @@ -152,6 +152,11 @@ fn ereach(a: &UpperColumns, k: usize, parent: &[usize], stack: &mut [usize], mar } /// The elimination tree and the entry count of each column of L, without computing values. +/// The entries of the Cholesky factor of `a`'s pattern (in its own order). +pub(crate) fn factor_entries(a: &UpperColumns) -> usize { + symbolic(a).1.iter().sum() +} + fn symbolic(a: &UpperColumns) -> (Vec, Vec) { let parent = etree(a); let mut counts = vec![1usize; a.n]; diff --git a/crates/solver/src/contact.rs b/crates/solver/src/contact.rs index 512a949..2c2671d 100644 --- a/crates/solver/src/contact.rs +++ b/crates/solver/src/contact.rs @@ -119,19 +119,14 @@ pub const PENALTY_FACTOR: f64 = 10.0; /// the element under it (and no open node keeps a multiplier). pub const GAP_TOLERANCE: f64 = 1e-5; -/// Once an increment's active set has chattered, a node open by less than this fraction of its -/// element's thickness keeps its penalty in the tangent for the rest of the step (its -/// pressure, and so the residual, is still zero): a node that touches with no pressure, the edge -/// of a contact that lifts, otherwise flips between open and closed from one Newton iteration to -/// the next, each state's correction pushing it into the other, and a smaller increment does not -/// help under proportional loading. -pub const ENGAGED_BAND: f64 = 1e-3; +/// Most updates of a held active set in one increment. +pub const MAX_ACTIVE_SET_UPDATES: usize = 8; /// Most augmentations of one increment. pub const MAX_AUGMENTATIONS: usize = 25; -/// An increment is cut back when the active set alternates between two states this many times -/// (and the step takes [`ENGAGED_BAND`] from then on). +/// An increment whose active set alternates between two states this many times holds it fixed +/// ([`Committed::active`]); if it still alternates once released, the increment is cut back. pub const CHATTER_LIMIT: usize = 3; /// Most faces over the contacts of a model. @@ -187,6 +182,10 @@ pub struct Committed { pub multipliers: Vec, pub traction: Vec, pub slipping: Vec, + /// An active set held fixed while an increment that chattered converges: per secondary + /// node, closed (its pressure λ − ε ĝ, tension allowed) or open (none). `None`: each + /// evaluation decides, p = ⟨λ − ε ĝ⟩. + pub active: Option>, } impl Committed { @@ -793,13 +792,13 @@ impl Interface { multipliers: vec![0.0; n], traction: if friction { vec![[0.0; 3]; n] } else { Vec::new() }, slipping: vec![false; if friction { n } else { 0 }], + active: None, } } /// The contact's forces at displacements `u` (added to `internal`), its tangent blocks when - /// `tangent`, and its state, from the last converged increment's state `committed` and - /// displacements `reference`; nodes open by less than `band` of their element's thickness - /// keep their penalty in the tangent ([`ENGAGED_BAND`]; 0: only closed nodes). + /// `tangent`, and its state, from the last converged increment's state `committed` (with + /// its active set when one is held fixed) and displacements `reference`. #[allow(clippy::too_many_arguments)] pub fn evaluate( &self, @@ -807,7 +806,6 @@ impl Interface { u: &[f64], reference: &[f64], committed: &Committed, - band: f64, tangent: bool, internal: &mut [f64], blocks: &mut Blocks, @@ -837,17 +835,21 @@ impl Interface { for v in &mut state.normal { *v = unit(*v).unwrap_or(*v); } - // Nodes closed, or open by less than `band` of their element's thickness: their penalty - // stays in the tangent. + // The closed nodes (their penalty in the tangent): those pressing, or those an active + // set held fixed closes. let mut engaged = vec![false; n]; for i in 0..self.nodes.len() { let (Some(Some(g)), Some(a), Some(e)) = (weighted.get(i), self.area.get(i), self.penalty.get(i)) else { continue }; let gap = g / a; - let lambda = multipliers.get(i).copied().unwrap_or(0.0); - if let (Some(slot), Some(p), Some(h)) = (state.gap.get_mut(i), state.pressure.get_mut(i), self.thickness.get(i)) { + let trial = multipliers.get(i).copied().unwrap_or(0.0) - e * gap; + let closed = match committed.active.as_ref().and_then(|s| s.get(i)) { + Some(held) => *held, + None => trial > 0.0, + }; + if let (Some(slot), Some(p)) = (state.gap.get_mut(i), state.pressure.get_mut(i)) { *slot = Some(gap); - *p = (lambda - e * gap).max(0.0); - engaged[i] = lambda - e * gap > -e * band * h; + *p = if closed { trial } else { 0.0 }; + engaged[i] = closed; } } // The pairs with a closed secondary corner, linearised. @@ -1172,6 +1174,29 @@ impl Interface { } } + /// The active set held at `committed`, checked at `state` (computed with it): `None` when it + /// is consistent there (no held-closed node pulls with more than ε_A times + /// [`GAP_TOLERANCE`] of its element's thickness, and no held-open node presses by more: a + /// node within that band touches with no pressure, and either status is the solution to the + /// tolerance), else the set with the offending nodes switched (the primal–dual active set + /// update; S. Hüeber, B. I. Wohlmuth, *A primal–dual active set strategy for non-linear + /// multibody contact problems*, CMAME 194:3147–3166, 2005). + pub fn update_held(&self, committed: &Committed, state: &State) -> Option> { + let held = committed.active.as_ref()?; + let mut next = held.clone(); + let mut changed = false; + for (i, closed) in next.iter_mut().enumerate() { + let (Some(Some(gap)), Some(e), Some(h)) = (state.gap.get(i), self.penalty.get(i), self.thickness.get(i)) else { continue }; + let trial = committed.multipliers.get(i).copied().unwrap_or(0.0) - e * gap; + let limit = e * GAP_TOLERANCE * h; + if (*closed && trial < -limit) || (!*closed && trial > limit) { + *closed = !*closed; + changed = true; + } + } + changed.then_some(next) + } + /// How far the augmentation is from converged at `state`, computed with `multipliers`: the /// largest gap of a closed node, or the gap equivalent λ_A / ε_A of the multiplier an open /// node keeps, over the thickness of the element under the node. @@ -1186,3 +1211,94 @@ impl Interface { worst } } + +#[cfg(test)] +mod tests { + use super::*; + use fieldcraft_model::{Elastic, Material, Section, SectionKind, Set}; + + /// Two separate quads facing each other, the upper one tilted and pressed into the lower + /// one, with general contact; `axisymmetric` off the axis. + fn pressed_quads(axisymmetric: bool) -> Model { + let mut m = Model::default(); + let x0 = if axisymmetric { 2.0 } else { 0.0 }; + let lower = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; + let upper = [[0.25, 1.03], [1.2, 1.01], [1.25, 1.9], [0.3, 1.95]]; + for (i, p) in lower.iter().chain(&upper).enumerate() { + m.set_node(i as u32 + 1, [x0 + p[0], p[1], 0.0]).unwrap(); + } + m.set_element(1, ElementKind::Quad4, vec![1, 2, 3, 4]).unwrap(); + m.set_element(2, ElementKind::Quad4, vec![5, 6, 7, 8]).unwrap(); + m.add_set("all", Set { elements: vec![1, 2], ..Set::default() }).unwrap(); + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: 10.0, poisson: 0.3 }), ..Material::default() }).unwrap(); + let kind = if axisymmetric { SectionKind::Axisymmetric {} } else { SectionKind::PlaneStrain { thickness: 1.5 } }; + m.add_section("s", Section { set: "all".into(), material: "Mat".into(), kind, formulation: Default::default() }).unwrap(); + m.add_contact("c", Contact::new(ContactSurfaces::General { surface: None })).unwrap(); + let step = fieldcraft_model::Step::new("Step", fieldcraft_model::Procedure::NonlinearStatic(Default::default())); + m.add_step(step, None).unwrap(); + m + } + + /// The contact's tangent blocks are the derivatives of its forces (central differences), in + /// plane strain and axisymmetric models. + #[test] + fn the_2d_contact_tangent_is_the_derivative_of_its_forces() { + for axisymmetric in [false, true] { + let m = pressed_quads(axisymmetric); + let resolved = m.resolve().unwrap(); + let mesh = Mesh::new(&m, &resolved, 1_000_000).unwrap(); + let interfaces = Interface::all(&m, &resolved, &mesh, true, &mut Vec::new()).unwrap(); + let n = mesh.dofs(); + // A small motion of every node, so no pair is symmetric by accident, and the upper + // quad pressed into the lower one. + let mut u: Vec = (0..n).map(|i| 0.01 * ((i * 7 % 5) as f64 - 2.0)).collect(); + for node in 4..8 { + if let Some(d) = translations(&mesh, node).ok().and_then(|d| d.get(1).copied()) { + u[d] -= 0.06; + } + } + let forces = |u: &[f64], blocks: &mut Blocks| { + let mut f = vec![0.0; n]; + for c in &interfaces { + let committed = c.committed(); + c.evaluate(&mesh, u, &vec![0.0; n], &committed, true, &mut f, blocks).unwrap(); + } + f + }; + let mut blocks = Blocks::new(); + forces(&u, &mut blocks); + // The dense tangent over (node, x/y) from the 3-per-node blocks. + let dof = |node: usize, k: usize| translations(&mesh, node).ok().and_then(|d| d.get(k).copied()); + let mut k = vec![0.0; n * n]; + for (nodes, data) in &blocks { + let w = 3 * nodes.len(); + for (a, na) in nodes.iter().enumerate() { + for (b, nb) in nodes.iter().enumerate() { + for r in 0..2 { + for c in 0..2 { + if let (Some(i), Some(j)) = (dof(*na, r), dof(*nb, c)) { + k[i * n + j] += data[(3 * a + r) * w + 3 * b + c]; + } + } + } + } + } + } + let scale = k.iter().fold(0.0f64, |a, v| a.max(v.abs())); + assert!(scale > 0.0, "the quads are in contact"); + let h = 1e-7; + for j in 0..n { + let shifted = |s: f64| { + let mut v = u.clone(); + v[j] += s; + forces(&v, &mut Blocks::new()) + }; + let (p, q) = (shifted(h), shifted(-h)); + for i in 0..n { + let fd = (p[i] - q[i]) / (2.0 * h); + assert!((fd - k[i * n + j]).abs() < 1e-5 * scale, "axisymmetric {axisymmetric}: ({i}, {j}) {fd} vs {}", k[i * n + j]); + } + } + } + } +} diff --git a/crates/solver/src/ldlt.rs b/crates/solver/src/ldlt.rs index 1c99f4b..53511fb 100644 --- a/crates/solver/src/ldlt.rs +++ b/crates/solver/src/ldlt.rs @@ -527,7 +527,9 @@ impl Symmetric { Kind::Definite => Factor::new(a, max_entries).map(Symmetric::Cholesky), Kind::Indefinite => LdltFactor::new(a, max_entries).map(Symmetric::Indefinite), Kind::Auto => { + let t0 = std::time::Instant::now(); let symbolic = Symbolic::new(a, max_entries)?; + if std::env::var("FC_TRACE").is_ok() { eprintln!("TRACE symbolic {:?}", t0.elapsed()); } match Factor::from_symbolic(&symbolic) { Ok(f) => Ok(Symmetric::Cholesky(f)), Err(FactorError::NotPositiveDefinite { .. }) => LdltFactor::from_symbolic(&symbolic).map(Symmetric::Indefinite), diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs index 91978bd..41aba74 100644 --- a/crates/solver/src/nonlinear.rs +++ b/crates/solver/src/nonlinear.rs @@ -217,6 +217,16 @@ pub struct History { pub iterations: usize, /// Arc-length steps: where the tangent's inertia changes along the path. pub critical_points: Vec, + /// Steps with stabilization: the energy it dissipated. + pub stabilization: Option, +} + +/// What a step's viscous stabilization did: the energy its damping dissipated, and the work of +/// the internal forces over the step (the energy the model stores, contact's included). +#[derive(Debug, Clone, Copy, PartialEq, Serialize)] +pub struct Stabilization { + pub dissipated: f64, + pub stored: f64, } impl History { @@ -229,6 +239,11 @@ impl History { /// Most displacement values the frames of one step keep (3 per node per frame). pub const MAX_FRAME_VALUES: usize = 30_000_000; + +/// A stabilization that dissipates more than this fraction of the work of the internal forces is +/// reported. +const STABILIZATION_WARNING: f64 = 0.05; + /// Line search: accept a step whose residual work is at most this fraction of the first. const LINE_SEARCH_ACCEPT: f64 = 0.8; /// Line search: most trial steps, and the range of the step. @@ -425,7 +440,6 @@ impl Problem<'_> { committed: &[Vec], contact: &[contact::Committed], reference: &[f64], - band: f64, ) -> Result { let layout = &self.mesh.layout; let n = u.len(); @@ -433,10 +447,13 @@ impl Problem<'_> { // Contact first: its face pairs extend the tangent's pattern. let mut blocks = contact::Blocks::new(); let mut states = Vec::with_capacity(self.contacts.len()); + let t0 = std::time::Instant::now(); for (i, c) in self.contacts.iter().enumerate() { let empty = contact::Committed::default(); - states.push(c.evaluate(self.mesh, u, reference, contact.get(i).unwrap_or(&empty), band, tangent, &mut internal, &mut blocks)?); + states.push(c.evaluate(self.mesh, u, reference, contact.get(i).unwrap_or(&empty), tangent, &mut internal, &mut blocks)?); } + TIMES[0].fetch_add(t0.elapsed().as_micros() as u64, std::sync::atomic::Ordering::Relaxed); + let t0 = std::time::Instant::now(); let mut k = tangent.then(|| { if blocks.is_empty() { self.pattern.clone() @@ -451,6 +468,8 @@ impl Problem<'_> { add_block(k, layout, nodes, kinds, &translation_block(data, nodes.len(), kinds))?; } } + TIMES[1].fetch_add(t0.elapsed().as_micros() as u64, std::sync::atomic::Ordering::Relaxed); + let t0 = std::time::Instant::now(); let mut history = Vec::with_capacity(self.mesh.elements.len()); for round in parallel::blocks(self.mesh.elements.len(), ROUND) { let start = round.start; @@ -473,6 +492,7 @@ impl Problem<'_> { } } } + TIMES[2].fetch_add(t0.elapsed().as_micros() as u64, std::sync::atomic::Ordering::Relaxed); if let Some(s) = &self.springs { let mut ku = vec![0.0; s.n]; s.mul(u.get(..s.n).unwrap_or_default(), &mut ku); @@ -585,14 +605,19 @@ struct Run<'a> { contact: Vec, /// The displacements of the last converged increment (all degrees of freedom). reference: Vec, - /// The band of the engaged contact nodes ([`contact::ENGAGED_BAND`] once an increment's - /// active set has chattered, for the rest of the step; else 0). - band: f64, augmentations: usize, /// The augmentation being solved in the current increment. augmentation: usize, /// The largest force of the converged states (contact's force scale, [`Run::residual`]). converged_scale: f64, + /// Viscous stabilization: the volume each free equation's node carries, and the damping + /// c V per free equation once the first increment has calibrated c; the current + /// increment's size as a fraction of the step (the pseudo-time step). + volumes: Option>, + damping: Option>, + pseudo_time: f64, + /// The energy the stabilization has dissipated over the converged increments. + dissipated: f64, } /// The outcome of a converged increment. @@ -624,6 +649,28 @@ impl Run<'_> { } } + /// Take the stabilization's damping forces (D/Δt)(q − q₀) out of the free residual `r` at + /// `q`, the increment starting from `q0`. + fn damp(&self, r: &mut [f64], q: &[f64], q0: &[f64]) { + if let Some(d) = &self.damping { + let c = 1.0 / self.pseudo_time; + for (((ri, di), a), b) in r.iter_mut().zip(d).zip(q).zip(q0) { + *ri -= c * di * (a - b); + } + } + } + + /// Add the stabilization's D/Δt to the free tangent. + fn damp_tangent(&self, kff: &mut Csr) -> std::result::Result<(), Failure> { + if let Some(d) = &self.damping { + let c = 1.0 / self.pseudo_time; + for (i, di) in d.iter().enumerate() { + add_entry(kff, i, i, c * di).map_err(Failure::Fatal)?; + } + } + Ok(()) + } + fn scale(&self) -> f64 { self.force_scale.max(self.force_floor.unwrap_or(0.0)).max(f64::MIN_POSITIVE) } @@ -636,7 +683,6 @@ impl Run<'_> { &self.committed, &self.contact, &self.reference, - self.band, )?) } @@ -647,7 +693,7 @@ impl Run<'_> { if let Some(u) = &self.undeformed { return Ok(u.clone()); } - let mut e = self.problem.evaluate(&vec![0.0; self.problem.mesh.dofs()], 0.0, true, &self.committed, &[], &[], 0.0)?; + let mut e = self.problem.evaluate(&vec![0.0; self.problem.mesh.dofs()], 0.0, true, &self.committed, &[], &[])?; let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; let [kff] = self.problem.red.reduce([&k])?; let u = factor_raw(&kff, self.limits).map(|_| ()); @@ -696,7 +742,7 @@ impl Run<'_> { /// factored (the increment then starts from `q0`). fn tangent_predictor(&mut self, q0: &[f64], from: f64, to: f64) -> std::result::Result, f64)>, Failure> { let red = self.problem.red; - let mut e = self.problem.evaluate(&red.expand_scaled(q0, from), to, true, &self.committed, &self.contact, &self.reference, self.band)?; + let mut e = self.problem.evaluate(&red.expand_scaled(q0, from), to, true, &self.committed, &self.contact, &self.reference)?; let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; let rf: Vec = free_residual(red, &e).iter().zip(prescribed_force(red, &k, to - from)).map(|(a, b)| a - b).collect(); let mut dq = rf.clone(); @@ -721,6 +767,12 @@ impl Run<'_> { ) -> std::result::Result { let red = self.problem.red; let tol = self.settings.tolerances; + // A held active set is the increment's own (`Committed::active`). + for c in &mut self.contact { + c.active = None; + } + let mut held = false; + let mut updates = 0usize; // Start from the extrapolated predictor when it is a valid state, else from q0. let start = predictor.map(|p| self.state(p, lambda, true).map(|e| (p.to_vec(), e))); let (mut q, mut eval) = match start { @@ -729,6 +781,7 @@ impl Run<'_> { Some(Err(Failure::Cutback(_))) | None => (q0.to_vec(), self.state(q0, lambda, true)?), }; let mut r = self.residual(&eval); + self.damp(&mut r, &q, q0); // With path-dependent materials, an increment that starts within the force tolerance (a // predictor that already solved it) has only round-off left for the energy criterion to // measure against: it is waived. @@ -743,9 +796,18 @@ impl Run<'_> { let mut sets = (None, active_sets(&eval)); let mut flips = 0usize; let max_iterations = usize::try_from(self.settings.max_iterations).unwrap_or(usize::MAX); - for _ in 0..max_iterations { + // Holding an active set starts the increment's iterations again (once). + let mut remaining = max_iterations; + while remaining > 0 { + remaining -= 1; // The round-off rule: nothing left to correct. if norm_inf(&r) <= ROUND_OFF * tol.force * self.scale() && self.force_floor.is_some() { + if self.release(&mut eval, &q, lambda, &mut updates)? { + r = self.residual(&eval); + self.damp(&mut r, &q, q0); + remaining = max_iterations; + continue; + } return Ok(Converged { q, evaluation: eval, energy: first_energy }); } if self.iterations >= self.limits.max_newton_iterations { @@ -756,14 +818,17 @@ impl Run<'_> { } self.iterations += 1; let k = eval.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; - let [kff] = red.reduce([&k])?; + let [mut kff] = red.reduce([&k])?; // The first tangent of the step: a singular one is a modelling error, not a cutback. let first = self.force_floor.is_none(); if first { self.seed_force_floor(&kff); } + self.damp_tangent(&mut kff)?; drop(k); + let t0 = std::time::Instant::now(); let factor = self.factor(&kff, first)?; + TIMES[3].fetch_add(t0.elapsed().as_micros() as u64, std::sync::atomic::Ordering::Relaxed); self.factor_entries = factor.entries(); drop(kff); let mut dq = r.clone(); @@ -772,24 +837,38 @@ impl Run<'_> { return Err(Failure::Cutback("the Newton correction is not finite".into())); } let g0 = dot(&dq, &r); - let (s, next) = self.line_search(&q, &dq, lambda, g0)?; + let (s, next) = self.line_search(&q, q0, &dq, lambda, g0)?; + let negative = factor.inertia().negative; for (x, d) in q.iter_mut().zip(&dq) { *x += s * d; } eval = next; - if !self.problem.contacts.is_empty() { + let holding = self.contact.iter().any(|c| c.active.is_some()); + if !self.problem.contacts.is_empty() && !holding { let set = active_sets(&eval); if sets.0 == Some(set) && set != sets.1 { flips += 1; if flips >= contact::CHATTER_LIMIT { - // From now on the nodes about to touch stay in the tangent. - self.band = contact::ENGAGED_BAND; - return Err(Failure::Cutback("the contact status alternates between two states (chattering)".into())); + if held { + return Err(Failure::Cutback("the contact status alternates between two states (chattering)".into())); + } + // Hold the current active set until the increment converges (then + // check it, `release`): the nodes touching with no pressure stop + // flipping, and the tangent stays the derivative of the forces. + held = true; + for (c, st) in self.contact.iter_mut().zip(&eval.contact) { + c.active = Some(st.pressure.iter().map(|p| *p > 0.0).collect()); + } + eval = self.state(&q, lambda, true)?; + flips = 0; + remaining = max_iterations; } } - sets = (Some(sets.1), set); + sets = (Some(sets.1), active_sets(&eval)); } - let r_old = std::mem::replace(&mut r, self.residual(&eval)); + let mut r_new = self.residual(&eval); + self.damp(&mut r_new, &q, q0); + let r_old = std::mem::replace(&mut r, r_new); let energy = (s * dot(&dq, &r_old)).abs(); let first = *first_energy.get_or_insert(energy); let moved = q.iter().zip(q0).fold(0.0f64, |m, (a, b)| m.max((a - b).abs())); @@ -799,7 +878,7 @@ impl Run<'_> { displacement: s * norm_inf(&dq) / reference, energy: if first > 0.0 { energy / first } else { 0.0 }, line_search: s, - negative_pivots: factor.inertia().negative, + negative_pivots: negative, augmentation: self.augmentation, }; iterations.push(it); @@ -808,6 +887,13 @@ impl Run<'_> { } let energy_met = it.energy <= tol.energy || started_converged; if (it.force <= tol.force && it.displacement <= tol.displacement && energy_met) || it.force <= ROUND_OFF * tol.force { + if self.release(&mut eval, &q, lambda, &mut updates)? { + r = self.residual(&eval); + self.damp(&mut r, &q, q0); + history.clear(); + remaining = max_iterations; + continue; + } return Ok(Converged { q, evaluation: eval, energy: first_energy }); } history.push(it.force); @@ -861,7 +947,8 @@ impl Run<'_> { self.augmentations += 1; self.augmentation = *count; for (l, s) in self.contact.iter_mut().zip(&c.evaluation.contact) { - l.multipliers.clone_from(&s.pressure); + // A held-closed node may have pulled (within the tolerance); multipliers press. + l.multipliers = s.pressure.iter().map(|p| p.max(0.0)).collect(); } // From the converged state, measured against the increment's start and its first // correction like the increment itself. @@ -872,7 +959,37 @@ impl Run<'_> { /// The line search along `dq` from `q`: the accepted step and the state there (with the /// tangent). - fn line_search(&mut self, q: &[f64], dq: &[f64], lambda: f64, g0: f64) -> std::result::Result<(f64, Evaluation), Failure> { + /// At a converged state `eval` (at `q`, load factor `lambda`) with active sets held: when + /// every one is consistent there ([`Interface::update_held`]) they are released and the + /// state is the solution (false); otherwise the offending nodes switch, at most + /// [`contact::MAX_ACTIVE_SET_UPDATES`] times, `eval` is evaluated again and the iterations + /// go on (true). + fn release(&mut self, eval: &mut Evaluation, q: &[f64], lambda: f64, updates: &mut usize) -> std::result::Result { + if self.contact.iter().all(|c| c.active.is_none()) { + return Ok(false); + } + let next: Vec>> = + self.problem.contacts.iter().zip(&self.contact).zip(&eval.contact).map(|((i, c), s)| i.update_held(c, s)).collect(); + if next.iter().all(Option::is_none) { + for c in &mut self.contact { + c.active = None; + } + return Ok(false); + } + *updates += 1; + if *updates > contact::MAX_ACTIVE_SET_UPDATES { + return Err(Failure::Cutback("the contact's active set does not settle".into())); + } + for (c, n) in self.contact.iter_mut().zip(next) { + if n.is_some() { + c.active = n; + } + } + *eval = self.state(q, lambda, true)?; + Ok(true) + } + + fn line_search(&mut self, q: &[f64], q0: &[f64], dq: &[f64], lambda: f64, g0: f64) -> std::result::Result<(f64, Evaluation), Failure> { let at = |s: f64| -> Vec { q.iter().zip(dq).map(|(a, b)| a + s * b).collect() }; // The full step, with its tangent (taken in most iterations). let full = self.state(&at(1.0), lambda, true); @@ -880,11 +997,15 @@ impl Run<'_> { return full.map(|e| (1.0, e)); } let red = self.problem.red; - let work = |e: &Evaluation| dot(dq, &free_residual(red, e)); + let work = |e: &Evaluation, s: f64| { + let mut r = free_residual(red, e); + self.damp(&mut r, &at(s), q0); + dot(dq, &r) + }; let mut s = 1.0; let mut g = match &full { Ok(e) => { - let g1 = work(e); + let g1 = work(e, 1.0); if g1.abs() <= LINE_SEARCH_ACCEPT * g0.abs() { return full.map(|e| (1.0, e)); } @@ -906,7 +1027,7 @@ impl Run<'_> { } s = next; g = match self.state(&at(s), lambda, false) { - Ok(e) => Some(work(&e)), + Ok(e) => Some(work(&e, s)), Err(Failure::Cutback(_)) => None, Err(f) => return Err(f), }; @@ -929,6 +1050,26 @@ impl Run<'_> { } } +/// Per free equation, the volume its node carries (the lumped mass of a unit density): the +/// stabilization's mass matrix. +fn nodal_volumes(problem: &Problem) -> Result> { + let mut v = vec![0.0; problem.mesh.dofs()]; + for (e, dofs) in problem.mesh.elements.iter().zip(&problem.element_dofs) { + let lumped = e.kernel()?.lumped_mass(1.0).map_err(|err| Error::element(e.label, err))?; + let per = dofs.len() / lumped.len().max(1); + for (n, m) in lumped.iter().enumerate() { + for d in dofs.iter().skip(n * per).take(per) { + if let Some(x) = v.get_mut(*d) { + *x += m; + } + } + } + } + Ok(problem.red.restrict(&v)) +} + +pub static TIMES: [std::sync::atomic::AtomicU64; 4] = [const { std::sync::atomic::AtomicU64::new(0) }; 4]; + /// Refuse what nonlinear steps cannot take yet, and warn about what they simplify. fn check_supported(model: &Model, mesh: &Mesh, step: &Step, settings: &NonlinearStatic, warnings: &mut Vec) -> Result<()> { let at = |what: String| Error::NotSupported(format!("{what} in nonlinear steps (step `{}`)", step.name)); @@ -1178,7 +1319,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R .collect::>>()?; // The undeformed state first: an element that is bad before anything moves is an error in // the model, not a reason to cut increments back. - problem.evaluate(&vec![0.0; mesh.dofs()], 0.0, false, &committed, &[], &[], 0.0)?; + problem.evaluate(&vec![0.0; mesh.dofs()], 0.0, false, &committed, &[], &[])?; let path_dependent = problem.materials.iter().any(MaterialModel::is_path_dependent); let mut run = Run { problem: &problem, @@ -1194,10 +1335,13 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R undeformed: None, contact: problem.contacts.iter().map(Interface::committed).collect(), reference: vec![0.0; mesh.dofs()], - band: 0.0, augmentations: 0, augmentation: 0, converged_scale: 0.0, + volumes: settings.stabilization.map(|_| nodal_volumes(&problem)).transpose()?, + damping: None, + pseudo_time: 1.0, + dissipated: 0.0, }; let mut history = History::default(); @@ -1228,6 +1372,9 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R let mut size = settings.initial_increment; let mut fast = 0usize; let mut last: Option = None; + // The work of the internal forces over the step (the energy the model stores, with + // contact's), the stabilization's measure. + let mut stored = 0.0f64; // The last converged increment's free motion and the change of the load factor over it. let mut previous: Option<(Vec, f64)> = None; let max_increments = usize::try_from(settings.max_increments).unwrap_or(usize::MAX); @@ -1262,6 +1409,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R } } let saved = run.contact.clone(); + run.pseudo_time = this; let mut augmentations = 0; run.augmentation = 0; let solved = run @@ -1270,8 +1418,36 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R if solved.is_err() { run.contact = saved; } + if std::env::var("FC_TRACE").is_ok() { + let f: Vec = iterations.iter().map(|i| format!("{:.0e}{}", i.force, if i.line_search < 1.0 { format!("/s{:.2}", i.line_search) } else { String::new() })).collect(); + let closed: usize = match &solved { Ok(c) => c.evaluation.contact.iter().map(|s| s.pressure.iter().filter(|p| **p > 0.0).count()).sum(), Err(_) => 0 }; + eprintln!("TRACE times contact {} pattern {} elements {} factor {} ms", TIMES[0].load(std::sync::atomic::Ordering::Relaxed) / 1000, TIMES[1].load(std::sync::atomic::Ordering::Relaxed) / 1000, TIMES[2].load(std::sync::atomic::Ordering::Relaxed) / 1000, TIMES[3].load(std::sync::atomic::Ordering::Relaxed) / 1000); + eprintln!("TRACE damp {:?} dissipated {:.3e} stored {stored:.3e}", run.damping.as_ref().map(|d| d.iter().fold(0.0f64, |m, v| m.max(*v))), run.dissipated); + eprintln!("TRACE inc {number} lam {target_lambda:.5} size {:.2e} closed {closed} neg {:?} {} [{}] {}", change.abs(), iterations.last().map(|i| i.negative_pivots), if solved.is_ok() { "ok" } else { "CUT" }, f.join(" "), match &solved { Err(Failure::Cutback(r)) => r.clone(), _ => String::new() }); + } match solved { Ok(c) => { + if let Some(d) = &run.damping { + run.dissipated += c.q.iter().zip(&q).zip(d).map(|((a, b), di)| di * (a - b) * (a - b)).sum::() / this; + } + // The work of the internal forces over the increment (trapezoidal; zero forces + // at the unloaded start). + let (u0, u1) = (red.expand_scaled(&q, lambda), red.expand_scaled(&c.q, target_lambda)); + let zero = vec![0.0; u0.len()]; + let before = last.as_ref().map_or(zero.as_slice(), |e| e.internal.as_slice()); + let work = 0.5 * before.iter().zip(&c.evaluation.internal).zip(u0.iter().zip(&u1)).map(|((f0, f1), (a, b))| (f0 + f1) * (b - a)).sum::(); + stored += work; + // The stabilization's damping, calibrated on the first increment (solved + // without it): moving on at that increment's rate Δq/Δt, c V would dissipate + // the fraction α of the work the step does at that rate, W₁/Δt² (the work + // grows with the square of the motion), so c = α W₁ / (Δqᵀ V Δq). + if let (Some(alpha), Some(v), None) = (settings.stabilization, &run.volumes, &run.damping) { + let motion: f64 = c.q.iter().zip(&q).zip(v).map(|((a, b), vi)| vi * (a - b) * (a - b)).sum(); + if motion > 0.0 && work.abs() > 0.0 { + let factor = alpha * work.abs() / motion; + run.damping = Some(v.iter().map(|vi| factor * vi).collect()); + } + } previous = Some((c.q.iter().zip(&q).map(|(a, b)| a - b).collect(), change)); // The increment's own Newton iterations size the next one (augmentations // solve the same increment again). @@ -1327,13 +1503,23 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R } } } + if run.damping.is_some() { + history.stabilization = Some(Stabilization { dissipated: run.dissipated, stored }); + if run.dissipated > STABILIZATION_WARNING * stored.abs() { + warnings.push(format!( + "the stabilization dissipated {:.3e}, {:.1} % of the work of the internal forces ({stored:.3e}): lower it, or the result is not quasi-static", + run.dissipated, + 100.0 * run.dissipated / stored.abs().max(f64::MIN_POSITIVE) + )); + } + } (q, lambda, last) }; history.iterations = run.iterations; let u = red.expand_scaled(&q, lambda); let evaluation = match last { Some(e) => e, - None => problem.evaluate(&u, lambda, false, &run.committed, &run.contact, &run.reference, 0.0)?, + None => problem.evaluate(&u, lambda, false, &run.committed, &run.contact, &run.reference)?, }; if let Some((at, negative)) = history .increments diff --git a/crates/solver/src/ordering.rs b/crates/solver/src/ordering.rs index 4e04340..4c0055b 100644 --- a/crates/solver/src/ordering.rs +++ b/crates/solver/src/ordering.rs @@ -86,6 +86,182 @@ pub fn nested_dissection(xyz: &[[f64; 3]], adjacency: &[Vec]) -> Vec]) -> Vec { + let n = adjacency.len(); + let mut order = Vec::with_capacity(n); + // The part each node belongs to (a token per split), and its level in the current search. + let mut owner = vec![0usize; n]; + let mut level = vec![usize::MAX; n]; + let mut token = 1usize; + enum Task { + Split(Vec), + Emit(Vec), + } + let mut stack = vec![Task::Split((0..n).collect())]; + // The level structure of the nodes of part `token` reached from `root`. + let levels = |root: usize, token: usize, owner: &[usize], level: &mut Vec| -> Vec> { + let mut out: Vec> = vec![vec![root]]; + if let Some(l) = level.get_mut(root) { + *l = 0; + } + loop { + let depth = out.len(); + let mut next = Vec::new(); + for &v in out.last().into_iter().flatten() { + for &w in adjacency.get(v).into_iter().flatten() { + if owner.get(w) == Some(&token) && level.get(w) == Some(&usize::MAX) { + if let Some(l) = level.get_mut(w) { + *l = depth; + } + next.push(w); + } + } + } + if next.is_empty() { + break; + } + out.push(next); + } + out + }; + let degree = |v: usize, token: usize, owner: &[usize]| adjacency.get(v).map_or(0, |row| row.iter().filter(|w| owner.get(**w) == Some(&token)).count()); + while let Some(task) = stack.pop() { + let part = match task { + Task::Emit(sep) => { + order.extend(sep); + continue; + } + Task::Split(p) => p, + }; + if part.len() <= LEAF { + order.extend(part); + continue; + } + token += 1; + for &v in &part { + if let Some(o) = owner.get_mut(v) { + *o = token; + } + } + let reset = |level: &mut Vec, nodes: &[usize]| { + for &v in nodes { + if let Some(l) = level.get_mut(v) { + *l = usize::MAX; + } + } + }; + // A pseudo-peripheral root: from the first node, the end of the deepest level structure + // with the fewest neighbours, while the depth grows. + let mut root = part.first().copied().unwrap_or(0); + let mut structure = levels(root, token, &owner, &mut level); + for _ in 0..PERIPHERAL_SEARCHES { + let reached: Vec = structure.iter().flatten().copied().collect(); + reset(&mut level, &reached); + let Some(candidate) = structure.last().and_then(|last| last.iter().copied().min_by_key(|v| (degree(*v, token, &owner), *v))) else { break }; + let deeper = levels(candidate, token, &owner, &mut level); + if deeper.len() <= structure.len() { + let reached: Vec = deeper.iter().flatten().copied().collect(); + reset(&mut level, &reached); + structure = levels(root, token, &owner, &mut level); + break; + } + (root, structure) = (candidate, deeper); + } + let reached: Vec = structure.iter().flatten().copied().collect(); + reset(&mut level, &reached); + if reached.len() < part.len() { + // Disconnected: this component, then the rest. + let mut inside = vec![false; 0]; + inside.resize(n, false); + for &v in &reached { + if let Some(x) = inside.get_mut(v) { + *x = true; + } + } + let rest: Vec = part.iter().copied().filter(|v| !inside.get(*v).copied().unwrap_or(false)).collect(); + let mut component = reached; + component.sort_unstable(); + stack.push(Task::Split(rest)); + stack.push(Task::Split(component)); + continue; + } + if structure.len() < 3 { + order.extend(part); + continue; + } + // The middle level: the first whose levels up to it hold half the part. + let mut count = 0; + let mut middle = 1; + for (i, l) in structure.iter().enumerate().skip(1) { + count += l.len(); + if 2 * count >= part.len() { + middle = i.min(structure.len() - 2); + break; + } + } + // Its nodes that touch the next level separate the levels before it from those after. + for (i, l) in structure.iter().enumerate() { + for &v in l { + if let Some(x) = level.get_mut(v) { + *x = i; + } + } + } + let touches_next = |v: usize, level: &[usize]| adjacency.get(v).is_some_and(|row| row.iter().any(|w| owner.get(*w) == Some(&token) && level.get(*w) == Some(&(middle + 1)))); + let mut separator = Vec::new(); + let mut before: Vec = structure.iter().take(middle).flatten().copied().collect(); + for &v in structure.get(middle).into_iter().flatten() { + if touches_next(v, &level) { separator.push(v) } else { before.push(v) } + } + let after: Vec = structure.iter().skip(middle + 1).flatten().copied().collect(); + reset(&mut level, &part); + before.sort_unstable(); + separator.sort_unstable(); + let mut after = after; + after.sort_unstable(); + stack.push(Task::Emit(separator)); + stack.push(Task::Split(after)); + stack.push(Task::Split(before)); + } + order +} + +/// Searches for a pseudo-peripheral node. +const PERIPHERAL_SEARCHES: usize = 8; + +/// The entries of the Cholesky factor of the node graph `adjacency` eliminated in `order` (one +/// per node pair): the fill an ordering gives. +pub fn fill(adjacency: &[Vec], order: &[usize]) -> usize { + let n = adjacency.len(); + let mut position = vec![0usize; n]; + for (k, &v) in order.iter().enumerate() { + if let Some(p) = position.get_mut(v) { + *p = k; + } + } + let mut ptr = vec![0usize]; + let mut idx = Vec::new(); + for &v in order { + let k = position.get(v).copied().unwrap_or(0); + let mut rows: Vec = adjacency.get(v).into_iter().flatten().filter_map(|w| position.get(*w).copied()).filter(|r| *r <= k).collect(); + rows.push(k); + rows.sort_unstable(); + rows.dedup(); + idx.extend(rows); + ptr.push(idx.len()); + } + let val = vec![1.0; idx.len()]; + crate::cholesky::factor_entries(&crate::cholesky::UpperColumns { n, ptr, idx, val }) +} + #[cfg(test)] mod tests { use super::*; diff --git a/crates/solver/tests/contact.rs b/crates/solver/tests/contact.rs index 9b8cb4c..cedc30c 100644 --- a/crates/solver/tests/contact.rs +++ b/crates/solver/tests/contact.rs @@ -686,14 +686,19 @@ fn a_rubber_bellows_folds_onto_itself() { }, ) .unwrap(); - let settings = NonlinearStatic { initial_increment: 0.05, max_increment: 0.1, max_increments: 400, ..NonlinearStatic::default() }; + let stabilization = std::env::var("STAB").ok().and_then(|v| v.parse().ok()); + let settings = NonlinearStatic { initial_increment: 0.05, max_increment: 0.1, max_increments: 1000, line_search: stabilization.is_some(), stabilization, ..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(); + let squeeze: f64 = std::env::var("SQUEEZE").ok().and_then(|v| v.parse().ok()).unwrap_or(0.8); + b.bc("squeeze", [None, None, Some(-squeeze)], move |x| near(x[2], height)); + let t0 = std::time::Instant::now(); + let r = solve_nonlinear(&b.model, Some("Step"), &Limits::default()); + eprintln!("bellows squeeze {squeeze}: {:?} {}", t0.elapsed(), r.as_ref().err().map(ToString::to_string).unwrap_or_default()); + let r = r.unwrap(); report(&r); let c = &r.contacts[0]; eprintln!("closed {} penetration {:.2e} augmentations {}", c.closed(), c.penetration, c.augmentations); @@ -806,3 +811,11 @@ fn what_2d_contact_refuses() { let e = solve_nonlinear(&b.model, Some("Step"), &Limits::default()).unwrap_err().to_string(); assert!(e.contains("arc-length"), "{e}"); } + +#[test] +#[ignore = "scratch"] +fn scratch_quadratic_2d() { + let (_, r2) = planar_c_ring(0.15, 1); + report(&r2); + quadratic(&r2); +}