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