From 2183fdbe810e8c733365bca8f3c56784c2596d28 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 03:19:31 -0400 Subject: [PATCH] =?UTF-8?q?solver:=20nonlinear=20static=20analysis=20by=20?= =?UTF-8?q?Newton=E2=80=93Raphson=20with=20automatic=20increments?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit solve.nonlinear follows a nonlinear static step in load increments: dead loads and prescribed displacements ramp with the load factor, a pressure follows the deformed faces (its load stiffness symmetrised so the supernodal Cholesky stays), the consistent tangent is refactored every iteration with an optional line search, and an increment converges on force, displacement and energy together. Increments start from the extrapolated last one, grow after fast convergence and are cut back on inversion, a non-finite or growing residual, an indefinite tangent or too many iterations; the step stops cleanly at the minimum increment, the increment cap or the Newton iteration limit. Elements are assembled in parallel in label order (same bits on any thread count). Results carry Cauchy stress, logarithmic strain, the increment history and a displacement frame per increment, which results.frame and the results panel show. Trusses, shells, beams, plane stress, incompatible modes, temperature and centrifugal loads are refused by name; couplings warn that they stay linear. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --- crates/engine/src/catalogue.rs | 8 +- crates/engine/src/cmd/mod.rs | 1 + crates/engine/src/cmd/results.rs | 15 +- crates/engine/src/cmd/solve.rs | 89 ++- crates/engine/src/results.rs | 15 +- crates/solver/src/analysis.rs | 32 +- crates/solver/src/constraints.rs | 44 ++ crates/solver/src/lib.rs | 12 +- crates/solver/src/nonlinear.rs | 897 ++++++++++++++++++++++++++++ crates/solver/src/system.rs | 37 +- crates/solver/tests/nonlinear.rs | 435 ++++++++++++++ crates/ui-egui/src/results.rs | 20 +- crates/ui-egui/src/results_panel.rs | 11 + docs/features.md | 6 +- 14 files changed, 1579 insertions(+), 43 deletions(-) create mode 100644 crates/solver/src/nonlinear.rs create mode 100644 crates/solver/tests/nonlinear.rs diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index 360ce63..420346a 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -301,7 +301,7 @@ pub static CATALOGUE: &[Capability] = &[ "Material", Properties, "Materials", - "Elastic (isotropic, orthotropic) material with density", + "Elastic (isotropic, orthotropic) material with density; neo-Hookean hyperelasticity for nonlinear steps", &["model.material.edit", "model.material.delete"], ), ready("model.material.plasticity", "Plasticity", Properties, "Materials", "Von Mises plasticity: bilinear, hardening table or Johnson–Cook"), @@ -343,7 +343,7 @@ pub static CATALOGUE: &[Capability] = &[ "Step", Steps, "Steps", - "A linear static or frequency step with its output", + "A linear or nonlinear static, frequency, heat or explicit step with its output", &["model.step.edit", "model.step.delete"], ), ready("model.step.activate", "Activate", Steps, "Steps", "Choose which loads and boundary conditions act in a step"), @@ -437,8 +437,8 @@ pub static CATALOGUE: &[Capability] = &[ "Static", Solve, "Implicit", - "Linear static (or coupled temperature–displacement) solve of a step with FieldCraft's own solver; heat transfer in the drop-down", - &["solve.heat"], + "Linear static (or coupled temperature–displacement) solve of a step with FieldCraft's own solver; nonlinear static and heat transfer in the drop-down", + &["solve.nonlinear", "solve.heat"], ), ready("solve.frequency", "Frequency", Solve, "Implicit", "Natural frequencies and mode shapes of a frequency step"), ready( diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index c12d67c..694c372 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -95,6 +95,7 @@ pub static ALL: &[CommandSpec] = &[ deck::WRITE, solve::RUN, solve::STATIC, + solve::NONLINEAR, solve::FREQUENCY, solve::HEAT, solve::EXPLICIT, diff --git a/crates/engine/src/cmd/results.rs b/crates/engine/src/cmd/results.rs index d8f15e6..975c368 100644 --- a/crates/engine/src/cmd/results.rs +++ b/crates/engine/src/cmd/results.rs @@ -43,7 +43,7 @@ pub const ANIMATE: CommandSpec = CommandSpec { pub const FRAME: CommandSpec = CommandSpec { id: "results.frame", label: "Field Frame", - params: r#"{"step": "Impact", "frame": 3, "time": 0.002} — shows a field frame of a solved explicit step, by index from 0 or the frame nearest a time (neither: the last); the deformed shape and contours follow it"#, + params: r#"{"step": "Impact", "frame": 3, "time": 0.002} — shows a field frame of a solved explicit step, by index from 0 or the frame nearest a time (neither: the last), or of a nonlinear static step (one per converged increment, `time` its load factor; the deformed shape follows it, its contours show the final state); the deformed shape and contours follow it"#, kind: CommandKind::Session, run: frame, }; @@ -79,21 +79,23 @@ fn solved_list(s: &Session) -> String { enum Want { Any, Frequency, - Explicit, + /// Steps with frames: explicit dynamics and nonlinear static (one per converged increment). + Frames, } impl Want { fn fits(self, s: &Session, step: &str) -> bool { s.latest_result(step).is_some_and(|(r, _)| match (self, r.kind()) { - (Want::Any, _) | (Want::Frequency, Kind::Frequency(_)) | (Want::Explicit, Kind::Explicit(_)) => true, - (Want::Frequency | Want::Explicit, _) => false, + (Want::Any, _) | (Want::Frequency, Kind::Frequency(_)) | (Want::Frames, Kind::Explicit(_)) => true, + (Want::Frames, Kind::Static(r)) => r.nonlinear.is_some(), + (Want::Frequency | Want::Frames, _) => false, }) } fn name(self) -> &'static str { match self { Want::Any => "solved", Want::Frequency => "frequency", - Want::Explicit => "explicit", + Want::Frames => "explicit or nonlinear static", } } } @@ -222,9 +224,10 @@ fn animate(s: &mut Session, p: &Value) -> Result { } fn frame(s: &mut Session, p: &Value) -> Result { - let step = step_for(s, p, Want::Explicit)?; + let step = step_for(s, p, Want::Frames)?; let times: Vec = match s.latest_result(&step).map(|(r, _)| r.kind()) { Some(Kind::Explicit(e)) => e.frames.iter().map(|f| f.time).collect(), + Some(Kind::Static(r)) => r.nonlinear.as_ref().map(|h| h.load_factors()).unwrap_or_default(), _ => Vec::new(), }; let last = times.len().checked_sub(1).ok_or_else(|| EngineError::Params(format!("step `{step}` has no field frames")))?; diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index a84f39c..e642569 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -7,7 +7,8 @@ use fieldcraft_mesh::surface::{BoundarySurface, SurfaceOptions}; use fieldcraft_model::FluidDomain; use fieldcraft_parallel as parallel; use fieldcraft_solver::{ - FrequencyOptions, Limits, Location, MassKind, StaticResult, StepResult as ImplicitResult, solve_frequency, solve_heat, solve_static, + FrequencyOptions, Limits, Location, MassKind, StaticResult, StepResult as ImplicitResult, solve_frequency, solve_heat, solve_nonlinear, + solve_static, }; use serde_json::{Value, json}; @@ -30,6 +31,14 @@ pub const STATIC: CommandSpec = CommandSpec { run: run_static, }; +pub const NONLINEAR: CommandSpec = CommandSpec { + id: "solve.nonlinear", + label: "Solve Nonlinear Static Step", + params: r#"{"step": "Inflate", "max_dofs": 6000000, "max_factor_entries": 400000000, "max_newton_iterations": 100000, "threads": 8} — solves a nonlinear static step (the only one when `step` is left out) in load increments with Newton–Raphson iterations and keeps the result in the session: the final displacements, reactions, Cauchy stresses and logarithmic strains, every increment's load factor, iterations and residuals (force, displacement and energy, relative to their tolerances' scales), and a displacement frame per converged increment; the caps may only be lowered (`max_newton_iterations` over the whole step); `threads` as for solve.static"#, + kind: CommandKind::Model, + run: run_nonlinear, +}; + pub const HEAT: CommandSpec = CommandSpec { id: "solve.heat", label: "Solve Heat Transfer Step", @@ -117,6 +126,7 @@ fn limits(p: &Value) -> Result { max_dofs: cap(p, "max_dofs", Limits::MAX_DOFS)?, max_factor_entries: cap(p, "max_factor_entries", Limits::MAX_FACTOR_ENTRIES)?, max_lanczos_vectors: cap(p, "max_iterations", Limits::MAX_LANCZOS_VECTORS)?, + max_newton_iterations: cap(p, "max_newton_iterations", Limits::MAX_NEWTON_ITERATIONS)?, }) } @@ -194,8 +204,34 @@ fn run_static(s: &mut Session, p: &Value) -> Result { Ok(summary) } -fn run_nonlinear(_: &mut Session, _: &Value) -> Result { - Err(EngineError::Params("not supported yet: nonlinear static steps (the nonlinear solver arrives with solve.nonlinear)".into())) +fn run_nonlinear(s: &mut Session, p: &Value) -> Result { + let limits = limits(p)?; + let step = opt::(p, "step")?; + let model = &s.model; + let r = parallel::with_threads(threads(s, p)?, || solve_nonlinear(model, step.as_deref(), &limits))?; + let history = r.nonlinear.as_ref(); + let summary = json!({ + "step": r.step, + "dofs": r.stats.dofs, + "equations": r.stats.equations, + "factor_entries": r.stats.factor_entries, + "increments": history.map(|h| &h.increments), + "iterations": history.map(|h| h.iterations), + "frames": history.map(|h| h.load_factors()), + "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))), + "applied_total": r.applied_total, + "reaction_total": r.reaction_total, + "reference_points": r.reference_points.iter().map(|p| json!({ + "name": p.name, + "label": p.label, + "displacement": p.displacement, + "reaction": p.reaction, + })).collect::>(), + "warnings": r.warnings, + }); + s.store_result(StepResult::Implicit(ImplicitResult::Static(r))); + Ok(summary) } fn run_heat(s: &mut Session, p: &Value) -> Result { @@ -1359,4 +1395,51 @@ pub(crate) mod tests { assert_eq!(fields, vec![crate::results::Field::Temperature, crate::results::Field::HeatFlux]); assert!(crate::results::displacements(heat, crate::results::Pick::default()).is_err()); } + + /// `solve.nonlinear`: a cantilever bent far by a dead tip load, with its increment history, + /// a frame per converged increment shown by `results.frame`, `solve.run` dispatching to it, + /// and the same bits on 1, 2 and 7 threads. + #[test] + fn nonlinear_steps_solve_with_their_history_and_frames() { + let solve = |threads: usize| -> (Value, u64, Session) { + let mut s = Session::default(); + cantilever(&mut s); + let script = json!([ + {"command": "model.section.edit", "params": {"name": "bar", "set": "all", "material": "Mat", "kind": {"type": "solid"}, "formulation": {}}}, + {"command": "model.step.create", "params": {"name": "Large", "procedure": {"type": "nonlinear_static", "initial_increment": 0.25}, + "boundary_conditions": ["clamp"], "loads": ["push"]}}, + {"command": "model.load.edit", "params": {"name": "push", "type": "force", "target": {"set": "tip"}, "vector": [0, 0, -15]}}, + ]); + s.run_script(&script).unwrap(); + let r = s.execute("solve.nonlinear", &json!({"step": "Large", "threads": threads})).unwrap(); + let hash = result_hash(&s, "Large"); + (r, hash, s) + }; + let (r, one, mut s) = solve(1); + for threads in [2, 7] { + assert_eq!(solve(threads).1, one, "{threads} threads"); + } + let increments = r["increments"].as_array().unwrap(); + assert!(increments.iter().any(|i| i["converged"] == true) && increments.last().unwrap()["load_factor"] == 1.0, "{r}"); + assert!(increments[0]["iterations"].as_array().is_some_and(|it| !it.is_empty() && it[0]["force"].is_number()), "{r}"); + let frames = r["frames"].as_array().unwrap().len(); + assert_eq!(frames, 1 + increments.iter().filter(|i| i["converged"] == true).count()); + // Far from linear: P L² / (E I) = 15·100/(1000/12) = 18, the tip swings well down and in. + let tip = r["max_displacement"]["value"].as_f64().unwrap(); + assert!(tip > 8.0 && tip < 12.0, "{r}"); + let view = s.execute("results.frame", &json!({"step": "Large", "frame": 1})).unwrap(); + assert_eq!(view["view"]["frame"], 1, "{view}"); + assert!(view["load_factor"].as_f64().is_some_and(|f| f > 0.0 && f < 1.0), "{view}"); + assert!(s.execute("results.frame", &json!({"step": "Large", "frame": frames})).is_err()); + let probe = s.execute("solve.probe", &json!({"step": "Large", "nodes": [11]})).unwrap(); + assert!(probe["nodes"][0]["displacement"][2].as_f64().is_some_and(|w| w < -5.0), "{probe}"); + // solve.run solves every step with its own solver; the hostile caps are refused. + assert!(s.execute("solve.nonlinear", &json!({"step": "Large", "max_newton_iterations": 0})).is_err()); + assert!(s.execute("solve.nonlinear", &json!({"step": "Large", "max_newton_iterations": 1e300})).is_err()); + let e = s.execute("solve.nonlinear", &json!({"step": "Large", "max_newton_iterations": 2})).unwrap_err().to_string(); + assert!(e.contains("Newton iterations"), "{e}"); + assert!(s.execute("solve.nonlinear", &json!({"step": "Bend"})).is_err(), "a linear step"); + let all = s.execute("solve.run", &json!({})).unwrap(); + assert_eq!(all["steps"].as_array().map(Vec::len), Some(3), "{all}"); + } } diff --git a/crates/engine/src/results.rs b/crates/engine/src/results.rs index ba9b328..22baa0b 100644 --- a/crates/engine/src/results.rs +++ b/crates/engine/src/results.rs @@ -268,7 +268,12 @@ impl ResultsView { /// of the picked mode or the picked frame's displacement, in [`StepResult::nodes`] order. pub fn displacements(r: &StepResult, pick: Pick) -> Result<&[[f64; 3]]> { match r.kind() { - Kind::Static(s) => Ok(&s.displacement), + Kind::Static(s) => match (pick.frame, &s.nonlinear) { + (Some(k), Some(h)) => h.frames.get(k).map(|f| f.displacement.as_slice()).ok_or_else(|| { + EngineError::Params(format!("step `{}` has frames 0 to {}; there is no frame {k}", s.step, h.frames.len().saturating_sub(1))) + }), + _ => Ok(&s.displacement), + }, Kind::Frequency(f) => pick .mode .checked_sub(1) @@ -418,7 +423,13 @@ impl Session { out["time"] = json!(t); } } - Kind::Static(_) => {} + Kind::Static(s) => { + if let Some(h) = &s.nonlinear { + out["frames"] = json!(h.load_factors()); + let shown = v.frame.and_then(|k| h.frames.get(k)).or(h.frames.last()); + out["load_factor"] = json!(shown.map(|f| f.load_factor)); + } + } } out } diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 1692d81..1cf6337 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -106,6 +106,9 @@ pub struct StaticResult { pub warnings: Vec, /// The temperatures of a coupled temperature–displacement step, per node. pub temperature: Option>, + /// The increments of a nonlinear static step (its stresses are Cauchy stresses and its + /// strains logarithmic strains). + pub nonlinear: Option, } /// One natural mode. @@ -167,10 +170,10 @@ struct Assembled { } /// Elements whose matrices are computed together before they are added. -const ROUND: usize = 4096; +pub(crate) const ROUND: usize = 4096; /// Elements per parallel task. -const GRAIN: usize = 16; +pub(crate) const GRAIN: usize = 16; /// Assemble the stiffness (and the mass when `mass` is set) on the mesh's node graph. fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option) -> Result { @@ -240,13 +243,13 @@ fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option) -> Result std::result::Result { +pub(crate) fn factor_raw(a: &Csr, limits: &Limits) -> std::result::Result { Factor::new(&UpperColumns::from_symmetric(a), limits.max_factor_entries) } /// A factorization failure in the user's terms: the node (or reference point) and direction of /// a singular pivot. -fn factor_error(e: FactorError, mesh: &Mesh, reduction: &Reduction) -> Error { +pub(crate) fn factor_error(e: FactorError, mesh: &Mesh, reduction: &Reduction) -> Error { match e { FactorError::TooLarge { entries, limit } => { Error::Limit(format!("the factorization needs {entries} entries; the limit is {limit} (refine less or raise the limit)")) @@ -270,7 +273,7 @@ fn directions(full: &[f64], layout: &DofLayout, slot: usize) -> /// A full vector split into the mesh nodes' translations and rotations (zero where a node has /// none), per node, and the reference points' translations and rotations, per point. -fn split(full: &[f64], mesh: &Mesh) -> (Vec<[f64; 3]>, Vec<[f64; 3]>, Vec<[f64; 6]>) { +pub(crate) fn split(full: &[f64], mesh: &Mesh) -> (Vec<[f64; 3]>, Vec<[f64; 3]>, Vec<[f64; 6]>) { let l = &mesh.layout; let (translation, rotation) = (0..l.nodes()) .map(|n| { @@ -283,7 +286,7 @@ fn split(full: &[f64], mesh: &Mesh) -> (Vec<[f64; 3]>, Vec<[f64; 3]>, Vec<[f64; /// The matrix pattern of the mesh (node blocks of its elements, trusses and discrete elements) /// and the order of its nodes for the factorization. -fn pattern(mesh: &Mesh) -> (NodeGraph, Vec) { +pub(crate) fn pattern(mesh: &Mesh) -> (NodeGraph, Vec) { let graph = NodeGraph::new(mesh.block_count(), mesh.blocks()); // The order covers the mesh nodes; reference point equations follow them. let n = mesh.labels.len(); @@ -297,7 +300,7 @@ fn pattern(mesh: &Mesh) -> (NodeGraph, Vec) { } /// Warnings for what a step without velocities ignores. -fn ignored_dashpots(mesh: &Mesh, warnings: &mut Vec) { +pub(crate) fn ignored_dashpots(mesh: &Mesh, warnings: &mut Vec) { let names = mesh.discrete.dashpots(); if !names.is_empty() { warnings.push(format!( @@ -309,20 +312,20 @@ fn ignored_dashpots(mesh: &Mesh, warnings: &mut Vec) { /// Constraint rows, prescribed values and their elimination for a step. /// The model, its resolved view, the mesh and the step being solved. -struct Setup<'a> { - model: &'a Model, - resolved: &'a Resolved, - mesh: &'a Mesh, - step: &'a Step, +pub(crate) struct Setup<'a> { + pub(crate) model: &'a Model, + pub(crate) resolved: &'a Resolved, + pub(crate) mesh: &'a Mesh, + pub(crate) step: &'a Step, } -fn reduction(s: &Setup, values: bool, order: &[usize], budget: Budget, warnings: &mut Vec) -> Result { +pub(crate) fn reduction(s: &Setup, values: bool, order: &[usize], budget: Budget, warnings: &mut Vec) -> Result { let rows = coupling_rows(s.model, s.resolved, s.mesh)?; Reduction::new(s.mesh, order, &rows, |active| prescribed(s.model, s.resolved, s.mesh, s.step, values, active, warnings), budget) } /// The reduced matrices' memory budget: the factorization's. -fn budget(limits: &Limits, mesh: &Mesh, graph: &NodeGraph, matrices: usize) -> Budget { +pub(crate) fn budget(limits: &Limits, mesh: &Mesh, graph: &NodeGraph, matrices: usize) -> Budget { Budget { max_bytes: limits.max_factor_entries.saturating_mul(8), pattern_entries: graph.entries(&mesh.layout), matrices } } @@ -499,6 +502,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul stats: Stats { dofs: mesh.dofs(), equations: red.neq(), factor_entries: chol.entries(), lanczos_vectors: 0 }, warnings, temperature, + nonlinear: None, }) } diff --git a/crates/solver/src/constraints.rs b/crates/solver/src/constraints.rs index ce49732..504eee4 100644 --- a/crates/solver/src/constraints.rs +++ b/crates/solver/src/constraints.rs @@ -650,6 +650,50 @@ impl Reduction { .collect() } + /// T_fᵀ r: a vector over every degree of freedom (a residual) restricted to the free + /// equations. + pub fn restrict(&self, r: &[f64]) -> Vec { + let mut out = vec![0.0; self.neq()]; + for (dof, v) in r.iter().enumerate() { + for (e, c) in self.free.row(dof) { + if let Some(x) = out.get_mut(e) { + *x += c * v; + } + } + } + out + } + + /// Every degree of freedom from the free equations' values `q` with the prescribed values + /// scaled by `scale` (a load factor): T_f q + scale T_p u_p. + pub fn expand_scaled(&self, q: &[f64], scale: f64) -> Vec { + (0..self.prescribed_motion.len()) + .map(|d| { + let free: f64 = self.free.row(d).map(|(e, c)| c * q.get(e).copied().unwrap_or(0.0)).sum(); + free + scale * self.prescribed_motion.get(d).copied().unwrap_or(0.0) + }) + .collect() + } + + /// The largest prescribed motion, |T_p u_p|∞. + pub fn largest_prescribed(&self) -> f64 { + self.prescribed_motion.iter().fold(0.0, |m, v| m.max(v.abs())) + } + + /// T_pᵀ g for a vector g over every degree of freedom (the reactions of the out-of-balance + /// forces g = f_int − f_ext), zero where nothing is prescribed. + pub fn reactions_of(&self, g: &[f64]) -> Vec { + let mut out = vec![0.0; self.prescribed_motion.len()]; + for (dof, v) in g.iter().enumerate() { + for (p, c) in self.fixed.row(dof) { + if let Some(x) = out.get_mut(p) { + *x += c * v; + } + } + } + out + } + /// The reactions at the prescribed degrees of freedom, T_pᵀ (K u − f), over every degree of /// freedom (zero where nothing is prescribed). pub fn reactions(&self, k: &Csr, u: &[f64], f: &[f64]) -> Vec { diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index 41e9737..b7aca2a 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -31,6 +31,7 @@ pub mod constraints; pub mod discrete; pub mod eigen; pub mod heat; +pub mod nonlinear; pub mod ordering; pub mod sparse; pub mod system; @@ -40,6 +41,7 @@ pub use analysis::{ solve_static, }; pub use heat::{HeatResult, solve_heat}; +pub use nonlinear::solve_nonlinear; #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { @@ -82,16 +84,24 @@ pub struct Limits { pub max_factor_entries: usize, /// Lanczos vectors (iterations × block size). pub max_lanczos_vectors: usize, + /// Newton iterations of a nonlinear step, over all its increments. + pub max_newton_iterations: usize, } impl Limits { pub const MAX_DOFS: usize = 6_000_000; pub const MAX_FACTOR_ENTRIES: usize = 400_000_000; pub const MAX_LANCZOS_VECTORS: usize = 3_000; + pub const MAX_NEWTON_ITERATIONS: usize = 100_000; } impl Default for Limits { fn default() -> Limits { - Limits { max_dofs: Limits::MAX_DOFS, max_factor_entries: Limits::MAX_FACTOR_ENTRIES, max_lanczos_vectors: Limits::MAX_LANCZOS_VECTORS } + Limits { + max_dofs: Limits::MAX_DOFS, + max_factor_entries: Limits::MAX_FACTOR_ENTRIES, + max_lanczos_vectors: Limits::MAX_LANCZOS_VECTORS, + max_newton_iterations: Limits::MAX_NEWTON_ITERATIONS, + } } } diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs new file mode 100644 index 0000000..7cad70a --- /dev/null +++ b/crates/solver/src/nonlinear.rs @@ -0,0 +1,897 @@ +//! Nonlinear static analysis: equilibrium followed in load increments with Newton–Raphson +//! iterations. +//! +//! **Load path.** The step's load factor λ runs from 0 to 1. The dead loads (forces, tractions, +//! gravity, body forces; pressure too without nonlinear geometry) are λ times their values, a +//! follower pressure λ p on the deformed faces, and every prescribed displacement λ times its +//! value (u = T_f q + λ T_p u_p, the constraint elimination of the linear solver, +//! [`crate::constraints`]). +//! +//! **Newton–Raphson** with the full consistent tangent, refactored every iteration (Bathe +//! 2014, §8.4; M. A. Crisfield, *Non-linear Finite Element Analysis of Solids and Structures*, +//! vol. 1, Wiley 1991, ch. 9): K_T δq = r with r = T_fᵀ (λ f_ext(u) − f_int(u)), then, when the +//! step asks for it (it is off by default: on large elastic rotations the full Newton step +//! converges in fewer iterations), a line search on the work of the residual along δq, G(s) = δq · r(q + s δq): s = 1 is taken when +//! |G(1)| ≤ 0.8 |G(0)| (always near the solution, so quadratic convergence is kept), otherwise +//! s is moved to the root of the line through (0, G(0)) and (s, G(s)), within [0.1, 1], at most +//! five times (Crisfield 1991, §9.3). An increment has converged when, at once, +//! - **force**: the largest residual component is at most `force` times the largest force the +//! step has seen (applied, internal, reaction), never less than a round-off floor of 10⁻⁸ of a +//! typical diagonal stiffness times the model's size; +//! - **displacement**: the largest correction is at most `displacement` times the largest +//! displacement of the increment; +//! - **energy**: the work of the correction against the residual, s δq · r, is at most `energy` +//! times that of the first iteration (K.-J. Bathe, A. P. Cimento, *Some practical procedures +//! for the solution of nonlinear finite element equations*, CMAME 22:59–85, 1980); +//! +//! or when the force residual is a thousandth of its tolerance (round-off: no correction can +//! improve it). +//! +//! **Automatic incrementation.** Each increment starts from the linear extrapolation of the +//! last converged one (the free motion scaled by the ratio of their sizes), which keeps the +//! inner nodes in step with the prescribed ones; from the last converged state when the +//! extrapolation inverts an element. The first increment is the step's `initial_increment`; two +//! increments in a row that converge in at most four iterations let the next grow by half (up +//! to `max_increment`), one that needs more than eight shrinks it by a quarter. An increment is +//! cut back to a quarter of its size when it fails: an element inverts, the residual is not +//! finite, the residual grows in two successive iterations, the tangent is not positive definite +//! (a limit point or an unstable material), or it does not converge in `max_iterations`. The +//! step fails cleanly, naming the load factor and the last reason, when an increment falls below +//! `min_increment`, when it has taken `max_increments` increments, or when the Newton iterations +//! of the whole step reach [`crate::Limits::max_newton_iterations`]. +//! +//! **Follower pressure** and the Cholesky factorization: the load stiffness of a pressure on the +//! deformed surface is not symmetric (Schweizerhof & Ramm 1984); its skew part is left out of the +//! tangent (the symmetric part is kept), so the factorization stays the supernodal Cholesky of the +//! linear solver. The residual is always exact, so the converged solution is unchanged; only the +//! rate of convergence can drop where the skew part matters. Over a closed surface, or one whose +//! boundary is held, the skew parts of neighbouring faces cancel and nothing is left out; the +//! result reports the skew part's size relative to the load stiffness when it is not small. +//! +//! **Elements.** Solids (tet4, tet10, wedge6, hex8, hex20 full and reduced, hybrid hex8, hex20 +//! and tet10) and 2D plane strain and axisymmetric elements, total Lagrangian +//! ([`fieldcraft_fem::finite`]); without nonlinear geometry the elements are the linear ones +//! (every material is elastic today, so such a step is the linear solution reached by Newton). +//! Springs act along their fixed directions. Trusses, shells and beams (large rotations), plane +//! stress, incompatible modes, temperature and centrifugal loads are refused by name. +//! Couplings stay the linear (small-rotation) constraints of the linear solver, with a warning. +//! +//! Element forces and tangents are computed in parallel and added in label order, so results are +//! the same bit for bit on any number of threads. + +use fieldcraft_fem::constitutive::MaterialModel; +use fieldcraft_fem::dense::Matrix; +use fieldcraft_fem::element::{Extrapolation, extrapolation}; +use fieldcraft_fem::quadrature::Rule; +use fieldcraft_fem::{DofKind, DofSet, ElementType, PointStress, face, von_mises}; +use fieldcraft_model::{Extent, Load, Model, NonlinearStatic, Procedure, Step}; +use fieldcraft_parallel as parallel; +use serde::Serialize; + +use crate::analysis::{ + ElementStress, GRAIN, PointResult, ROUND, Setup, StaticResult, Stats, budget, factor_error, factor_raw, ignored_dashpots, pattern, reduction, + split, +}; +use crate::cholesky::FactorError; +use crate::constraints::Reduction; +use crate::sparse::{Csr, add_block}; +use crate::system::{Mesh, inelastic_warnings, loads, translations}; +use crate::{Error, Limits, Result}; + +/// One Newton iteration of an increment. +#[derive(Debug, Clone, Copy, PartialEq, Serialize)] +pub struct Iteration { + /// Largest residual component over the force scale. + pub force: f64, + /// Largest correction over the largest displacement of the increment. + pub displacement: f64, + /// Work of the correction against the residual over the first iteration's. + pub energy: f64, + /// The line search's step along the correction (1: the full Newton step). + pub line_search: f64, +} + +/// One increment (converged, or cut back). +#[derive(Debug, Clone, PartialEq, Serialize)] +pub struct Increment { + /// From 1, attempts included. + pub number: usize, + /// The load factor reached (when converged) or attempted. + pub load_factor: f64, + /// Its size, as a fraction of the step. + pub size: f64, + pub iterations: Vec, + pub converged: bool, + /// Why it was cut back. + #[serde(skip_serializing_if = "Option::is_none")] + pub cutback: Option, +} + +/// The displacements at the end of a converged increment. +#[derive(Debug, Clone, PartialEq)] +pub struct Frame { + pub load_factor: f64, + /// Per node, in [`StaticResult::nodes`] order. + pub displacement: Vec<[f64; 3]>, +} + +/// The increments of a nonlinear static step and the displacements after each converged one. +#[derive(Debug, Clone, PartialEq, Default)] +pub struct History { + pub increments: Vec, + /// The undeformed state (load factor 0), then one per converged increment while they fit + /// [`MAX_FRAME_VALUES`]. + pub frames: Vec, + /// Newton iterations over the step. + pub iterations: usize, +} + +impl History { + /// The load factor of every frame. + pub fn load_factors(&self) -> Vec { + self.frames.iter().map(|f| f.load_factor).collect() + } +} + +/// Most displacement values the frames of one step keep (3 per node per frame). +pub const MAX_FRAME_VALUES: usize = 30_000_000; + +/// 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. +const LINE_SEARCH_TRIALS: usize = 5; +const LINE_SEARCH_MIN: f64 = 0.1; +/// Cutback factor of a failed increment, growth after fast ones and shrinking after slow ones. +const CUTBACK: f64 = 0.25; +const GROWTH: f64 = 1.5; +const SHRINK: f64 = 0.75; +const FAST_ITERATIONS: usize = 4; +const SLOW_ITERATIONS: usize = 8; +/// A force residual this fraction of its tolerance is round-off. +const ROUND_OFF: f64 = 1e-3; + +/// A pressure on element faces (or 2D edges) that follows the deformation. +struct Follower { + element: u32, + /// Layout slots of the face's nodes, and their reference positions. + nodes: Vec, + xyz: Vec<[f64; 3]>, + magnitude: f64, + /// The out-of-plane extent of a 2D edge; `None` for a solid face. + extent: Option, + /// The translations of its nodes (x, y, z, or x, y per node). + dofs: Vec, +} + +/// Everything an evaluation of the out-of-balance forces needs. +struct Problem<'a> { + mesh: &'a Mesh, + red: &'a Reduction, + /// Dead loads at λ = 1. + dead: Vec, + followers: Vec, + /// Springs (linear along their fixed directions), on the pattern. + springs: Option, + /// Large-deformation kinematics. + finite: bool, + /// Each element's global degrees of freedom, in its matrices' order, and the kinds per node. + element_dofs: Vec>, + element_kinds: Vec, + /// Each element's material model. + materials: Vec, + /// The tangent's pattern, zero. + pattern: Csr, +} + +/// The internal and external forces at a state, and the tangent when asked for. +struct Evaluation { + internal: Vec, + external: Vec, + tangent: Option, + /// Every element's material history at this state (committed when the increment + /// converges). + history: Vec>, +} + +/// Why an attempt at an increment failed, when that calls for a cutback rather than an error. +enum Failure { + Cutback(String), + Fatal(Error), +} + +impl From for Failure { + fn from(e: Error) -> Failure { + match e { + // An inverted element or a material that cannot take the deformation: a smaller + // increment may avoid it. + Error::Element { label, source: fieldcraft_fem::Error::Degenerate(m) } => Failure::Cutback(format!("element {label}: {m}")), + Error::NotConverged(m) => Failure::Cutback(m), + e => Failure::Fatal(e), + } + } +} + +fn norm_inf(v: &[f64]) -> f64 { + v.iter().fold(0.0, |m, x| m.max(x.abs())) +} + +fn dot(a: &[f64], b: &[f64]) -> f64 { + a.iter().zip(b).map(|(x, y)| x * y).sum() +} + +/// ½(K + sign Kᵀ) of a row-major m × m matrix: its symmetric (sign 1) or skew (−1) part. +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() +} + +/// 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(); + red.restrict(&r) +} + +impl Problem<'_> { + /// Nodal displacements of an element (3 per node for solids, 2 for 2D elements). + fn element_vector(&self, e: usize, u: &[f64]) -> Vec { + self.element_dofs.get(e).map(|dofs| dofs.iter().map(|d| u.get(*d).copied().unwrap_or(0.0)).collect()).unwrap_or_default() + } + + /// Internal forces and the tangent of element `e` at displacements `u` (all degrees of + /// freedom). + fn element(&self, e: usize, u: &[f64], tangent: bool, committed: &[f64]) -> Result<(Vec, Option, Vec)> { + let el = self.mesh.elements.get(e).ok_or_else(|| Error::Internal("element index".into()))?; + let at = |err| Error::element(el.label, err); + let ue = self.element_vector(e, u); + let kernel = el.kernel()?; + if !self.finite { + let k = kernel.stiffness(&el.d).map_err(at)?; + let f = k.mul_vec(&ue); + return Ok((f, tangent.then_some(k), committed.to_vec())); + } + let material = self.materials.get(e).ok_or_else(|| Error::Internal("element material".into()))?; + let mut updated = committed.to_vec(); + let r = kernel.finite_strain(material, &ue, committed, &mut updated, tangent).map_err(at)?; + Ok((r.force, r.tangent, updated)) + } + + /// The forces of a follower pressure at λ (on [`Follower::dofs`]) and its load stiffness + /// (full, not symmetric). + fn follower(&self, f: &Follower, u: &[f64], lambda: f64) -> Result<(Vec, Vec)> { + let at = |err| Error::element(f.element, err); + let ue: Vec = f.dofs.iter().map(|d| u.get(*d).copied().unwrap_or(0.0)).collect(); + let p = lambda * f.magnitude; + let (force, k) = match f.extent { + None => { + let (force, k) = face::follower_pressure(&f.xyz, ue.as_chunks::<3>().0, p).map_err(at)?; + (force.into_iter().flatten().collect(), k) + } + Some(extent) => { + let (force, k) = face::follower_edge_pressure(&f.xyz, ue.as_chunks::<2>().0, extent, p).map_err(at)?; + (force.into_iter().flatten().collect(), k) + } + }; + Ok((force, k)) + } + + /// The forces at displacements `u` and load factor λ, and the tangent (symmetric part of the + /// load stiffness included) when `tangent`. + fn evaluate(&self, u: &[f64], lambda: f64, tangent: bool, committed: &[Vec]) -> Result { + let layout = &self.mesh.layout; + let n = u.len(); + let mut internal = vec![0.0; n]; + let mut k = tangent.then(|| self.pattern.clone()); + let mut history = Vec::with_capacity(self.mesh.elements.len()); + for round in parallel::blocks(self.mesh.elements.len(), ROUND) { + let start = round.start; + let computed = parallel::map_grained(round.len(), GRAIN, |i| { + self.element(start + i, u, tangent, committed.get(start + i).map(Vec::as_slice).unwrap_or_default()) + }); + for (i, c) in computed.into_iter().enumerate() { + let (f, ke, updated) = c?; + history.push(updated); + let e = start + i; + let el = self.mesh.elements.get(e).ok_or_else(|| Error::Internal("element index".into()))?; + for (d, v) in self.element_dofs.get(e).into_iter().flatten().zip(&f) { + if let Some(x) = internal.get_mut(*d) { + *x += v; + } + } + if let (Some(k), Some(ke)) = (k.as_mut(), ke) { + let kinds = self.element_kinds.get(e).copied().ok_or_else(|| Error::Internal("element kinds".into()))?; + add_block(k, layout, &el.nodes, kinds, &ke.data)?; + } + } + } + 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); + for (x, v) in internal.iter_mut().zip(&ku) { + *x += v; + } + // Both are on the mesh's pattern. + if let Some(k) = k.as_mut() { + k.val.iter_mut().zip(&s.val).for_each(|(x, v)| *x += v); + } + } + let mut external: Vec = self.dead.iter().map(|f| lambda * f).collect(); + for f in &self.followers { + let (force, kl) = self.follower(f, u, lambda)?; + for (d, v) in f.dofs.iter().zip(&force) { + if let Some(x) = external.get_mut(*d) { + *x += v; + } + } + if let Some(k) = k.as_mut() { + add_block(k, layout, &f.nodes, self.mesh.translation_kinds(), &half(&kl, f.dofs.len(), 1.0))?; + } + } + if internal.iter().chain(&external).any(|v| !v.is_finite()) { + return Err(Error::NotConverged("the forces are not finite".into())); + } + Ok(Evaluation { internal, external, tangent: k, history }) + } + + /// 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> { + if self.followers.is_empty() { + return Ok(None); + } + let layout = &self.mesh.layout; + let (mut skew, mut whole) = (self.pattern.clone(), self.pattern.clone()); + for f in &self.followers { + let (_, kl) = self.follower(f, u, lambda)?; + add_block(&mut skew, layout, &f.nodes, self.mesh.translation_kinds(), &half(&kl, f.dofs.len(), -1.0))?; + add_block(&mut whole, layout, &f.nodes, self.mesh.translation_kinds(), &kl)?; + } + let mut free = vec![false; skew.n]; + for d in &self.red.eq_dof { + if let Some(x) = free.get_mut(*d) { + *x = true; + } + } + let norm = |a: &Csr| -> f64 { + let mut sum = 0.0; + for r in (0..a.n).filter(|r| free.get(*r).copied().unwrap_or(false)) { + let (cols, vals) = a.row(r); + for (c, v) in cols.iter().zip(vals) { + if free.get(*c).copied().unwrap_or(false) { + sum += v * v; + } + } + } + sum.sqrt() + }; + let total = norm(&whole); + Ok((total > 0.0).then(|| norm(&skew) / total)) + } +} + +/// Mutable state of the Newton loop over the step. +struct Run<'a> { + problem: &'a Problem<'a>, + settings: &'a NonlinearStatic, + limits: &'a Limits, + /// The largest force seen so far (∞-norm), and the round-off floor of the force scale. + force_scale: f64, + force_floor: Option, + length: f64, + iterations: usize, + factor_entries: usize, + /// Every element's material history at the end of the last converged increment + /// (integration points × the material's history length). + committed: Vec>, +} + +/// The outcome of a converged increment. +struct Converged { + q: Vec, + evaluation: Evaluation, +} + +impl Run<'_> { + /// The residual over the free equations, T_fᵀ (f_ext − f_int), and the force scale. + fn residual(&mut self, e: &Evaluation) -> Vec { + self.force_scale = self.force_scale.max(norm_inf(&e.external)).max(norm_inf(&e.internal)); + free_residual(self.problem.red, e) + } + + fn scale(&self) -> f64 { + self.force_scale.max(self.force_floor.unwrap_or(0.0)).max(f64::MIN_POSITIVE) + } + + fn state(&self, q: &[f64], lambda: f64, tangent: bool) -> std::result::Result { + Ok(self.problem.evaluate(&self.problem.red.expand_scaled(q, lambda), lambda, tangent, &self.committed)?) + } + + /// Newton iterations from the converged `q0` (at the previous load factor) to load factor + /// `lambda`. + fn increment( + &mut self, + q0: &[f64], + predictor: Option<&[f64]>, + lambda: f64, + size: f64, + iterations: &mut Vec, + ) -> std::result::Result { + let red = self.problem.red; + let tol = self.settings.tolerances; + // 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 { + Some(Ok(s)) => s, + Some(Err(Failure::Fatal(e))) => return Err(Failure::Fatal(e)), + Some(Err(Failure::Cutback(_))) | None => (q0.to_vec(), self.state(q0, lambda, true)?), + }; + let mut r = self.residual(&eval); + let prescribed = size * red.largest_prescribed(); + let mut first_energy: Option = None; + let mut history: Vec = Vec::new(); + let max_iterations = usize::try_from(self.settings.max_iterations).unwrap_or(usize::MAX); + for _ in 0..max_iterations { + // The round-off rule: nothing left to correct. + if norm_inf(&r) <= ROUND_OFF * tol.force * self.scale() && self.force_floor.is_some() { + return Ok(Converged { q, evaluation: eval }); + } + 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 k = eval.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; + let [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 { + 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); + } + drop(k); + let chol = match factor_raw(&kff, self.limits) { + Ok(c) => c, + Err(FactorError::NotPositiveDefinite { .. }) if !first => { + return Err(Failure::Cutback( + "the tangent stiffness is not positive definite (a limit point, a bifurcation or an unstable material)".into(), + )); + } + Err(e) => return Err(Failure::Fatal(factor_error(e, self.problem.mesh, red))), + }; + self.factor_entries = chol.entries(); + drop(kff); + let mut dq = r.clone(); + chol.solve(&mut dq)?; + if dq.iter().any(|v| !v.is_finite()) { + 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)?; + for (x, d) in q.iter_mut().zip(&dq) { + *x += s * d; + } + eval = next; + let r_old = std::mem::replace(&mut r, self.residual(&eval)); + 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())); + let reference = moved.max(prescribed).max(1e-12 * self.length); + let it = Iteration { + force: norm_inf(&r) / self.scale(), + displacement: s * norm_inf(&dq) / reference, + energy: if first > 0.0 { energy / first } else { 0.0 }, + line_search: s, + }; + iterations.push(it); + if !it.force.is_finite() { + return Err(Failure::Cutback("the residual is not finite".into())); + } + if (it.force <= tol.force && it.displacement <= tol.displacement && it.energy <= tol.energy) || it.force <= ROUND_OFF * tol.force { + return Ok(Converged { q, evaluation: eval }); + } + history.push(it.force); + if let [.., a, b, c] = history[..] + && c > b + && b > a + { + return Err(Failure::Cutback(format!("the residual grows ({a:.2e}, {b:.2e}, {c:.2e} of the force scale)"))); + } + } + Err(Failure::Cutback(format!("no convergence in {} iterations", self.settings.max_iterations))) + } + + /// 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> { + let at = |s: f64| -> Vec { q.iter().zip(dq).map(|(a, b)| a + s * b).collect() }; + // The full step, with its tangent (taken in most iterations). + let full = self.state(&at(1.0), lambda, true); + if !self.settings.line_search { + return full.map(|e| (1.0, e)); + } + let red = self.problem.red; + let work = |e: &Evaluation| dot(dq, &free_residual(red, e)); + let mut s = 1.0; + let mut g = match &full { + Ok(e) => { + let g1 = work(e); + if g1.abs() <= LINE_SEARCH_ACCEPT * g0.abs() { + return full.map(|e| (1.0, e)); + } + Some(g1) + } + Err(Failure::Cutback(_)) => None, + Err(Failure::Fatal(_)) => return full.map(|e| (1.0, e)), + }; + let mut best: Option<(f64, f64)> = g.map(|g| (1.0, g.abs())); + for _ in 0..LINE_SEARCH_TRIALS { + // Root of the line through (0, g0) and (s, g), or half the step after a failure. + let next = match g { + Some(gs) if (g0 - gs).abs() > 0.0 => (s * g0 / (g0 - gs)).clamp(LINE_SEARCH_MIN, 1.0), + Some(_) => break, + None => (0.5 * s).max(LINE_SEARCH_MIN), + }; + if next == s { + break; + } + s = next; + g = match self.state(&at(s), lambda, false) { + Ok(e) => Some(work(&e)), + Err(Failure::Cutback(_)) => None, + Err(f) => return Err(f), + }; + if let Some(gs) = g { + if best.is_none_or(|(_, b)| gs.abs() < b) { + best = Some((s, gs.abs())); + } + if gs.abs() <= LINE_SEARCH_ACCEPT * g0.abs() { + break; + } + } + } + let Some((s, _)) = best else { + return Err(Failure::Cutback("every step along the Newton correction inverts an element".into())); + }; + if s == 1.0 { + return full.map(|e| (1.0, e)); + } + self.state(&at(s), lambda, true).map(|e| (s, e)) + } +} + +/// Refuse what nonlinear steps cannot take yet, and warn about what they simplify. +fn check_supported(model: &Model, mesh: &Mesh, step: &Step, settings: &NonlinearStatic, warnings: &mut Vec) -> Result<()> { + let at = |what: String| Error::NotSupported(format!("{what} in nonlinear steps (step `{}`)", step.name)); + if let Some(t) = mesh.trusses.first() { + return Err(at(format!("trusses (element {})", t.label))); + } + if let Some(s) = mesh.structural.first() { + return Err(at(format!("shells and beams, whose large rotations are not implemented yet (element {})", s.label))); + } + for e in &mesh.elements { + if !settings.nonlinear_geometry && e.hyperelastic.is_some() { + return Err(at(format!("a hyperelastic material without nonlinear geometry (element {})", e.label))); + } + } + for name in &step.loads { + match model.load(name) { + Some(Load::Temperature { .. }) => return Err(at(format!("temperature loads (load `{name}`): thermal strain"))), + Some(Load::Centrifugal { .. }) => return Err(at(format!("centrifugal loads (load `{name}`), which depend on the deformed position"))), + _ => {} + } + } + if settings.nonlinear_geometry + && let Some(name) = model.couplings().keys().next() + { + warnings.push(format!( + "coupling `{name}` and the others are linear constraints (small rotations of their reference points) in this nonlinear step" + )); + } + Ok(()) +} + +/// The follower pressures of a step (with nonlinear geometry), and the step without them for +/// the dead loads. +fn followers(model: &Model, resolved: &fieldcraft_model::Resolved, mesh: &Mesh, step: &Step) -> Result<(Vec, Step)> { + let mut out = Vec::new(); + let mut dead = step.clone(); + for name in &step.loads { + let Some(Load::Pressure { set, magnitude }) = model.load(name) else { continue }; + dead.loads.remove(name); + let slots = |labels: Vec| -> Result<(Vec, Vec<[f64; 3]>, Vec)> { + let idx: Vec = labels + .iter() + .map(|l| mesh.node(*l).ok_or_else(|| Error::Invalid(format!("load `{name}`: node {l} belongs to no solid element")))) + .collect::>()?; + let xyz = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); + let mut dofs = Vec::with_capacity(3 * idx.len()); + for n in &idx { + dofs.extend(translations(mesh, *n)?); + } + Ok((idx, xyz, dofs)) + }; + if let Some(theory) = mesh.theory { + if resolved.sets.get(set).is_some_and(|s| !s.faces.is_empty()) { + return Err(Error::Invalid(format!("load `{name}`: a 2D model's pressures act on element edges, and set `{set}` has element faces"))); + } + for g in resolved.planar_edges(set) { + let labels = + model.edge_nodes(g).ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} edge {} is missing", g.element, g.edge)))?; + let (nodes, xyz, dofs) = slots(labels)?; + out.push(Follower { element: g.element, nodes, xyz, magnitude: *magnitude, extent: Some(theory.extent()), dofs }); + } + } else { + for fc in resolved.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default() { + let labels = model + .face_nodes(*fc) + .ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} face {} is missing", fc.element, fc.face)))?; + let (nodes, xyz, dofs) = slots(labels)?; + out.push(Follower { element: fc.element, nodes, xyz, magnitude: *magnitude, extent: None, dofs }); + } + } + } + Ok((out, dead)) +} + +/// Solve a nonlinear static step: the step named `step`, or the model's only nonlinear static +/// step. +pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Result { + let resolved = model.resolve()?; + let step = model.select_step(step, |p| matches!(p, Procedure::NonlinearStatic(_)), "nonlinear static")?; + let Procedure::NonlinearStatic(settings) = &step.procedure else { + return Err(Error::Internal("nonlinear static step".into())); + }; + settings.validate()?; + let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; + let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).chain(mesh.warnings.iter().cloned()).collect(); + warnings.extend(inelastic_warnings(model)); + check_supported(model, &mesh, step, settings, &mut warnings)?; + let (graph, order) = pattern(&mesh); + ignored_dashpots(&mesh, &mut warnings); + let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; + let red = reduction(&setup, true, &order, budget(limits, &mesh, &graph, 1), &mut warnings)?; + if red.neq() == 0 { + return Err(Error::Invalid("every degree of freedom is prescribed; there is nothing to solve".into())); + } + let (follower_loads, dead_step) = + if settings.nonlinear_geometry { followers(model, &resolved, &mesh, step)? } else { (Vec::new(), step.clone()) }; + let dead = loads(model, &resolved, &mesh, &dead_step, &red.active)?; + let springs = if mesh.discrete.springs.iter().any(|s| !s.dashpot) { + let mut s = graph.zero_matrix(&mesh.layout); + mesh.discrete.add_stiffness(&mut s)?; + Some(s) + } else { + None + }; + let element_kinds = mesh.elements.iter().map(|e| Ok(e.kernel()?.dofs().node)).collect::>>()?; + let element_dofs = mesh + .elements + .iter() + .zip(&element_kinds) + .map(|(e, kinds)| { + mesh.layout.element_dofs(&e.nodes, *kinds).ok_or_else(|| Error::Internal(format!("element {} degrees of freedom", e.label))) + }) + .collect::>>()?; + let problem = Problem { + mesh: &mesh, + red: &red, + dead, + followers: follower_loads, + springs, + finite: settings.nonlinear_geometry, + element_dofs, + element_kinds, + materials: mesh.elements.iter().map(|e| MaterialModel::from_elasticity(&e.d, e.hyperelastic.as_ref())).collect(), + pattern: graph.zero_matrix(&mesh.layout), + }; + let (lo, hi) = mesh.xyz.iter().fold(([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]), |(lo, hi), p| { + (std::array::from_fn(|k| lo[k].min(p[k])), std::array::from_fn(|k| hi[k].max(p[k]))) + }); + let length = (0..3).map(|k| (hi[k] - lo[k]).powi(2)).sum::().sqrt().max(f64::MIN_POSITIVE); + let committed = mesh + .elements + .iter() + .zip(&problem.materials) + .map(|(e, m)| Ok(vec![0.0; e.kernel()?.rule().points().len() * m.history_len()])) + .collect::>>()?; + // The undeformed state first: an element that is bad before anything moves is an error in + // the model, not a reason to cut increments back. + problem.evaluate(&vec![0.0; mesh.dofs()], 0.0, false, &committed)?; + let mut run = + Run { problem: &problem, settings, limits, force_scale: 0.0, force_floor: None, length, iterations: 0, factor_entries: 0, committed }; + + let mut history = History::default(); + let (mut frame_values, mut frames_full) = (0usize, false); + let values = 3 * mesh.labels.len(); + let mut keep_frame = |history: &mut History, lambda: f64, u: &[f64], warnings: &mut Vec| { + if frames_full { + return; + } + if frame_values.saturating_add(values) > MAX_FRAME_VALUES { + warnings.push(format!("the displacement frames stop at load factor {lambda}: the step keeps at most {MAX_FRAME_VALUES} values")); + frames_full = true; + return; + } + frame_values += values; + 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()]; + let mut lambda = 0.0f64; + let mut size = settings.initial_increment; + let mut fast = 0usize; + let mut last: Option = None; + let mut previous: Option<(Vec, f64)> = None; + let max_increments = usize::try_from(settings.max_increments).unwrap_or(usize::MAX); + while lambda < 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 remaining = 1.0 - lambda; + let (target, this) = if size >= remaining * (1.0 - 1e-9) { (1.0, remaining) } else { (lambda + size, size) }; + let mut iterations = Vec::new(); + let number = history.increments.len() + 1; + // Linear extrapolation of the last converged increment's free motion. + let predictor: Option> = + previous.as_ref().map(|(dq, prev): &(Vec, f64)| q.iter().zip(dq).map(|(a, d)| a + this / prev * d).collect()); + match run.increment(&q, predictor.as_deref(), target, this, &mut iterations) { + Ok(c) => { + previous = Some((c.q.iter().zip(&q).map(|(a, b)| a - b).collect(), this)); + let count = iterations.len(); + history.increments.push(Increment { number, load_factor: target, size: this, iterations, converged: true, cutback: None }); + lambda = target; + 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, + 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}: the increment to {target:.6} failed ({reason}) and a smaller one would be below the minimum increment {}", + step.name, settings.min_increment + ))); + } + } + } + } + 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)?, + }; + if let Some(a) = problem.follower_asymmetry(&u, lambda)? + && a > 1e-6 + { + warnings.push(format!( + "follower pressure: the skew part of its load stiffness ({a:.1e} of it, over the free degrees of freedom) is left out of the tangent; the solution is exact, Newton may converge more slowly" + )); + } + finish(&mesh, step, &red, &problem, u, evaluation, &run.committed, history, run.factor_entries, warnings) +} + +/// The final state as a static result: displacements, reactions, Cauchy stresses and +/// logarithmic strains. +#[allow(clippy::too_many_arguments)] +fn finish( + mesh: &Mesh, + step: &Step, + red: &Reduction, + problem: &Problem, + u: Vec, + e: Evaluation, + committed: &[Vec], + history: History, + factor_entries: usize, + mut warnings: Vec, +) -> Result { + let (displacement, rotation, point_motion) = split(&u, mesh); + let out_of_balance: Vec = e.internal.iter().zip(&e.external).map(|(a, b)| a - b).collect(); + let reactions = red.reactions_of(&out_of_balance); + let (reaction, reaction_moment, point_reaction) = split(&reactions, mesh); + let (mut applied_total, mut reaction_total) = ([0.0; 3], [0.0; 3]); + for (dof, (fv, rv)) in e.external.iter().zip(&reactions).enumerate() { + let Some(d) = mesh.kind(dof).and_then(DofKind::direction) else { continue }; + if let (Some(a), Some(r)) = (applied_total.get_mut(d), reaction_total.get_mut(d)) { + *a += fv; + *r += rv; + } + } + let reference_points = mesh + .points + .iter() + .zip(point_motion.iter().zip(&point_reaction)) + .map(|(p, (d, r))| PointResult { name: p.name.clone(), label: p.label, displacement: *d, reaction: *r }) + .collect(); + let all_points = parallel::map_grained(mesh.elements.len(), GRAIN, |i| -> Result<(Vec, &'static ElementType, Rule)> { + let el = mesh.elements.get(i).ok_or_else(|| Error::Internal("element index".into()))?; + let at = |err| Error::element(el.label, err); + let ue = problem.element_vector(i, &u); + let kernel = el.kernel()?; + let points = if problem.finite { + let material = problem.materials.get(i).ok_or_else(|| Error::Internal("element material".into()))?; + kernel.finite_strain_stresses(material, &ue, committed.get(i).map(Vec::as_slice).unwrap_or_default()) + } else { + let per_node: Vec<[f64; 3]> = el.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); + kernel.stresses_with(&el.d, &per_node, None) + } + .map_err(at)?; + Ok((points, kernel.ty(), kernel.rule())) + }); + let mut extrapolations: Vec<(&'static ElementType, Rule, Extrapolation)> = Vec::new(); + let mut elements = Vec::with_capacity(mesh.elements.len()); + let mut sum = vec![[0.0; 6]; mesh.labels.len()]; + let mut count = vec![0u32; mesh.labels.len()]; + for (el, computed) in mesh.elements.iter().zip(all_points) { + let (points, ty, rule) = computed?; + let values: Vec<[f64; 6]> = points.iter().map(|p| p.stress).collect(); + let index = match extrapolations.iter().position(|(t, r, _)| std::ptr::eq(*t, ty) && *r == rule) { + Some(i) => i, + None => { + extrapolations.push((ty, rule, extrapolation(ty, rule).map_err(|err| Error::element(el.label, err))?)); + extrapolations.len() - 1 + } + }; + let nodal = extrapolations.get(index).map(|(_, _, x)| x.apply(&values)).unwrap_or_default(); + for (n, s) in el.nodes.iter().zip(&nodal) { + if let (Some(acc), Some(c)) = (sum.get_mut(*n), count.get_mut(*n)) { + for k in 0..6 { + acc[k] += s[k]; + } + *c += 1; + } + } + elements.push(ElementStress { label: el.label, kind: el.kind, points, nodal }); + } + let stress: Vec<[f64; 6]> = sum.iter().zip(&count).map(|(s, c)| s.map(|v| if *c > 0 { v / f64::from(*c) } else { 0.0 })).collect(); + let von_mises = stress.iter().map(von_mises).collect(); + let scale = applied_total.iter().chain(&reaction_total).fold(0.0f64, |m, v| m.max(v.abs())); + let ground = mesh.discrete.ground_force(&u, mesh); + let imbalance = (0..3).map(|c| (applied_total[c] + reaction_total[c] + ground[c]).abs()).fold(0.0, f64::max); + if scale > 0.0 && imbalance > 1e-6 * scale { + warnings.push(format!("loads and reactions differ by {imbalance:.3e} (relative {:.1e})", imbalance / scale)); + } + Ok(StaticResult { + step: step.name.clone(), + nodes: mesh.labels.clone(), + displacement, + rotation, + reaction, + reaction_moment, + stress, + von_mises, + elements, + shells: Vec::new(), + beams: Vec::new(), + reference_points, + applied_total, + reaction_total, + stats: Stats { dofs: mesh.dofs(), equations: red.neq(), factor_entries, lanczos_vectors: 0 }, + warnings, + temperature: None, + nonlinear: Some(history), + }) +} diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index b7557cc..103c795 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -8,6 +8,7 @@ //! name): translations along X, Y, Z and rotations about them. use fieldcraft_fem::beam::{Beam, BeamSection}; +use fieldcraft_fem::constitutive::MaterialModel; use fieldcraft_fem::dense::Matrix; use fieldcraft_fem::element::{Formulation, Solid, ThermalStrain}; use fieldcraft_fem::plane::{Plane, Theory}; @@ -17,7 +18,7 @@ use fieldcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, ElementDofs, PointS use fieldcraft_model::nodal::edge_pressure; use std::collections::BTreeSet; -use fieldcraft_model::{Elastic, ElementKind, Integration, Load, Model, Resolved, SectionKind, Step, Target, centrifugal_acceleration}; +use fieldcraft_model::{Elastic, ElementKind, Hyperelastic, Integration, Load, Model, Resolved, SectionKind, Step, Target, centrifugal_acceleration}; use crate::analysis::MassKind; use crate::{Error, Result}; @@ -38,6 +39,8 @@ pub struct Elem { /// The expansion coefficients along X, Y, Z and the reference temperature, when the /// material expands. pub expansion: Option<([f64; 3], f64)>, + /// The large-strain option of the element's material (nonlinear steps). + pub hyperelastic: Option, } /// An element's kernel: a solid or a 2D element. @@ -110,6 +113,30 @@ impl Kernel<'_> { } } + /// Finite-strain internal forces and tangent ([`fieldcraft_fem::finite`]) for the element's + /// displacements `u` (its degrees of freedom in order: 3 per node, or 2 for 2D elements). + pub fn finite_strain( + &self, + m: &MaterialModel, + u: &[f64], + committed: &[f64], + updated: &mut [f64], + tangent: bool, + ) -> fieldcraft_fem::Result { + match self { + Kernel::Solid(s) => s.finite_strain(m, u.as_chunks::<3>().0, committed, updated, tangent), + Kernel::Plane(p) => p.finite_strain(m, u.as_chunks::<2>().0, committed, updated, tangent), + } + } + + /// Cauchy stress and logarithmic strain at the integration points, as [`Self::finite_strain`]. + pub fn finite_strain_stresses(&self, m: &MaterialModel, u: &[f64], committed: &[f64]) -> fieldcraft_fem::Result> { + match self { + Kernel::Solid(s) => s.finite_strain_stresses(m, u.as_chunks::<3>().0, committed), + Kernel::Plane(p) => p.finite_strain_stresses(m, u.as_chunks::<2>().0, committed), + } + } + pub fn rule(&self) -> fieldcraft_fem::quadrature::Rule { match self { Kernel::Solid(s) => s.rule(), @@ -449,7 +476,7 @@ impl Mesh { } .map_err(|e| Error::element(label, e))?; let expansion = material.thermal.and_then(|t| t.expansion.map(|a| (a.axes(), t.reference_temperature))); - elements.push((label, el.kind, el.nodes.clone(), formulation, d, material.density, theory, dofs, expansion)); + elements.push((label, el.kind, el.nodes.clone(), formulation, d, material.density, theory, dofs, expansion, material.hyperelastic)); } } if let (Some(_), Some(name)) = (planar, model.couplings().keys().next()) { @@ -515,10 +542,10 @@ impl Mesh { warnings, discrete: crate::discrete::Discrete::default(), }; - for (label, kind, nodes, formulation, d, density, theory, _, expansion) in elements { + for (label, kind, nodes, formulation, d, density, theory, _, expansion, hyperelastic) in elements { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); - mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density, theory, expansion }); + mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density, theory, expansion, hyperelastic }); } for (label, kind, nodes, young, area, density, expansion) in trusses { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); @@ -648,7 +675,7 @@ fn vector_dofs(mesh: &Mesh, slot: usize, rotation: bool) -> Option> { } /// The translations of mesh node `node`. -fn translations(mesh: &Mesh, node: usize) -> Result> { +pub(crate) fn translations(mesh: &Mesh, node: usize) -> Result> { vector_dofs(mesh, node, false).ok_or_else(|| Error::Internal(format!("mesh node {node} has no translations"))) } diff --git a/crates/solver/tests/nonlinear.rs b/crates/solver/tests/nonlinear.rs new file mode 100644 index 0000000..99a3578 --- /dev/null +++ b/crates/solver/tests/nonlinear.rs @@ -0,0 +1,435 @@ +//! Nonlinear static analysis against closed forms: uniaxial stretch of neo-Hookean and +//! St Venant–Kirchhoff blocks (3D, plane strain, axisymmetric), the nonlinear patch test, a +//! large rigid rotation, Newton's quadratic convergence, the linear limit, and what nonlinear +//! steps refuse or stop at. Run with `--nocapture` to see the residuals. + +mod common; + +use common::*; +use fieldcraft_model::*; +use fieldcraft_solver::{Limits, StaticResult, solve_nonlinear, solve_static}; + +fn nonlinear() -> Procedure { + Procedure::NonlinearStatic(NonlinearStatic::default()) +} + +fn neo_hookean(g: &mut Grid, e: f64, nu: f64) { + let m = Material { elastic: Some(iso(e, nu)), hyperelastic: Some(Hyperelastic::NeoHookean {}), ..Material::default() }; + g.model.replace_material("Mat", m).unwrap(); +} + +/// Prescribe every displacement component of each node where `f` holds from its position: +/// one set and one boundary condition per node. +fn prescribe(g: &mut Grid, f: impl Fn([f64; 3]) -> bool, u: impl Fn([f64; 3]) -> [f64; 3]) { + for l in g.nodes_where(f) { + let x = g.model.node(l).unwrap().xyz; + let v = u(x); + let name = format!("n{l}"); + g.model.add_set(&name, Set { nodes: vec![l], ..Set::default() }).unwrap(); + let values = [Some(v[0]), Some(v[1]), Some(v[2]), None, None, None]; + g.bc(&name, BoundaryCondition::Displacement { target: Target::Set(name.clone()), values }); + } +} + +fn at(r: &StaticResult, label: u32) -> usize { + r.nodes.binary_search(&label).unwrap() +} + +fn solve(g: &Grid) -> StaticResult { + solve_nonlinear(&g.model, Some("Step"), &Limits::default()).unwrap() +} + +/// The lateral stretch of a neo-Hookean bar stretched by λ in uniaxial stress (σ_yy = σ_zz = 0) +/// and its axial Cauchy stress, by bisection on σ_yy(l) = μ J^(−5/3) (l² − I₁/3) + κ (J − 1). +fn neo_hookean_uniaxial(mu: f64, kappa: f64, lam: f64) -> (f64, f64) { + let stress = |l: f64| { + let (j, i1) = (lam * l * l, lam * lam + 2.0 * l * l); + let c = mu * j.powf(-5.0 / 3.0); + (c * (lam * lam - i1 / 3.0) + kappa * (j - 1.0), c * (l * l - i1 / 3.0) + kappa * (j - 1.0)) + }; + let (mut lo, mut hi) = (0.1, 1.5); + for _ in 0..200 { + let mid = 0.5 * (lo + hi); + if stress(mid).1 > 0.0 { hi = mid } else { lo = mid } + } + let l = 0.5 * (lo + hi); + (l, stress(l).0) +} + +/// A unit block stretched to twice its length (x), free laterally (symmetry planes), in each +/// element kind: the homogeneous closed form at every node, and the logarithmic strain ln 2. +#[test] +fn neo_hookean_block_in_uniaxial_stretch() { + let (e, nu, lam) = (10.0, 0.3, 2.0); + let (mu, kappa) = (e / (2.0 * (1.0 + nu)), e / (3.0 * (1.0 - 2.0 * nu))); + let (l, sigma) = neo_hookean_uniaxial(mu, kappa, lam); + for (kind, formulation, name) in variants().into_iter().filter(|v| v.1 != incompatible()) { + let mut g = grid([2, 2, 2], kind, |x| x); + g.define(iso(e, nu), None, formulation, nonlinear()); + neo_hookean(&mut g, e, nu); + g.symmetry("x0", Symmetry::X, |x| near(x[0], 0.0)); + g.symmetry("y0", Symmetry::Y, |x| near(x[1], 0.0)); + g.symmetry("z0", Symmetry::Z, |x| near(x[2], 0.0)); + g.node_set("x1", |x| near(x[0], 1.0)); + g.bc("pull", BoundaryCondition::Displacement { target: Target::Set("x1".into()), values: [Some(lam - 1.0), None, None, None, None, None] }); + let r = solve(&g); + for (i, label) in r.nodes.iter().enumerate() { + let x = g.model.node(*label).unwrap().xyz; + let want = [(lam - 1.0) * x[0], (l - 1.0) * x[1], (l - 1.0) * x[2]]; + for k in 0..3 { + assert!((r.displacement[i][k] - want[k]).abs() < 1e-8, "{name} node {label}: {:?} vs {want:?}", r.displacement[i]); + } + } + for s in &r.stress { + assert!((s[0] - sigma).abs() < 1e-8 * sigma && s[1].abs() < 1e-8 * sigma && s[3].abs() < 1e-8 * sigma, "{name}: {s:?} vs {sigma}"); + } + let strain = r.elements[0].points[0].strain; + assert!((strain[0] - lam.ln()).abs() < 1e-9 && (strain[1] - l.ln()).abs() < 1e-9, "{name}: {strain:?}"); + // The reactions on the pulled face balance the axial force σ times the deformed area l². + let pulled: f64 = g.nodes_where(|x| near(x[0], 1.0)).iter().map(|n| r.reaction[at(&r, *n)][0]).sum(); + assert!((pulled - sigma * l * l).abs() < 1e-8 * sigma, "{name}: {pulled} vs {}", sigma * l * l); + let h = r.nonlinear.as_ref().unwrap(); + assert!(h.increments.last().is_some_and(|i| i.converged && i.load_factor == 1.0), "{name}: {:?}", h.increments); + assert_eq!(h.frames.first().map(|f| f.load_factor), Some(0.0)); + assert_eq!(h.frames.last().map(|f| f.load_factor), Some(1.0)); + eprintln!("{name}: {} increments, {} iterations", h.increments.len(), h.iterations); + } +} + +/// St Venant–Kirchhoff in uniaxial stress: S₁₁ = E E₁₁ with E₁₁ = (λ² − 1)/2, lateral +/// E₂₂ = −ν E₁₁; the Cauchy stress σ = λ² S₁₁ / J. +#[test] +fn st_venant_kirchhoff_uniaxial_stretch() { + let (e, nu, lam): (f64, f64, f64) = (1000.0, 0.25, 1.3); + let e11 = 0.5 * (lam * lam - 1.0); + let l = (1.0 - 2.0 * nu * e11).sqrt(); + let sigma = lam * lam * e * e11 / (lam * l * l); + let mut g = grid([1, 1, 1], ElementKind::Hex20, |x| x); + g.define(iso(e, nu), None, full(), nonlinear()); + g.symmetry("x0", Symmetry::X, |x| near(x[0], 0.0)); + g.symmetry("y0", Symmetry::Y, |x| near(x[1], 0.0)); + g.symmetry("z0", Symmetry::Z, |x| near(x[2], 0.0)); + g.node_set("x1", |x| near(x[0], 1.0)); + g.bc("pull", BoundaryCondition::Displacement { target: Target::Set("x1".into()), values: [Some(lam - 1.0), None, None, None, None, None] }); + let r = solve(&g); + let corner = g.nodes_where(|x| near(x[0], 1.0) && near(x[1], 1.0) && near(x[2], 1.0))[0]; + assert!((r.displacement[at(&r, corner)][1] - (l - 1.0)).abs() < 1e-9, "{:?}", r.displacement[at(&r, corner)]); + assert!(r.stress.iter().all(|s| (s[0] - sigma).abs() < 1e-8 * sigma), "{:?} vs {sigma}", r.stress[0]); +} + +/// The nonlinear patch test: a distorted mesh whose boundary nodes follow a homogeneous finite +/// deformation x = F X (stretch, shear and rotation) reproduces it at the inner nodes and gives +/// the uniform stress of F, for every element kind. +#[test] +fn nonlinear_patch_test() { + let f = [[1.3, 0.2, -0.1], [-0.15, 0.9, 0.25], [0.05, -0.2, 1.1]]; + let affine = |x: [f64; 3]| -> [f64; 3] { std::array::from_fn(|i| (0..3).map(|j| f[i][j] * x[j]).sum::() - x[i]) }; + let distort = |q: [f64; 3]| { + let bump = (std::f64::consts::PI * q[0]).sin() * (std::f64::consts::PI * q[1]).sin() * (std::f64::consts::PI * q[2]).sin(); + [q[0] + 0.12 * bump, q[1] - 0.08 * bump, q[2] + 0.1 * bump] + }; + let (e, nu) = (100.0, 0.3); + let neo = fieldcraft_fem::constitutive::MaterialModel::from_model(&iso(e, nu), Some(&Hyperelastic::NeoHookean {})).unwrap(); + let tau = neo.response(&f, &[], &mut []).unwrap().stress; + let j = fieldcraft_fem::dense::det3(&f); + let hybrid = Formulation { hybrid: true, ..Formulation::default() }; + let mut cases: Vec<(ElementKind, Formulation, &str)> = variants().into_iter().filter(|v| v.1 != incompatible()).collect(); + cases.extend([ + (ElementKind::Hex8, hybrid, "hex8 hybrid"), + (ElementKind::Hex20, hybrid, "hex20 hybrid"), + (ElementKind::Tet10, hybrid, "tet10 hybrid"), + ]); + for (kind, formulation, name) in cases { + let mut g = grid([2, 2, 2], kind, distort); + g.define(iso(e, nu), None, formulation, nonlinear()); + neo_hookean(&mut g, e, nu); + let boundary = |q: [f64; 3]| q.iter().any(|c| near(*c, 0.0) || near(*c, 1.0)); + prescribe(&mut g, boundary, affine); + let r = solve(&g); + for (i, label) in r.nodes.iter().enumerate() { + let want = affine(g.model.node(*label).unwrap().xyz); + assert!((0..3).all(|k| (r.displacement[i][k] - want[k]).abs() < 1e-9), "{name} node {label}: {:?} vs {want:?}", r.displacement[i]); + } + for el in &r.elements { + for p in &el.points { + assert!((0..6).all(|k| (p.stress[k] - tau[k] / j).abs() < 1e-8 * e), "{name}: {:?} vs {:?}", p.stress, tau.map(|t| t / j)); + } + } + } +} + +/// A block turned rigidly by 90° about an inclined axis, its boundary prescribed: no strain, no +/// stress, no reaction, whatever the path (objectivity of the formulation). +#[test] +fn rigid_rotation_strains_and_stresses_nothing() { + let axis = [1.0 / 3f64.sqrt(); 3]; + let k = fieldcraft_fem::rotation::skew(&axis); + let k2 = fieldcraft_fem::rotation::mul(&k, &k); + let rot: [[f64; 3]; 3] = std::array::from_fn(|i| std::array::from_fn(|j| if i == j { 1.0 } else { 0.0 } + k[i][j] + k2[i][j])); + for (kind, formulation, name) in + [(ElementKind::Hex20, full(), "hex20"), (ElementKind::Tet10, full(), "tet10"), (ElementKind::Hex8, full(), "hex8")] + { + let mut g = grid([2, 2, 2], kind, |x| x); + g.define(iso(1000.0, 0.3), None, formulation, nonlinear()); + neo_hookean(&mut g, 1000.0, 0.3); + let boundary = |q: [f64; 3]| q.iter().any(|c| near(*c, 0.0) || near(*c, 1.0)); + prescribe(&mut g, boundary, |x| std::array::from_fn(|i| (0..3).map(|j| rot[i][j] * x[j]).sum::() - x[i])); + let r = solve(&g); + let worst = r.elements.iter().flat_map(|e| &e.points).flat_map(|p| p.stress.iter().chain(&p.strain)).fold(0.0f64, |m, v| m.max(v.abs())); + assert!(worst < 1e-10 * 1000.0, "{name}: largest stress or strain {worst:e}"); + assert!(r.reaction.iter().flatten().all(|v| v.abs() < 1e-9), "{name}"); + } +} + +/// Plane strain and axisymmetric 2D elements reproduce the closed forms: a plane strain block +/// stretched in x (σ_yy = 0, λ_z = 1) and a solid cylinder stretched along its axis (the 3D +/// uniaxial case). +#[test] +fn plane_strain_and_axisymmetric_stretch() { + let (e, nu, lam) = (10.0, 0.3, 1.6); + let (mu, kappa) = (e / (2.0 * (1.0 + nu)), e / (3.0 * (1.0 - 2.0 * nu))); + // Plane strain: σ_yy(l) = μ J^(−5/3) (l² − I₁/3) + κ (J − 1) = 0 with J = λ l, I₁ = λ² + l² + 1. + let syy = |l: f64| { + let (j, i1) = (lam * l, lam * lam + l * l + 1.0); + mu * j.powf(-5.0 / 3.0) * (l * l - i1 / 3.0) + kappa * (j - 1.0) + }; + let (mut lo, mut hi) = (0.1, 1.5); + for _ in 0..200 { + let mid = 0.5 * (lo + hi); + if syy(mid) > 0.0 { hi = mid } else { lo = mid } + } + let l_plane = 0.5 * (lo + hi); + let (l_axi, sigma_axi) = neo_hookean_uniaxial(mu, kappa, lam); + for (section, kind) in [ + (SectionKind::PlaneStrain { thickness: 1.0 }, ElementKind::Quad8), + (SectionKind::PlaneStrain { thickness: 1.0 }, ElementKind::Tri3), + (SectionKind::Axisymmetric {}, ElementKind::Quad8), + (SectionKind::Axisymmetric {}, ElementKind::Quad4), + ] { + let mut m = Model::default(); + let n = 3usize; + let label = |i: usize, j: usize| (1 + i + (2 * n + 1) * j) as u32; + for j in 0..=2 * n { + for i in 0..=2 * n { + m.set_node(label(i, j), [i as f64 / (2 * n) as f64, j as f64 / (2 * n) as f64, 0.0]).unwrap(); + } + } + let mut next = 1; + for cj in 0..n { + for ci in 0..n { + let (i, j) = (2 * ci, 2 * cj); + let c = [label(i, j), label(i + 2, j), label(i + 2, j + 2), label(i, j + 2)]; + match kind { + ElementKind::Quad8 => { + let mid = [label(i + 1, j), label(i + 2, j + 1), label(i + 1, j + 2), label(i, j + 1)]; + m.set_element(next, kind, c.into_iter().chain(mid).collect()).unwrap(); + next += 1; + } + ElementKind::Quad4 => { + m.set_element(next, kind, c.to_vec()).unwrap(); + next += 1; + } + _ => { + m.set_element(next, kind, vec![c[0], c[1], c[2]]).unwrap(); + m.set_element(next + 1, kind, vec![c[0], c[2], c[3]]).unwrap(); + next += 2; + } + } + } + } + let all: Vec = m.elements().keys().copied().collect(); + m.add_set("all", Set { elements: all, ..Set::default() }).unwrap(); + let mat = Material { elastic: Some(iso(e, nu)), hyperelastic: Some(Hyperelastic::NeoHookean {}), ..Material::default() }; + m.add_material("Mat", mat).unwrap(); + m.add_section("s", Section { set: "all".into(), material: "Mat".into(), kind: section, formulation: Formulation::default() }).unwrap(); + let axisymmetric = matches!(section, SectionKind::Axisymmetric {}); + // Stretch along x (plane strain) or along the axis y (axisymmetric). + let (pulled, axis) = if axisymmetric { (1, 1) } else { (0, 0) }; + let on = |k: usize, v: f64| -> Vec { m.nodes().iter().filter(|(_, n)| near(n.xyz[k], v)).map(|(l, _)| *l).collect() }; + let (low, high, side) = (on(pulled, 0.0), on(pulled, 1.0), on(1 - pulled, 0.0)); + m.add_set("low", Set { nodes: low, ..Set::default() }).unwrap(); + m.add_set("high", Set { nodes: high, ..Set::default() }).unwrap(); + m.add_set("side", Set { nodes: side, ..Set::default() }).unwrap(); + let hold = |k: usize, v: Option| { + let mut values = [None; 6]; + values[k] = v; + values + }; + m.add_boundary_condition("low", BoundaryCondition::Displacement { target: Target::Set("low".into()), values: hold(axis, Some(0.0)) }) + .unwrap(); + m.add_boundary_condition("high", BoundaryCondition::Displacement { target: Target::Set("high".into()), values: hold(axis, Some(lam - 1.0)) }) + .unwrap(); + m.add_boundary_condition("side", BoundaryCondition::Displacement { target: Target::Set("side".into()), values: hold(1 - axis, Some(0.0)) }) + .unwrap(); + let mut step = Step::new("Step", nonlinear()); + step.boundary_conditions = ["low", "high", "side"].iter().map(|s| s.to_string()).collect(); + m.add_step(step, None).unwrap(); + let r = solve_nonlinear(&m, None, &Limits::default()).unwrap(); + let lateral = if axisymmetric { l_axi } else { l_plane }; + for (i, label) in r.nodes.iter().enumerate() { + let x = m.node(*label).unwrap().xyz; + let want = [(lateral - 1.0) * x[1 - pulled], (lam - 1.0) * x[pulled]]; + let (got_axial, got_lateral) = (r.displacement[i][pulled], r.displacement[i][1 - pulled]); + assert!( + (got_axial - want[1]).abs() < 1e-8 && (got_lateral - want[0]).abs() < 1e-8, + "{section:?} {kind:?} node {label}: {:?}", + r.displacement[i] + ); + } + if axisymmetric { + assert!(r.stress.iter().all(|s| (s[1] - sigma_axi).abs() < 1e-7 * sigma_axi && s[0].abs() < 1e-7 * sigma_axi), "{:?}", r.stress[0]); + } + } +} + +/// Newton converges quadratically near the solution: in the last iterations of each increment +/// the residual falls at least like the square of the previous one. +#[test] +fn newton_converges_quadratically() { + let mut g = grid([8, 1, 2], ElementKind::Hex20, |q| [10.0 * q[0], q[1], q[2]]); + let mut step = NonlinearStatic { initial_increment: 0.5, line_search: false, ..NonlinearStatic::default() }; + step.tolerances.force = 1e-12; + step.tolerances.displacement = 1e-12; + step.tolerances.energy = 1e-20; + g.define(iso(1000.0, 0.0), None, full(), Procedure::NonlinearStatic(step)); + g.fixed("root", |q| near(q[0], 0.0)); + g.face_set("tip", |q| near(q[0], 1.0)); + // P L² / (E I) = 1: a tip rotation of about 0.46 rad. + let p = 1000.0 * (2.0f64.powi(3) / 12.0) / 100.0; + g.load("push", Load::Force { target: Target::Set("tip".into()), vector: [0.0, 0.0, -p], distribution: Distribution::Total }); + let r = solve(&g); + let h = r.nonlinear.unwrap(); + // The order of convergence ln(r₂/r₁) / ln(r₁/r₀) over successive residuals in the + // asymptotic range (below 10⁻³, above round-off): 2 for Newton with the consistent tangent. + let mut best = 0.0f64; + for inc in h.increments.iter().filter(|i| i.converged) { + let res: Vec = inc.iterations.iter().map(|i| i.force).collect(); + let orders: Vec = res.windows(3).filter(|w| w[0] < 1e-3 && w[2] > 1e-15).map(|w| (w[2] / w[1]).ln() / (w[1] / w[0]).ln()).collect(); + eprintln!("increment {} (λ = {}): residuals {res:?}, orders {orders:?}", inc.number, inc.load_factor); + best = orders.iter().fold(best, |m, o| m.max(*o)); + } + assert!(best > 1.8, "no quadratic convergence (best order {best})"); +} + +/// Without nonlinear geometry, and with a small load with it, the step is the linear solution. +#[test] +fn small_loads_and_small_kinematics_give_the_linear_solution() { + let build = |procedure: Procedure, scale: f64| { + let mut g = grid([6, 1, 2], ElementKind::Hex20, |q| [6.0 * q[0], q[1], q[2]]); + g.define(iso(1000.0, 0.3), None, full(), procedure); + g.fixed("root", |q| near(q[0], 0.0)); + g.face_set("tip", |q| near(q[0], 1.0)); + g.load("push", Load::Force { target: Target::Set("tip".into()), vector: [0.0, 0.0, -scale], distribution: Distribution::Total }); + g.face_set("top", |q| near(q[2], 1.0)); + g.load("press", Load::Pressure { set: "top".into(), magnitude: 0.1 * scale }); + g + }; + let linear = solve_static(&build(Procedure::Static {}, 1e-3).model, None, &Limits::default()).unwrap(); + let small = NonlinearStatic { nonlinear_geometry: false, ..NonlinearStatic::default() }; + let r = solve(&build(Procedure::NonlinearStatic(small), 1e-3)); + let tip = linear.displacement.iter().fold(0.0f64, |m, u| m.max(u[2].abs())); + for (a, b) in r.displacement.iter().zip(&linear.displacement) { + assert!((0..3).all(|k| (a[k] - b[k]).abs() < 1e-9 * tip), "{a:?} vs {b:?}"); + } + let r = solve(&build(nonlinear(), 1e-3)); + for (a, b) in r.displacement.iter().zip(&linear.displacement) { + assert!((0..3).all(|k| (a[k] - b[k]).abs() < 1e-3 * tip), "geometric effects of order u/L: {a:?} vs {b:?}"); + } +} + +/// Hostile and impossible input: the step's controls are refused when defined; an element +/// that inverts makes the step cut back and stop at the minimum increment with the element +/// named; the caps stop a step; unsupported elements and loads are refused by name. +#[test] +fn nonlinear_steps_fail_cleanly() { + let base = || { + let mut g = grid([2, 1, 1], ElementKind::Hex8, |q| q); + let step = NonlinearStatic { min_increment: 1e-3, ..NonlinearStatic::default() }; + g.define(iso(100.0, 0.3), None, full(), Procedure::NonlinearStatic(step)); + neo_hookean(&mut g, 100.0, 0.3); + g.fixed("root", |q| near(q[0], 0.0)); + g + }; + // Squash the bar through itself: past −100 % every element inverts. + let mut g = base(); + g.node_set("x1", |q| near(q[0], 1.0)); + g.bc("crush", BoundaryCondition::Displacement { target: Target::Set("x1".into()), values: [Some(-1.5), Some(0.0), Some(0.0), None, None, None] }); + let err = solve_nonlinear(&g.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(err.contains("stopped at load factor") && err.contains("element"), "{err}"); + // A tight cap on the Newton iterations. + let err = solve_nonlinear(&g.model, None, &Limits { max_newton_iterations: 3, ..Limits::default() }).unwrap_err().to_string(); + assert!(err.contains("Newton iterations"), "{err}"); + // The increment cap. + let mut capped = base(); + capped.node_set("x1", |q| near(q[0], 1.0)); + capped.bc("pull", BoundaryCondition::Displacement { target: Target::Set("x1".into()), values: [Some(0.5), None, None, None, None, None] }); + let mut step = capped.model.step("Step").unwrap().clone(); + step.set_procedure(Procedure::NonlinearStatic(NonlinearStatic { + initial_increment: 0.01, + max_increment: 0.01, + min_increment: 0.01, + max_increments: 3, + ..NonlinearStatic::default() + })); + capped.model.replace_step(step).unwrap(); + let err = solve_nonlinear(&capped.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(err.contains("increments reached the step's cap"), "{err}"); + // Hostile controls never reach the solver. + let d = NonlinearStatic::default(); + for bad in [ + NonlinearStatic { initial_increment: f64::NAN, ..d.clone() }, + NonlinearStatic { initial_increment: 0.0, ..d.clone() }, + NonlinearStatic { max_iterations: u32::MAX, ..d.clone() }, + NonlinearStatic { max_increments: 0, ..d.clone() }, + NonlinearStatic { tolerances: Tolerances { force: -1.0, ..Tolerances::default() }, ..d.clone() }, + ] { + let mut step = capped.model.step("Step").unwrap().clone(); + step.set_procedure(Procedure::NonlinearStatic(bad.clone())); + assert!(capped.model.replace_step(step).is_err(), "{bad:?}"); + } + // Refused by name: incompatible modes, a hybrid St Venant–Kirchhoff element, temperature. + let mut im = grid([2, 1, 1], ElementKind::Hex8, |q| q); + im.define(iso(100.0, 0.3), None, incompatible(), nonlinear()); + im.fixed("root", |q| near(q[0], 0.0)); + im.face_set("end", |q| near(q[0], 1.0)); + im.load("p", Load::Pressure { set: "end".into(), magnitude: 1.0 }); + let err = solve_nonlinear(&im.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(err.contains("incompatible modes in nonlinear steps"), "{err}"); + let mut hy = grid([2, 1, 1], ElementKind::Hex8, |q| q); + hy.define(iso(100.0, 0.3), None, Formulation { hybrid: true, ..Formulation::default() }, nonlinear()); + hy.fixed("root", |q| near(q[0], 0.0)); + hy.face_set("end", |q| near(q[0], 1.0)); + hy.load("p", Load::Pressure { set: "end".into(), magnitude: 1.0 }); + let err = solve_nonlinear(&hy.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(err.contains("hybrid elements of a St Venant–Kirchhoff material"), "{err}"); + hy.load("hot", Load::Temperature { target: Target::Set("root".into()), value: 50.0 }); + neo_hookean(&mut hy, 100.0, 0.3); + let err = solve_nonlinear(&hy.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(err.contains("temperature loads"), "{err}"); + // A model that can move freely: the first tangent is singular, named like the linear solver. + let mut free = grid([2, 1, 1], ElementKind::Hex8, |q| q); + free.define(iso(100.0, 0.3), None, full(), nonlinear()); + free.face_set("end", |q| near(q[0], 1.0)); + free.load("p", Load::Pressure { set: "end".into(), magnitude: 1.0 }); + let err = solve_nonlinear(&free.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(err.contains("singular at node"), "{err}"); + // A linear static step is not a nonlinear one. + let mut lin = grid([1, 1, 1], ElementKind::Hex8, |q| q); + lin.define(iso(100.0, 0.3), None, full(), Procedure::Static {}); + assert!(solve_nonlinear(&lin.model, None, &Limits::default()).is_err()); +} + +/// The same model gives the same numbers bit for bit, whatever the thread count. +#[test] +fn nonlinear_results_are_deterministic() { + let mut g = grid([3, 1, 1], ElementKind::Tet10, |q| [4.0 * q[0], q[1], q[2]]); + g.define(iso(100.0, 0.3), None, full(), nonlinear()); + neo_hookean(&mut g, 100.0, 0.3); + g.fixed("root", |q| near(q[0], 0.0)); + g.face_set("top", |q| near(q[2], 1.0)); + g.load("press", Load::Pressure { set: "top".into(), magnitude: 0.2 }); + let a = fieldcraft_parallel::with_threads(1, || solve(&g)); + for threads in [2, 7] { + assert_eq!(fieldcraft_parallel::with_threads(threads, || solve(&g)), a, "{threads} threads"); + } +} diff --git a/crates/ui-egui/src/results.rs b/crates/ui-egui/src/results.rs index cb3e8e2..b47ed96 100644 --- a/crates/ui-egui/src/results.rs +++ b/crates/ui-egui/src/results.rs @@ -73,10 +73,15 @@ impl Plot { } } - /// An explicit or transient heat step's frame shown, its time and the number of frames. + /// An explicit or transient heat step's frame shown, its time (a nonlinear static step's: + /// its load factor) and the number of frames. pub fn frame(&self) -> Option<(usize, f64, usize)> { let (n, time): (usize, &dyn Fn(usize) -> Option) = match self.result.kind() { Kind::Explicit(e) => (e.frames.len(), &|k| e.frames.get(k).map(|f| f.time)), + Kind::Static(r) if r.nonlinear.is_some() => { + let frames = r.nonlinear.as_ref().map(|h| h.frames.as_slice()).unwrap_or_default(); + (frames.len(), &|k| frames.get(k).map(|f| f.load_factor)) + } Kind::Heat(h) if h.times.len() > 1 => (h.times.len(), &|k| h.times.get(k).copied()), _ => return None, }; @@ -109,6 +114,8 @@ pub fn swing(result: &StepResult, frame: Option) -> f64 { let a = std::f64::consts::TAU * k as f64 / FRAMES as f64; match result.kind() { Kind::Frequency(_) => a.sin(), + // A nonlinear step plays its increments' frames. + Kind::Static(r) if r.nonlinear.is_some() => 1.0, Kind::Static(_) => 0.5 * (1.0 - a.cos()), Kind::Explicit(_) | Kind::Fluid(_) | Kind::Heat(_) => 1.0, } @@ -118,10 +125,13 @@ pub fn swing(result: &StepResult, frame: Option) -> f64 { /// step running through its field frames once per cycle. pub fn pick(view: &ResultsView, result: &StepResult, frame: Option) -> Pick { let mut p = view.pick(); - if let (Some(k), Kind::Explicit(e)) = (frame, result.kind()) - && let Some(last) = e.frames.len().checked_sub(1) - { - p.frame = Some((k as usize * e.frames.len() / FRAMES as usize).min(last)); + let frames = match result.kind() { + Kind::Explicit(e) => e.frames.len(), + Kind::Static(r) => r.nonlinear.as_ref().map_or(0, |h| h.frames.len()), + _ => 0, + }; + if let (Some(k), Some(last)) = (frame, frames.checked_sub(1)) { + p.frame = Some((k as usize * frames / FRAMES as usize).min(last)); } p } diff --git a/crates/ui-egui/src/results_panel.rs b/crates/ui-egui/src/results_panel.rs index 6ab1c86..9ee918a 100644 --- a/crates/ui-egui/src/results_panel.rs +++ b/crates/ui-egui/src/results_panel.rs @@ -162,6 +162,17 @@ pub fn controls(app: &mut ModelApp, plot: &Plot, ui: &mut egui::Ui, t: &Tokens) if let Kind::Explicit(e) = plot.result.kind() { ui.separator(); explicit(ui, e, plot, t, &mut run); + } else if let Some((k, factor, n)) = plot.frame().filter(|_| matches!(plot.result.kind(), Kind::Static(_))) { + // A nonlinear static step: its increments, the deformed shape following the slider. + ui.separator(); + ui.horizontal(|ui| { + ui.label("Increment"); + let mut frame = k; + if ui.add(egui::Slider::new(&mut frame, 0..=n.saturating_sub(1))).changed() { + run.push(("results.frame", json!({"frame": frame}))); + } + }); + ui.label(RichText::new(format!("load factor {}", number(factor))).color(t.text_dim)); } else if let Kind::Frequency(f) = plot.result.kind() { ui.separator(); ui.label(RichText::new("Modes").color(t.text_dim)); diff --git a/docs/features.md b/docs/features.md index 8c32ea8..0cab20d 100644 --- a/docs/features.md +++ b/docs/features.md @@ -56,7 +56,7 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands. | Panel | Capability | Status | Command | |---|---|---|---| -| Materials | Material — Elastic (isotropic, orthotropic) material with density | ready | `model.material.create`, `model.material.edit`, `model.material.delete` | +| Materials | Material — Elastic (isotropic, orthotropic) material with density; neo-Hookean hyperelasticity for nonlinear steps | ready | `model.material.create`, `model.material.edit`, `model.material.delete` | | Materials | Plasticity — Von Mises plasticity: bilinear, hardening table or Johnson–Cook | ready | `model.material.plasticity` | | Materials | Hyperelastic — Hyperelastic material models | coming soon (M8) | | | Materials | Damage — Johnson–Cook damage and element erosion | ready | `model.material.damage` | @@ -78,7 +78,7 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands. | Panel | Capability | Status | Command | |---|---|---|---| -| Steps | Step — A linear static or frequency step with its output | ready | `model.step.create`, `model.step.edit`, `model.step.delete` | +| Steps | Step — A linear or nonlinear static, frequency, heat or explicit step with its output | ready | `model.step.create`, `model.step.edit`, `model.step.delete` | | Steps | Activate — Choose which loads and boundary conditions act in a step | ready | `model.step.activate` | | Steps | Step Table — The step by load and boundary condition table | ready | `model.activity` | | Procedures | Amplitude — Time and frequency amplitudes | coming soon (M7) | | @@ -122,7 +122,7 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands. | Decks | Nastran Deck — Write Nastran bulk data (SOL 101, 103) | ready | `deck.write {"format":"nastran"}` | | Decks | LS-DYNA Deck — Write an LS-DYNA keyword file | ready | `deck.write {"format":"lsdyna"}` | | Decks | Read Deck — Read a deck, keeping keywords FieldCraft does not know | coming soon (M9) | | -| Implicit | Static — Linear static (or coupled temperature–displacement) solve of a step with FieldCraft's own solver; heat transfer in the drop-down | ready | `solve.static`, `solve.heat` | +| Implicit | Static — Linear static (or coupled temperature–displacement) solve of a step with FieldCraft's own solver; nonlinear static and heat transfer in the drop-down | ready | `solve.static`, `solve.nonlinear`, `solve.heat` | | Implicit | Frequency — Natural frequencies and mode shapes of a frequency step | ready | `solve.frequency` | | Explicit | Explicit — Transient dynamics of an explicit step by central differences (contact and erosion coming in M10) | ready | `solve.explicit` | | Explicit | Blast & Impact — EOS-driven fluids with fluid–structure contact (ALE or SPH) | coming soon (M14) | |