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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

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 {

crates/decks/src/lsdyna/mod.rs

@@ -290,6 +290,9 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result<Plan<'a>> {
             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)",

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();

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");

crates/engine/src/cmd/define.rs

@@ -292,7 +292,7 @@ fn fluid_boundary(s: &mut Session, p: &Value) -> Result<Value> {
 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,

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] =

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,

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<u32> {
+    let reversed: Vec<u32> = 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<Value> {
+    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<u32>)> = 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<u32> = {
+        let mut v: Vec<u32> = 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<u32, usize> = 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}");
+        }
+    }
+}

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::<usize>().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::<String>());
+        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::<String>(), Err(e) => e.to_string() });
+        if r.is_ok() {
+            let labels: Vec<u32> = 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<Value> = 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();
+        }
+    }
 }

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<Vec<Piece>> {
+    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<Value>) {
+    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<Value>; 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::<Vec<_>>()}}));
+    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<Value>, 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)
+}

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<ArcLength>,
+    /// 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<f64>,
 }
 
 /// 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;
 

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<usize>) {
     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<Vec<usize>> = 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<Vec<usize>> = 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)
 }
 

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<usize>, Vec<usize>) {
     let parent = etree(a);
     let mut counts = vec![1usize; a.n];

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<f64>,
     pub traction: Vec<Vec3>,
     pub slipping: Vec<bool>,
+    /// 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<Vec<bool>>,
 }
 
 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<Vec<bool>> {
+        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<f64> = (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]);
+                }
+            }
+        }
+    }
+}

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),

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<CriticalPoint>,
+    /// Steps with stabilization: the energy it dissipated.
+    pub stabilization: Option<Stabilization>,
+}
+
+/// 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<f64>],
         contact: &[contact::Committed],
         reference: &[f64],
-        band: f64,
     ) -> Result<Evaluation> {
         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<contact::Committed>,
     /// The displacements of the last converged increment (all degrees of freedom).
     reference: Vec<f64>,
-    /// 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<Vec<f64>>,
+    damping: Option<Vec<f64>>,
+    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<Option<(Vec<f64>, 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<f64> = 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<Converged, Failure> {
         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<bool, Failure> {
+        if self.contact.iter().all(|c| c.active.is_none()) {
+            return Ok(false);
+        }
+        let next: Vec<Option<Vec<bool>>> =
+            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<f64> { 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<Vec<f64>> {
+    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<String>) -> 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<String>) -> R
         .collect::<Result<Vec<_>>>()?;
     // 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<String>) -> 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<String>) -> R
         let mut size = settings.initial_increment;
         let mut fast = 0usize;
         let mut last: Option<Evaluation> = 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>, 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<String>) -> 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<String>) -> R
             if solved.is_err() {
                 run.contact = saved;
             }
+            if std::env::var("FC_TRACE").is_ok() {
+                let f: Vec<String> = 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::<f64>() / 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::<f64>();
+                    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<String>) -> 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

crates/solver/src/ordering.rs

@@ -86,6 +86,182 @@ pub fn nested_dissection(xyz: &[[f64; 3]], adjacency: &[Vec<usize>]) -> Vec<usiz
     order
 }
 
+/// A permutation of the nodes by nested dissection on the graph alone: each part's separator is
+/// a middle level of the breadth-first level structure rooted at a pseudo-peripheral node, so
+/// separators follow the mesh rather than a coordinate plane (A. George, J. W. H. Liu, *An
+/// automatic nested dissection algorithm for irregular finite element problems*, SIAM J. Numer.
+/// Anal. 15(5):1053–1069, 1978; the pseudo-peripheral node of N. E. Gibbs, W. G. Poole,
+/// P. K. Stockmeyer, SIAM J. Numer. Anal. 13(2):236–250, 1976). A thin curved wall (a shell, a
+/// bellows) is cut across its surface instead of many times by a plane. Disconnected parts are
+/// ordered one after the other; ties go by node index, so the ordering is deterministic.
+pub fn graph_dissection(adjacency: &[Vec<usize>]) -> Vec<usize> {
+    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<usize>),
+        Emit(Vec<usize>),
+    }
+    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<usize>| -> Vec<Vec<usize>> {
+        let mut out: Vec<Vec<usize>> = 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<usize>, 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<usize> = 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<usize> = deeper.iter().flatten().copied().collect();
+                reset(&mut level, &reached);
+                structure = levels(root, token, &owner, &mut level);
+                break;
+            }
+            (root, structure) = (candidate, deeper);
+        }
+        let reached: Vec<usize> = 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<usize> = 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<usize> = 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<usize> = 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<usize>], 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<usize> = 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::*;

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);
+}