diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 2c19212..e48fafb 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -1202,8 +1202,27 @@ impl Deck { } let nlgeom = if n.nonlinear_geometry { "YES" } else { "NO" }; self.line(&format!("*STEP, NAME={}, NLGEOM={nlgeom}, INC={}", step.name, n.max_increments)); - self.line("*STATIC"); - self.numbers_line("", &[n.initial_increment, 1.0, n.min_increment, n.max_increment])?; + if let Some(arc) = &n.arc_length { + // Riks: initial arc, period, minimum and maximum arcs, the load + // proportionality factor and the displacement that end the step. + if arc.constraint != fieldcraft_model::ArcConstraint::Cylindrical { + return Err(Error::Unsupported { + format: FORMAT, + what: format!("step `{}`: a spherical arc-length constraint", step.name), + }); + } + self.line("*STATIC, RIKS"); + self.numbers_line("", &[n.initial_increment, 1.0, n.min_increment, n.max_increment, arc.max_load_factor])?; + if let Some(l) = &arc.displacement_limit { + self.out.pop(); + let _ = write!(self.out, ", {}, {}, ", l.node, l.direction); + push_number(&mut self.out, l.displacement)?; + self.out.push('\n'); + } + } else { + self.line("*STATIC"); + self.numbers_line("", &[n.initial_increment, 1.0, n.min_increment, n.max_increment])?; + } if n.line_search { self.line("*CONTROLS, PARAMETERS=LINE SEARCH"); self.line("5,"); diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 8d86150..c2516f7 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -874,6 +874,9 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste d.card([Field::from(1)]); } Procedure::NonlinearStatic(n) => { + if n.arc_length.is_some() { + return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: arc-length control", step.name) }); + } if n.imperfection.is_some() { return Err(Error::Unsupported { format: FORMAT, diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs index 6da5b7b..74ebf37 100644 --- a/crates/decks/tests/abaqus.rs +++ b/crates/decks/tests/abaqus.rs @@ -845,3 +845,25 @@ fn buckling_steps() { assert!(abaqus::write(&m).unwrap_err().to_string().contains("after a general step")); assert!(fieldcraft_decks::lsdyna::write(&m).is_err() && fieldcraft_decks::nastran::write(&m).is_err()); } + +/// An arc-length step is `*STATIC, RIKS` with the arcs, the maximum load factor and the +/// displacement that end it; a spherical constraint is refused, and LS-DYNA refuses the step. +#[test] +fn arc_length_steps() { + use fieldcraft_model::{ArcConstraint, ArcLength, DisplacementLimit, NonlinearStatic}; + let mut m = cantilever_mesh(); + m.add_load("end", Load::Force { target: Target::Set("tip".into()), vector: [0.0, 0.0, -100.0], distribution: Distribution::Total }).unwrap(); + let arc = ArcLength { + max_load_factor: 3.0, + displacement_limit: Some(DisplacementLimit { node: 2, direction: 3, displacement: -5.0 }), + ..ArcLength::default() + }; + let settings = NonlinearStatic { arc_length: Some(arc.clone()), ..NonlinearStatic::default() }; + m.add_step(step("Snap", Procedure::NonlinearStatic(settings), &["clamp"], &["end"]), None).unwrap(); + let deck = abaqus::write(&m).unwrap(); + assert!(deck.contains("*STATIC, RIKS\n0.1, 1.0, 1.0E-5, 1.0, 3.0, 2, 3, -5.0\n"), "{deck}"); + assert!(fieldcraft_decks::lsdyna::write(&m).is_err()); + let spherical = NonlinearStatic { arc_length: Some(ArcLength { constraint: ArcConstraint::Spherical, ..arc }), ..NonlinearStatic::default() }; + m.replace_step(step("Snap", Procedure::NonlinearStatic(spherical), &["clamp"], &["end"])).unwrap(); + assert!(abaqus::write(&m).is_err()); +} diff --git a/crates/engine/src/cmd/forms.rs b/crates/engine/src/cmd/forms.rs index c8301e7..65d0abd 100644 --- a/crates/engine/src/cmd/forms.rs +++ b/crates/engine/src/cmd/forms.rs @@ -448,6 +448,15 @@ const PROCEDURE: Field = variants( num("scale", "Largest component", D::Length).req(), ], ), + group( + "arc_length", + "Arc length", + &[ + choice("constraint", "Constraint", &[("cylindrical", "Cylindrical"), ("spherical", "Spherical")]).def("\"cylindrical\""), + positive("max_load_factor", "Max load factor", D::None).def("1"), + json("displacement_limit", "Displacement limit").hint("{\"node\": 1, \"direction\": 3, \"displacement\": -10}"), + ], + ), ], ), variant("buckle", "Buckling", &[int("modes", "Modes", 1, 100_000).req().def("4"), num("shift", "Shift", D::None)]), diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 2b424d2..929dad0 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -237,6 +237,8 @@ fn run_nonlinear(s: &mut Session, p: &Value) -> Result { "increments": history.map(|h| &h.increments), "iterations": history.map(|h| h.iterations), "frames": history.map(|h| h.load_factors()), + "load_factor": history.and_then(|h| h.frames.last()).map(|f| f.load_factor), + "critical_points": history.map(|h| &h.critical_points), "max_displacement": largest(r.nodes.iter().zip(&r.displacement).map(|(l, u)| (*l, crate::results::magnitude(*u)))), "max_von_mises": largest(r.nodes.iter().zip(&r.von_mises).map(|(l, v)| (*l, *v))), "max_plastic_strain": largest(r.elements.iter().flat_map(|e| e.plastic_strain.iter().map(|v| (e.label, *v)))), diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index 361e868..f5ba1de 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -67,8 +67,9 @@ pub use mesh::{ pub use sections::{BeamProfile, Extent, Formulation, Integration, Section, SectionKind}; pub use sets::Set; pub use steps::{ - BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, Imperfection, InitialVelocity, MAX_LOAD_PATH_POINTS, MAX_NONLINEAR_INCREMENTS, - MAX_NONLINEAR_ITERATIONS, NonlinearStatic, OutputVariable, Procedure, RigidWall, Step, StressRate, Tolerances, Transient, + ArcConstraint, ArcLength, BulkViscosity, DisplacementLimit, ExplicitDynamic, Hourglass, HourglassKind, Imperfection, InitialVelocity, + MAX_LOAD_PATH_POINTS, MAX_NONLINEAR_INCREMENTS, MAX_NONLINEAR_ITERATIONS, NonlinearStatic, OutputVariable, Procedure, RigidWall, Step, + StressRate, Tolerances, Transient, }; pub use units::UnitSystem; diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index 2fb30bd..10bd7ce 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -149,6 +149,51 @@ pub struct NonlinearStatic { /// coordinates before this step is solved (the step and its results are on that geometry). #[serde(skip_serializing_if = "Option::is_none")] pub imperfection: Option, + /// Arc-length (Riks–Crisfield) control: the load factor becomes an unknown, so the step can + /// follow its equilibrium path through limit points (snap-through, snap-back). The first + /// increment's load factor change is `initial_increment`; the increment sizes become arc + /// lengths in proportion. `None`: load control. + #[serde(skip_serializing_if = "Option::is_none")] + pub arc_length: Option, +} + +/// The constraint of an arc-length step. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum ArcConstraint { + /// ‖Δu‖ = Δl (Crisfield's cylindrical arc). + #[default] + Cylindrical, + /// ‖Δu‖² + Δλ² ‖u_f‖² = Δl², u_f the first increment's tangent displacement of the unit + /// load (the spherical arc, scaled to displacement units). + Spherical, +} + +/// When an arc-length step ends: at `max_load_factor` (default 1), when the displacement of +/// `node` along `direction` (1, 2, 3) reaches `displacement` in magnitude, or at the step's +/// increment cap, whichever comes first. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(default, deny_unknown_fields)] +pub struct ArcLength { + pub constraint: ArcConstraint, + pub max_load_factor: f64, + #[serde(skip_serializing_if = "Option::is_none")] + pub displacement_limit: Option, +} + +impl Default for ArcLength { + fn default() -> ArcLength { + ArcLength { constraint: ArcConstraint::Cylindrical, max_load_factor: 1.0, displacement_limit: None } + } +} + +/// A displacement that ends an arc-length step. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct DisplacementLimit { + pub node: u32, + pub direction: u8, + pub displacement: f64, } /// Mode `mode` (from 1) of buckling step `step`, scaled so its largest component is `scale` @@ -196,6 +241,7 @@ impl Default for NonlinearStatic { line_search: false, load_path: None, imperfection: None, + arc_length: None, } } } @@ -205,6 +251,22 @@ impl NonlinearStatic { let fraction = |what: &str, v: f64| { if v > 0.0 && v <= 1.0 { Ok(v) } else { Err(Error::Invalid(format!("the {what} is a fraction of the step in (0, 1], got {v}"))) } }; + if let Some(arc) = &self.arc_length { + if !(arc.max_load_factor.is_finite() && arc.max_load_factor > 0.0) { + return Err(Error::Invalid(format!("an arc-length step's maximum load factor is a positive number, got {}", arc.max_load_factor))); + } + if let Some(l) = &arc.displacement_limit { + if !(1..=3).contains(&l.direction) { + return Err(Error::Invalid(format!("a displacement limit's direction is 1, 2 or 3, got {}", l.direction))); + } + if !(l.displacement.is_finite() && l.displacement != 0.0) { + return Err(Error::Invalid(format!("a displacement limit is a nonzero number, got {}", l.displacement))); + } + } + if self.load_path.is_some() { + return Err(Error::Invalid("an arc-length step finds its own load factor; it takes no load path".into())); + } + } if let Some(imp) = &self.imperfection { if !(1..=MAX_MODES).contains(&imp.mode) { return Err(Error::Invalid(format!("an imperfection takes mode 1 to {MAX_MODES}, got {}", imp.mode))); diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs index 4651d3f..3f9cc31 100644 --- a/crates/solver/src/nonlinear.rs +++ b/crates/solver/src/nonlinear.rs @@ -98,7 +98,9 @@ use fieldcraft_fem::finite::{Kinematics, Temperature}; use fieldcraft_fem::plane::Theory; use fieldcraft_fem::quadrature::Rule; use fieldcraft_fem::{DofKind, DofSet, ElementType, PointStress, face, von_mises}; -use fieldcraft_model::{BoundaryCondition, Extent, Hardening, Imperfection, Load, Model, NonlinearStatic, Plastic, Procedure, SectionKind, Step}; +use fieldcraft_model::{ + ArcConstraint, ArcLength, BoundaryCondition, Extent, Hardening, Imperfection, Load, Model, NonlinearStatic, Plastic, Procedure, SectionKind, Step, +}; use fieldcraft_parallel as parallel; use serde::Serialize; @@ -149,6 +151,35 @@ pub struct Increment { /// Why it was cut back. #[serde(skip_serializing_if = "Option::is_none")] pub cutback: Option, + /// Negative eigenvalues of the tangent at the converged state (arc-length steps, which + /// factor it for the next predictor). + #[serde(skip_serializing_if = "is_zero")] + pub negative_pivots: usize, +} + +/// Where the tangent's inertia changes along an arc-length path. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)] +#[serde(rename_all = "snake_case")] +pub enum CriticalKind { + /// The load factor turns (snap-through). + Limit, + /// The load factor goes on: another branch crosses this one. + Bifurcation, +} + +/// A critical point of an arc-length path: the tangent's count of negative eigenvalues changes +/// between two converged states. +#[derive(Debug, Clone, Copy, PartialEq, Serialize)] +pub struct CriticalPoint { + pub kind: CriticalKind, + /// The increment whose converged state has the new count. + pub increment: usize, + /// For a limit point the extreme load factor of the converged states around it, for a + /// bifurcation the last one before it. + pub load_factor: f64, + /// Negative eigenvalues before and after. + pub negative_before: usize, + pub negative_after: usize, } /// The displacements at the end of a converged increment. @@ -168,6 +199,8 @@ pub struct History { pub frames: Vec, /// Newton iterations over the step. pub iterations: usize, + /// Arc-length steps: where the tangent's inertia changes along the path. + pub critical_points: Vec, } impl History { @@ -271,6 +304,14 @@ fn half(k: &[f64], m: usize, sign: f64) -> Vec { (0..m * m).map(|i| 0.5 * (k.get(i).copied().unwrap_or(0.0) + sign * k.get((i % m) * m + i / m).copied().unwrap_or(0.0))).collect() } +/// T_fᵀ K T_p (s u_p): the free equations' forces of the prescribed motions scaled by `s`. +fn prescribed_force(red: &Reduction, k: &Csr, s: f64) -> Vec { + let jump = red.expand_scaled(&vec![0.0; red.neq()], s); + let mut kj = vec![0.0; jump.len()]; + k.mul(&jump, &mut kj); + red.restrict(&kj) +} + /// The out-of-balance forces over the free equations, T_fᵀ (f_ext − f_int). fn free_residual(red: &Reduction, e: &Evaluation) -> Vec { let r: Vec = e.external.iter().zip(&e.internal).map(|(a, b)| a - b).collect(); @@ -402,6 +443,21 @@ impl Problem<'_> { Ok(Evaluation { internal, external, tangent: k, history }) } + /// The loads at load factor 1 at displacements `u` (dead loads and follower pressures on + /// the current faces): ∂f_ext/∂λ. + fn reference_load(&self, u: &[f64]) -> Result> { + let mut external = self.dead.clone(); + for f in &self.followers { + let (force, _) = self.follower(f, u, 1.0)?; + for (d, v) in f.dofs.iter().zip(&force) { + if let Some(x) = external.get_mut(*d) { + *x += v; + } + } + } + Ok(external) + } + /// The size of the follower load stiffness's skew part relative to the whole, over the free /// degrees of freedom (Frobenius norms), at `u` and λ; `None` without followers. fn follower_asymmetry(&self, u: &[f64], lambda: f64) -> Result> { @@ -471,6 +527,15 @@ impl Run<'_> { free_residual(self.problem.red, e) } + /// The round-off floor of the force scale, from the first tangent's mean diagonal. + fn seed_force_floor(&mut self, kff: &Csr) { + if self.force_floor.is_none() { + let diagonal = kff.diagonal(); + let mean = diagonal.iter().map(|d| d.abs()).sum::() / diagonal.len().max(1) as f64; + self.force_floor = Some(1e-8 * mean * self.length); + } + } + fn scale(&self) -> f64 { self.force_scale.max(self.force_floor.unwrap_or(0.0)).max(f64::MIN_POSITIVE) } @@ -537,10 +602,7 @@ impl Run<'_> { let red = self.problem.red; let mut e = self.problem.evaluate(&red.expand_scaled(q0, from), to, true, &self.committed)?; let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; - let jump = red.expand_scaled(&vec![0.0; red.neq()], to - from); - let mut kj = vec![0.0; jump.len()]; - k.mul(&jump, &mut kj); - let rf: Vec = free_residual(red, &e).iter().zip(red.restrict(&kj)).map(|(a, b)| a - b).collect(); + 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(); let [kff] = red.reduce([&k])?; drop(k); @@ -598,9 +660,7 @@ impl Run<'_> { // 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 { - let diagonal = kff.diagonal(); - let mean = diagonal.iter().map(|d| d.abs()).sum::() / diagonal.len().max(1) as f64; - self.force_floor = Some(1e-8 * mean * self.length); + self.seed_force_floor(&kff); } drop(k); let factor = self.factor(&kff, first)?; @@ -980,84 +1040,101 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R history.frames.push(Frame { load_factor: lambda, displacement: split(u, &mesh).0 }); }; keep_frame(&mut history, 0.0, &vec![0.0; mesh.dofs()], &mut warnings); - let mut q = vec![0.0; red.neq()]; - // The step fraction t and the load factor λ(t) of the step's load path (λ = t without one). - let (mut t, mut lambda) = (0.0f64, 0.0f64); - let mut size = settings.initial_increment; - let mut fast = 0usize; - let mut last: Option = None; - // 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); - while t < 1.0 { - if history.increments.len() >= max_increments { - return Err(Error::Limit(format!( - "step `{}`: {} increments reached the step's cap at load factor {lambda:.6} (raise max_increments)", - step.name, settings.max_increments - ))); + let (q, lambda, last) = if let Some(arc) = &settings.arc_length { + if problem.field.is_some() { + return Err(Error::NotSupported(format!("temperatures in arc-length steps (step `{}`)", step.name))); } - // An increment ends at the next point of the load path at the latest. - let turn = settings.next_turn(t); - let remaining = turn - t; - let (target, this) = if size >= remaining * (1.0 - 1e-9) { (turn, remaining) } else { (t + size, size) }; - let target_lambda = settings.load_factor(target); - let change = if settings.load_path.is_some() { target_lambda - lambda } else { this }; - let mut iterations = Vec::new(); - let number = history.increments.len() + 1; - // Linear extrapolation of the last converged increment's free motion, in proportion to - // the change of the load factor (so it turns with the load path). - let mut predictor: Option> = previous - .as_ref() - .filter(|(_, prev)| *prev != 0.0) - .map(|(dq, prev): &(Vec, f64)| q.iter().zip(dq).map(|(a, d)| a + change / prev * d).collect()); - // Path-dependent materials start without one from the tangent predictor. - let mut predicted_energy = None; - if predictor.is_none() && path_dependent { - match run.tangent_predictor(&q, lambda, target_lambda) { - Ok(Some((p, energy))) => (predictor, predicted_energy) = (Some(p), Some(energy)), - Ok(None) | Err(Failure::Cutback(_)) => {} + arc_length(&mut run, arc, step, &mut history, &mut warnings, &mut keep_frame)? + } else { + let mut q = vec![0.0; red.neq()]; + // The step fraction t and the load factor λ(t) of the step's load path (λ = t without one). + let (mut t, mut lambda) = (0.0f64, 0.0f64); + let mut size = settings.initial_increment; + let mut fast = 0usize; + let mut last: Option = None; + // 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); + while t < 1.0 { + if history.increments.len() >= max_increments { + return Err(Error::Limit(format!( + "step `{}`: {} increments reached the step's cap at load factor {lambda:.6} (raise max_increments)", + step.name, settings.max_increments + ))); + } + // An increment ends at the next point of the load path at the latest. + let turn = settings.next_turn(t); + let remaining = turn - t; + let (target, this) = if size >= remaining * (1.0 - 1e-9) { (turn, remaining) } else { (t + size, size) }; + let target_lambda = settings.load_factor(target); + let change = if settings.load_path.is_some() { target_lambda - lambda } else { this }; + let mut iterations = Vec::new(); + let number = history.increments.len() + 1; + // Linear extrapolation of the last converged increment's free motion, in proportion to + // the change of the load factor (so it turns with the load path). + let mut predictor: Option> = previous + .as_ref() + .filter(|(_, prev)| *prev != 0.0) + .map(|(dq, prev): &(Vec, f64)| q.iter().zip(dq).map(|(a, d)| a + change / prev * d).collect()); + // Path-dependent materials start without one from the tangent predictor. + let mut predicted_energy = None; + if predictor.is_none() && path_dependent { + match run.tangent_predictor(&q, lambda, target_lambda) { + Ok(Some((p, energy))) => (predictor, predicted_energy) = (Some(p), Some(energy)), + Ok(None) | Err(Failure::Cutback(_)) => {} + Err(Failure::Fatal(e)) => return Err(e), + } + } + match run.increment(&q, predictor.as_deref(), predicted_energy, target_lambda, change.abs(), &mut iterations) { + Ok(c) => { + previous = Some((c.q.iter().zip(&q).map(|(a, b)| a - b).collect(), change)); + let count = iterations.len(); + history.increments.push(Increment { + number, + load_factor: target_lambda, + size: this, + iterations, + converged: true, + cutback: None, + negative_pivots: 0, + }); + (t, lambda) = (target, target_lambda); + q = c.q; + keep_frame(&mut history, lambda, &red.expand_scaled(&q, lambda), &mut warnings); + run.committed.clone_from(&c.evaluation.history); + last = Some(c.evaluation); + fast = if count <= FAST_ITERATIONS { fast + 1 } else { 0 }; + if fast >= 2 { + size *= GROWTH; + } else if count > SLOW_ITERATIONS { + size *= SHRINK; + } + size = size.clamp(settings.min_increment, settings.max_increment); + } Err(Failure::Fatal(e)) => return Err(e), - } - } - match run.increment(&q, predictor.as_deref(), predicted_energy, target_lambda, change.abs(), &mut iterations) { - Ok(c) => { - previous = Some((c.q.iter().zip(&q).map(|(a, b)| a - b).collect(), change)); - let count = iterations.len(); - history.increments.push(Increment { number, load_factor: target_lambda, size: this, iterations, converged: true, cutback: None }); - (t, lambda) = (target, target_lambda); - q = c.q; - keep_frame(&mut history, lambda, &red.expand_scaled(&q, lambda), &mut warnings); - run.committed.clone_from(&c.evaluation.history); - last = Some(c.evaluation); - fast = if count <= FAST_ITERATIONS { fast + 1 } else { 0 }; - if fast >= 2 { - size *= GROWTH; - } else if count > SLOW_ITERATIONS { - size *= SHRINK; - } - size = size.clamp(settings.min_increment, settings.max_increment); - } - Err(Failure::Fatal(e)) => return Err(e), - Err(Failure::Cutback(reason)) => { - history.increments.push(Increment { - number, - load_factor: target_lambda, - size: this, - iterations, - converged: false, - cutback: Some(reason.clone()), - }); - fast = 0; - size = this * CUTBACK; - if size < settings.min_increment * (1.0 - 1e-9) { - return Err(Error::NotConverged(format!( - "step `{}` stopped at load factor {lambda:.6} (step fraction {t:.6}): the increment to {target_lambda:.6} failed ({reason}) and a smaller one would be below the minimum increment {}", - step.name, settings.min_increment - ))); + Err(Failure::Cutback(reason)) => { + history.increments.push(Increment { + number, + load_factor: target_lambda, + size: this, + iterations, + converged: false, + cutback: Some(reason.clone()), + negative_pivots: 0, + }); + fast = 0; + size = this * CUTBACK; + if size < settings.min_increment * (1.0 - 1e-9) { + return Err(Error::NotConverged(format!( + "step `{}` stopped at load factor {lambda:.6} (step fraction {t:.6}): the increment to {target_lambda:.6} failed ({reason}) and a smaller one would be below the minimum increment {}", + step.name, settings.min_increment + ))); + } } } } - } + (q, lambda, last) + }; history.iterations = run.iterations; let u = red.expand_scaled(&q, lambda); let evaluation = match last { @@ -1067,6 +1144,7 @@ fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> R if let Some((at, negative)) = history .increments .iter() + .filter(|_| settings.arc_length.is_none()) .filter(|i| i.converged) .filter_map(|i| i.iterations.last().map(|it| (i.load_factor, it.negative_pivots))) .find(|(_, n)| *n > 0) @@ -1197,3 +1275,285 @@ fn finish( nonlinear: Some(history), }) } + +/// What one arc-length increment starts from: the tangent at the converged state, factored, and +/// the tangent displacement of the unit load. +struct ArcStart { + factor: Symmetric, + load_displacement: Vec, +} + +impl Run<'_> { + /// The tangent at (q, λ) factored (Cholesky, or LDLᵀ past a critical point), the unit load + /// ∂r/∂λ at that state over the free equations (the loads at factor 1, less K T_p u_p for + /// prescribed motions), and the evaluation's residual. + fn arc_tangent(&mut self, e: &mut Evaluation, u: &[f64], first: bool) -> std::result::Result<(Symmetric, Vec), Failure> { + let red = self.problem.red; + let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; + let mut g = red.restrict(&self.problem.reference_load(u)?); + if red.largest_prescribed() > 0.0 { + g.iter_mut().zip(prescribed_force(red, &k, 1.0)).for_each(|(a, b)| *a -= b); + } + let [kff] = red.reduce([&k])?; + drop(k); + if first { + self.seed_force_floor(&kff); + } + let factor = self.factor(&kff, first)?; + self.factor_entries = factor.entries(); + Ok((factor, g)) + } + + /// The converged state's tangent (from its evaluation `e`, which keeps its forces) and the + /// unit load's tangent displacement. + fn arc_start(&mut self, e: &mut Evaluation, q: &[f64], lambda: f64, first: bool) -> std::result::Result { + let (factor, mut g) = self.arc_tangent(e, &self.problem.red.expand_scaled(q, lambda), first)?; + factor.solve(&mut g)?; + if g.iter().any(|v| !v.is_finite()) || norm_inf(&g) == 0.0 { + return Err(Failure::Fatal(Error::Invalid("the step's loads do not move the model, so an arc length cannot be measured".into()))); + } + Ok(ArcStart { factor, load_displacement: g }) + } + + /// One arc-length increment from the converged (q0, λ0) (M. A. Crisfield, *A fast + /// incremental/iterative solution procedure that handles "snap-through"*, Computers & + /// Structures 13:55–62, 1981): the constraint Δqᵀ Δq + w Δλ² = Δl² (w = 0 cylindrical, w = + /// ‖q_f‖² spherical) with the load factor an unknown. Predictor along the tangent, its sign + /// by the angle with the last increment (M. Feng, D. Perić, D. R. J. Owen, *A new criterion for + /// determination of initial loading parameter in arc-length methods*, Computers & Structures + /// 58(3):479–485, 1996); each correction δq = δq_r + δλ δq_g solves the quadratic for δλ and + /// takes the root whose increment is closest in direction to the previous iterate's (the + /// angle criterion, Crisfield 1981). + #[allow(clippy::too_many_arguments)] + fn arc_increment( + &mut self, + q0: &[f64], + lambda0: f64, + start: &ArcStart, + arc: f64, + weight: f64, + previous: Option<&(Vec, f64)>, + iterations: &mut Vec, + ) -> std::result::Result<(Vec, f64, Evaluation), Failure> { + let tol = self.settings.tolerances; + let qg = &start.load_displacement; + let sign = match previous { + Some((dq, dl)) if dot(dq, qg) + weight * dl < 0.0 => -1.0, + _ => 1.0, + }; + let mut dlam = sign * arc / (dot(qg, qg) + weight).sqrt(); + let mut dq: Vec = qg.iter().map(|g| dlam * g).collect(); + let mut correction = f64::INFINITY; + let mut first_energy: Option = None; + let max_iterations = usize::try_from(self.settings.max_iterations).unwrap_or(usize::MAX); + let at = |dq: &[f64]| -> Vec { q0.iter().zip(dq).map(|(a, b)| a + b).collect() }; + for it in 0..=max_iterations { + let q = at(&dq); + let mut e = self.state(&q, lambda0 + dlam, true)?; + let r = self.residual(&e); + let force = norm_inf(&r) / self.scale(); + if !force.is_finite() { + return Err(Failure::Cutback("the residual is not finite".into())); + } + if it > 0 && ((force <= tol.force && correction <= tol.displacement) || force <= ROUND_OFF * tol.force) { + return Ok((q, lambda0 + dlam, e)); + } + if it == max_iterations { + break; + } + if self.iterations >= self.limits.max_newton_iterations { + return Err(Failure::Fatal(Error::Limit(format!( + "the step took {} Newton iterations, the limit (raise max_newton_iterations or loosen the tolerances)", + self.iterations + )))); + } + self.iterations += 1; + let u = self.problem.red.expand_scaled(&q, lambda0 + dlam); + let (factor, mut dqg) = self.arc_tangent(&mut e, &u, false)?; + let mut dqr = r.clone(); + factor.solve(&mut dqr)?; + factor.solve(&mut dqg)?; + // (Δq + δq_r + δλ δq_g)² + w (Δλ + δλ)² = Δl². + let base: Vec = dq.iter().zip(&dqr).map(|(a, b)| a + b).collect(); + let a = dot(&dqg, &dqg) + weight; + let b = 2.0 * (dot(&dqg, &base) + weight * dlam); + let c = dot(&base, &base) + weight * dlam * dlam - arc * arc; + let disc = b * b - 4.0 * a * c; + if !(disc >= 0.0) || !(a > 0.0) { + return Err(Failure::Cutback("the arc-length constraint has no real root (the arc is too long here)".into())); + } + let root = disc.sqrt(); + // Stable roots of a x² + b x + c (no cancellation). + let qr = -0.5 * (b + b.signum() * root); + let roots = if qr == 0.0 { [0.0, 0.0] } else { [qr / a, c / qr] }; + let cosine = |dl: f64| dot(&dq, &base) + dl * dot(&dq, &dqg) + weight * dlam * (dlam + dl); + let dl = if cosine(roots[0]) >= cosine(roots[1]) { roots[0] } else { roots[1] }; + let step: Vec = dqr.iter().zip(&dqg).map(|(r, g)| r + dl * g).collect(); + if step.iter().any(|v| !v.is_finite()) { + return Err(Failure::Cutback("the arc-length correction is not finite".into())); + } + dq.iter_mut().zip(&step).for_each(|(a, s)| *a += s); + dlam += dl; + correction = norm_inf(&step) / norm_inf(&dq).max(1e-12 * self.length); + let energy = dot(&step, &r).abs(); + let first = *first_energy.get_or_insert(energy); + iterations.push(Iteration { + force, + displacement: correction, + energy: if first > 0.0 { energy / first } else { 0.0 }, + line_search: 1.0, + negative_pivots: factor.inertia().negative, + }); + } + Err(Failure::Cutback(format!("no convergence in {} iterations", self.settings.max_iterations))) + } +} + +/// Follow a step's equilibrium path by arc length (Riks 1979, Crisfield 1981): the load factor +/// is an unknown, so the path goes through limit points. The first increment's load factor +/// change is the step's initial increment, which sets the arc length; increments grow and +/// shrink as under load control (the arc length in proportion, within the minimum and maximum +/// increments' arcs). The tangent at every converged state is factored for the next predictor +/// (Cholesky, LDLᵀ when indefinite); its negative eigenvalues are recorded per increment, and +/// where their count changes a critical point is reported: a limit point when the load factor +/// turns there, a bifurcation otherwise. The step ends at the maximum load factor, at the +/// displacement limit, or at its increment cap (with a warning). +#[allow(clippy::type_complexity)] +fn arc_length( + run: &mut Run, + arc: &ArcLength, + step: &Step, + history: &mut History, + warnings: &mut Vec, + keep_frame: &mut dyn FnMut(&mut History, f64, &[f64], &mut Vec), +) -> Result<(Vec, f64, Option)> { + let settings = run.settings; + let red = run.problem.red; + let mesh = run.problem.mesh; + let limit = match &arc.displacement_limit { + Some(l) => { + let slot = + mesh.node(l.node).ok_or_else(|| Error::Invalid(format!("the displacement limit's node {} has no degrees of freedom", l.node)))?; + let dofs = translations(mesh, slot)?; + let dof = dofs + .get(usize::from(l.direction).saturating_sub(1)) + .copied() + .ok_or_else(|| Error::Invalid(format!("the displacement limit's node {} has no direction {}", l.node, l.direction)))?; + Some((dof, l.displacement.abs())) + } + None => None, + }; + let fatal = |f: Failure| match f { + Failure::Fatal(e) => e, + Failure::Cutback(m) => Error::NotConverged(format!("step `{}`: {m}", step.name)), + }; + let mut q = vec![0.0; red.neq()]; + let mut lambda = 0.0f64; + let mut initial = run.state(&q, lambda, true).map_err(fatal)?; + let mut start = run.arc_start(&mut initial, &q, lambda, true).map_err(fatal)?; + // States: load factor and negative eigenvalues, from the undeformed one. + let mut states: Vec<(usize, f64, usize)> = vec![(0, 0.0, start.factor.inertia().negative)]; + let scale = (dot(&start.load_displacement, &start.load_displacement)).sqrt(); + let weight = if arc.constraint == ArcConstraint::Spherical { scale * scale } else { 0.0 }; + let unit = (scale * scale + weight).sqrt(); + let mut size = settings.initial_increment * unit; + let (smallest, largest) = (settings.min_increment * unit, settings.max_increment * unit); + let mut previous: Option<(Vec, f64)> = None; + let mut fast = 0usize; + let mut last: Option = None; + let max_increments = usize::try_from(settings.max_increments).unwrap_or(usize::MAX); + loop { + if history.increments.len() >= max_increments { + warnings.push(format!( + "the arc-length step `{}` ended at its cap of {} increments, at load factor {lambda:.6}", + step.name, settings.max_increments + )); + break; + } + let number = history.increments.len() + 1; + let mut iterations = Vec::new(); + match run.arc_increment(&q, lambda, &start, size, weight, previous.as_ref(), &mut iterations) { + Ok((q1, lambda1, mut evaluation)) => { + let change = lambda1 - lambda; + previous = Some((q1.iter().zip(&q).map(|(a, b)| a - b).collect(), change)); + let count = iterations.len(); + (q, lambda) = (q1, lambda1); + run.committed.clone_from(&evaluation.history); + let u = red.expand_scaled(&q, lambda); + keep_frame(history, lambda, &u, warnings); + start = run.arc_start(&mut evaluation, &q, lambda, false).map_err(fatal)?; + let negative = start.factor.inertia().negative; + history.increments.push(Increment { + number, + load_factor: lambda, + size: change, + iterations, + converged: true, + cutback: None, + negative_pivots: negative, + }); + states.push((number, lambda, negative)); + last = Some(evaluation); + fast = if count <= FAST_ITERATIONS { fast + 1 } else { 0 }; + if fast >= 2 { + size *= GROWTH; + } else if count > SLOW_ITERATIONS { + size *= SHRINK; + } + size = size.clamp(smallest, largest); + let reached = limit.is_some_and(|(dof, value)| u.get(dof).is_some_and(|x| x.abs() >= value)); + if lambda >= arc.max_load_factor || reached { + break; + } + } + Err(Failure::Fatal(e)) => return Err(e), + Err(Failure::Cutback(reason)) => { + history.increments.push(Increment { + number, + load_factor: lambda, + size: 0.0, + iterations, + converged: false, + cutback: Some(reason.clone()), + negative_pivots: 0, + }); + fast = 0; + size *= CUTBACK; + if size < smallest * (1.0 - 1e-9) { + return Err(Error::NotConverged(format!( + "arc-length step `{}` stopped at load factor {lambda:.6}: an increment failed ({reason}) and a shorter arc would be below the minimum increment", + step.name + ))); + } + } + } + } + history.critical_points = critical_points(&states); + Ok((q, lambda, last)) +} + +/// Critical points from the converged states' (increment, load factor, negative eigenvalues): +/// where the count changes, a limit point if the load factor turns within a state of it, a +/// bifurcation otherwise. +fn critical_points(states: &[(usize, f64, usize)]) -> Vec { + let mut out = Vec::new(); + for k in 1..states.len() { + let (Some(&(_, l0, n0)), Some(&(increment, _, n1))) = (states.get(k - 1), states.get(k)) else { continue }; + if n0 == n1 { + continue; + } + let around: Vec = states.get(k.saturating_sub(2)..(k + 2).min(states.len())).unwrap_or_default().iter().map(|s| s.1).collect(); + let turns = around.windows(3).any(|w| (w[1] - w[0]) * (w[2] - w[1]) < 0.0); + let (kind, load_factor) = if turns { + // The extreme load factor next to the change. + let rising = around.windows(2).next().is_some_and(|w| w[1] >= w[0]); + let extreme = + if rising { around.iter().copied().fold(f64::NEG_INFINITY, f64::max) } else { around.iter().copied().fold(f64::INFINITY, f64::min) }; + (CriticalKind::Limit, extreme) + } else { + (CriticalKind::Bifurcation, l0) + }; + out.push(CriticalPoint { kind, increment, load_factor, negative_before: n0, negative_after: n1 }); + } + out +} diff --git a/crates/solver/tests/arc_length.rs b/crates/solver/tests/arc_length.rs new file mode 100644 index 0000000..a273bf4 --- /dev/null +++ b/crates/solver/tests/arc_length.rs @@ -0,0 +1,137 @@ +//! Arc-length (Riks–Crisfield) control through limit points: the snap-through of a shallow +//! sinusoidal arch of solids under a sinusoidal dead load, against the one-mode closed form of +//! shallow arch theory (exact for this load and shape on the symmetric path), the limit points' +//! load factors and the tangent's inertia along the path, and what arc-length steps refuse. +//! +//! Closed form: an arch y₀ = a sin(πx/L) on pins at a fixed span, bending stiffness EI, axial EA, +//! r² = I/A, under q = q₀ sin(πx/L) deflects in the same shape, w = c sin(πx/L), with +//! q₀ = EI k⁴ [c + (a − c)(2ac − c²)/(4r²)], k = π/L (the energy of bending and of the axial +//! shortening ½∫(y₀′² − y′²)dx, S. P. Timoshenko, J. M. Gere, *Theory of Elastic Stability*, 2nd +//! ed., 1961, §7.6, and Y. C. Fung, A. Kaplan, NACA TN 2840, 1952). Its limit points are at +//! c = a ∓ r √((a²/r² − 4)/3); for a/r = 3.5 the antisymmetric mode never bifurcates first (the +//! axial force stays below the second Euler load 4EIk², reached only at a/r = 4). + +mod common; + +use common::*; +use fieldcraft_model::*; +use fieldcraft_solver::nonlinear::CriticalKind; +use fieldcraft_solver::{Limits, solve_nonlinear}; + +const E: f64 = 1.0e4; +const SPAN: f64 = 100.0; +const H: f64 = 1.0; + +fn rise() -> f64 { + 3.5 * H / 12f64.sqrt() +} + +/// q₀ for a midspan deflection c (positive down). +fn load(c: f64) -> f64 { + let (a, r2, k) = (rise(), H * H / 12.0, std::f64::consts::PI / SPAN); + E * H.powi(3) / 12.0 * k.powi(4) * (c + (a - c) * (2.0 * a * c - c * c) / (4.0 * r2)) +} + +fn limit_points() -> (f64, f64) { + let (a, r) = (rise(), H / 12f64.sqrt()); + let d = r * ((a * a / (r * r) - 4.0) / 3.0).sqrt(); + (a - d, a + d) +} + +/// The arch (40 × 2 × 1 hex20, unit width), its load scaled to the closed-form limit load, and +/// the label of the midspan node at mid-height. +fn arch(arc: Option) -> (Grid, u32, f64) { + let a = rise(); + let mut g = grid([40, 2, 1], ElementKind::Hex20, |p| { + let x = SPAN * p[0]; + [x, a * (std::f64::consts::PI * x / SPAN).sin() + H * (p[1] - 0.5), p[2]] + }); + let reference = load(limit_points().0); + let settings = NonlinearStatic { arc_length: arc, initial_increment: 0.1, max_increments: 200, ..NonlinearStatic::default() }; + g.define(iso(E, 0.0), None, Formulation::default(), Procedure::NonlinearStatic(settings)); + g.node_set("pins", |p| (near(p[0], 0.0) || near(p[0], 1.0)) && near(p[1], 0.5)); + g.bc("pins", BoundaryCondition::Displacement { target: Target::Set("pins".into()), values: [Some(0.0), Some(0.0), None, None, None, None] }); + g.node_set("width", |_| true); + g.bc("width", BoundaryCondition::Displacement { target: Target::Set("width".into()), values: [None, None, Some(0.0), None, None, None] }); + // q₀ sin(πx/L) per unit length, piecewise constant per element on the top face (dead). + let top = g.faces_where(|p| near(p[1], 1.0)); + for f in top { + let nodes = &g.model.element(f.element).unwrap().nodes; + let xc = nodes.iter().take(8).map(|n| g.model.node(*n).unwrap().xyz[0]).sum::() / 8.0; + let name = format!("q{}", f.element); + g.model.add_set(&name, Set { faces: vec![f], ..Set::default() }).unwrap(); + let total = -reference * (std::f64::consts::PI * xc / SPAN).sin() * SPAN / 40.0; + g.load(&name, Load::Force { target: Target::Set(name.clone()), vector: [0.0, total, 0.0], distribution: Distribution::Total }); + } + let mid = g.nodes_where(|p| near(p[0], 0.5) && near(p[1], 0.5) && near(p[2], 0.0))[0]; + (g, mid, reference) +} + +#[test] +fn shallow_arch_snaps_through() { + let a = rise(); + let (_, mid, _) = arch(None); + let arc = ArcLength { + displacement_limit: Some(DisplacementLimit { node: mid, direction: 2, displacement: -2.2 * a }), + max_load_factor: 10.0, + ..ArcLength::default() + }; + let (g, _, reference) = arch(Some(arc)); + let r = solve_nonlinear(&g.model, Some("Step"), &Limits::default()).unwrap(); + let h = r.nonlinear.as_ref().unwrap(); + let at = r.nodes.binary_search(&mid).unwrap(); + // Every converged state on the closed-form path. + let mut worst = 0.0f64; + for frame in &h.frames { + let c = -frame.displacement[at][1]; + let err = (frame.load_factor * reference - load(c)).abs() / reference; + worst = worst.max(err); + } + let (_, c_valley) = limit_points(); + eprintln!("{} increments; worst path error {:.3} % of the limit load", h.increments.len(), 100.0 * worst); + for p in &h.critical_points { + eprintln!(" {:?} at increment {}: λ = {:.5} ({} → {} negative)", p.kind, p.increment, p.load_factor, p.negative_before, p.negative_after); + } + assert!(worst < 0.02, "{worst}"); + assert!(-r.displacement[at][1] >= 2.2 * a, "ended at the displacement limit"); + // The limit points: the peak (λ = 1 by construction) and the valley, from the inertia. + let [peak, valley] = h.critical_points.as_slice() else { panic!("{:?}", h.critical_points) }; + assert_eq!((peak.kind, peak.negative_before, peak.negative_after), (CriticalKind::Limit, 0, 1)); + assert_eq!((valley.kind, valley.negative_before, valley.negative_after), (CriticalKind::Limit, 1, 0)); + assert!((peak.load_factor - 1.0).abs() < 0.02, "{}", peak.load_factor); + assert!((valley.load_factor - load(c_valley) / reference).abs() < 0.02, "{} vs {}", valley.load_factor, load(c_valley) / reference); + // Negative pivots exactly between the limit points: every converged state's count follows + // the sign of dλ along the path. + let states: Vec<(f64, usize)> = h.increments.iter().filter(|i| i.converged).map(|i| (i.load_factor, i.negative_pivots)).collect(); + for w in states.windows(3) { + let falling = w[2].0 < w[1].0 && w[1].0 < w[0].0; + let rising = w[2].0 > w[1].0 && w[1].0 > w[0].0; + if falling { + assert_eq!(w[1].1, 1, "{states:?}"); + } + if rising { + assert_eq!(w[1].1, 0, "{states:?}"); + } + } +} + +#[test] +fn arc_length_steps_refuse_what_they_cannot_do() { + let arc = ArcLength::default(); + let mut settings = NonlinearStatic { arc_length: Some(arc.clone()), load_path: Some(vec![[0.0, 0.0], [1.0, 1.0]]), ..NonlinearStatic::default() }; + assert!(settings.validate().is_err(), "a load path"); + settings.load_path = None; + settings.arc_length = Some(ArcLength { max_load_factor: f64::NAN, ..arc.clone() }); + assert!(settings.validate().is_err()); + settings.arc_length = Some(ArcLength { displacement_limit: Some(DisplacementLimit { node: 1, direction: 7, displacement: 1.0 }), ..arc.clone() }); + assert!(settings.validate().is_err()); + // A node that does not exist; no load at all. + let (mut g, _, _) = + arch(Some(ArcLength { displacement_limit: Some(DisplacementLimit { node: 999_999, direction: 2, displacement: 1.0 }), ..arc })); + assert!(solve_nonlinear(&g.model, Some("Step"), &Limits::default()).is_err()); + let loads: Vec = g.model.step("Step").unwrap().loads.iter().cloned().collect(); + for l in loads { + g.model.set_load_active("Step", &l, false).unwrap(); + } + assert!(solve_nonlinear(&g.model, Some("Step"), &Limits::default()).is_err()); +} diff --git a/docs/fem.md b/docs/fem.md index 882720a..4e11f45 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -313,6 +313,26 @@ displacement (at ν = 0.4999 every increment used to be cut back to the cap). On tangent of a model that can move freely is an error. The tangent predictor of path-dependent materials keeps Cholesky and is skipped when the tangent is indefinite. +**Arc length.** `"arc_length": {"constraint": "cylindrical"|"spherical", "max_load_factor": 1, +"displacement_limit": {"node", "direction", "displacement"}}` makes the load factor an unknown +(Riks 1979, Crisfield 1981): each increment keeps ‖Δq‖² + w Δλ² = Δl² (w = 0 cylindrical, w the +first tangent displacement's ‖q_f‖² for the spherical arc), predicts along the tangent with the +sign of its angle to the last increment (Feng, Perić & Owen 1996), and corrects with +δq = δq_r + δλ δq_g, δλ the root of the constraint's quadratic whose increment turns least +(Crisfield's angle criterion). The first increment's load factor change is `initial_increment`, +which sets the arc; arcs grow, shrink and cut back as increments do, within the minimum and +maximum increments' arcs. Prescribed motions scale with λ (∂r/∂λ includes −K T_p u_p) and +follower pressures with their current faces. The tangent at every converged state is factored +for the next predictor (LDLᵀ when indefinite); its negative eigenvalues are recorded per +increment (`negative_pivots`), and where their count changes the step reports a **critical +point**: a limit point where the load factor turns, a bifurcation otherwise +(`History::critical_points`, in the `solve.nonlinear` summary). The step ends at the maximum load +factor, at the displacement limit, or at its increment cap (a warning). Refused: a load path, +temperatures. Decks: Abaqus `*STATIC, RIKS` (cylindrical only); LS-DYNA refuses arc-length steps +by name. Verified by the snap-through of a shallow sinusoidal arch of hex20 solids +(`crates/solver/tests/arc_length.rs`, `tests/oracle/arc-shallow-arch`): the path within 1.9 % of +the one-mode closed form, both limit points, and the inertia 0 → 1 → 0 exactly between them. + **Results**: the final displacements and reactions, the **Cauchy stress** and the **logarithmic strain** (ln V, engineering shears; the small strain without nonlinear geometry) at the integration points and averaged at the nodes, the **equivalent plastic strain** at the diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index e59ea68..f7b961a 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -1244,6 +1244,14 @@ update are written from (cited again in the module docs): indefinite sparse symmetric linear equations*, ACM TOMS 9(3):302–325, 1983; J. W. H. Liu, *The multifrontal method for sparse matrix solution: theory and practice*, SIAM Review 34(1):82–109, 1992. +- Arc-length: E. Riks, *An incremental approach to the solution of snapping and buckling + problems*, Int. J. Solids Struct. 15(7):529–551, 1979; M. A. Crisfield, *A fast + incremental/iterative solution procedure that handles "snap-through"*, Computers & Structures + 13:55–62, 1981; E. Ramm, *Strategies for tracing the nonlinear response near limit points*, in + *Nonlinear Finite Element Analysis in Structural Mechanics*, Springer 1981, 63–89; the + predictor's sign, M. Feng, D. Perić, D. R. J. Owen, Computers & Structures 58(3):479–485, 1996. + Snap-through of a shallow sinusoidal arch: Y. C. Fung, A. Kaplan, NACA TN 2840, 1952; + Timoshenko & Gere 1961, §7.6. - Iterative refinement: N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM 2002, §12.2. - Linear buckling and the initial-stress stiffness: K.-J. Bathe, *Finite Element Procedures*,