From 1b57e2e672212fb63d4e8dcccbd1e1d1016bceb6 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 13:09:39 -0400 Subject: [PATCH] solver, model, engine: linear buckling steps MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Procedure::Buckle {modes, shift}: the step's loads as a linear preload, K_σ from its stresses, and (K + λ K_σ) φ = 0 by block Lanczos in the K inner product: K⁻¹(−K_σ) on Cholesky without a shift, the buckling transformation (K + σ K_σ)⁻¹ K on LDLᵀ with one, whose inertia counts the load factors below σ (Grimes–Lewis–Simon 1994). Lanczos gains an operand and an eigenvalue map (defaults keep frequency steps bit for bit) and reports converged negative eigenvalues. Results are frequency results with a `buckling` part; modes display with their load factor. solve.buckle (in the Static drop-down); Abaqus *BUCKLE, Nastran SOL 105 and LS-DYNA refuse it by name. Euler columns of hex20 within 0.18 % of Engesser's load. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/decks/src/abaqus.rs | 28 ++++ crates/decks/src/lsdyna/mod.rs | 2 +- crates/decks/src/nastran/mod.rs | 3 + crates/decks/tests/abaqus.rs | 17 +++ crates/engine/src/catalogue.rs | 4 +- crates/engine/src/cmd/mod.rs | 1 + crates/engine/src/cmd/solve.rs | 31 +++- crates/model/src/check.rs | 3 +- crates/model/src/steps.rs | 19 ++- crates/solver/src/analysis.rs | 26 ++-- crates/solver/src/buckle.rs | 213 ++++++++++++++++++++++++++++ crates/solver/src/eigen.rs | 36 ++++- crates/solver/src/lib.rs | 2 + crates/solver/src/system.rs | 8 ++ crates/solver/tests/buckle.rs | 93 ++++++++++++ crates/ui-egui/src/results.rs | 13 +- crates/ui-egui/src/results_panel.rs | 12 +- docs/features.md | 2 +- docs/fem.md | 20 +++ docs/solvers/literature.md | 5 + 20 files changed, 511 insertions(+), 27 deletions(-) create mode 100644 crates/solver/src/buckle.rs create mode 100644 crates/solver/tests/buckle.rs diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index e049cee..79a28e3 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -1216,6 +1216,34 @@ impl Deck { base_bcs.clone_from(&step.boundary_conditions); base_loads.clone_from(&step.loads); } + Procedure::Buckle { modes, shift } => { + // Abaqus's buckling base state is the end of the previous general step; + // FieldCraft's steps all start from the undeformed, unloaded model. + if !base_loads.is_empty() + || steps + .iter() + .take_while(|s| !std::ptr::eq(*s, step)) + .any(|s| !matches!(s.procedure, Procedure::Frequency { .. } | Procedure::Buckle { .. })) + { + return Err(Error::Unsupported { + format: FORMAT, + what: format!( + "step `{}`: a buckling step after a general step (Abaqus would buckle from that step's loaded state)", + step.name + ), + }); + } + if shift.is_some() { + return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: a buckling shift", step.name) }); + } + self.line(&format!("*STEP, NAME={}", step.name)); + self.line("*BUCKLE, EIGENSOLVER=LANCZOS"); + self.line(&format!("{modes}")); + self.step_boundaries(model, r, &base_bcs, &step.boundary_conditions)?; + let loads = StepLoads::new(model, r, step, &everything); + self.step_loads(&loads, written)?; + written = loads.written(); + } Procedure::Frequency { modes, min_hz, max_hz } => { self.line(&format!("*STEP, NAME={}", step.name)); self.line("*FREQUENCY, EIGENSOLVER=LANCZOS"); diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 6da44cf..43bb79a 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -899,7 +899,7 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste d.card([Field::Blank, Field::Blank, Field::Blank, Field::Blank, Field::from(1)]); } } - Procedure::ExplicitDynamic(_) | Procedure::HeatTransfer { .. } | Procedure::CoupledTemperatureDisplacement {} => { + Procedure::ExplicitDynamic(_) | Procedure::HeatTransfer { .. } | Procedure::CoupledTemperatureDisplacement {} | Procedure::Buckle { .. } => { return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: {} steps", step.name, step.procedure.name()) }); } } diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index 3e53e12..7669618 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -280,6 +280,9 @@ fn solution(steps: &[&Step]) -> std::result::Result { if let Some(s) = steps.iter().find(|s| matches!(s.procedure, Procedure::ExplicitDynamic(_))) { return Err(format!("step `{}`: explicit dynamic steps (Nastran's linear solutions have none)", s.name)); } + if let Some(s) = steps.iter().find(|s| matches!(s.procedure, Procedure::Buckle { .. })) { + return Err(format!("step `{}`: buckling steps (Nastran's SOL 105; this writer writes SOL 101 and 103 so far)", s.name)); + } if let Some(s) = steps.iter().find(|s| matches!(s.procedure, Procedure::NonlinearStatic(_))) { return Err(format!("step `{}`: nonlinear static steps (Nastran's SOL 400; this writer writes SOL 101 and 103 so far)", s.name)); } diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs index 54c694f..6da5b7b 100644 --- a/crates/decks/tests/abaqus.rs +++ b/crates/decks/tests/abaqus.rs @@ -828,3 +828,20 @@ fn nonlinear_plastic_steps_and_what_is_refused() { m.add_step(step("Again", Procedure::NonlinearStatic(NonlinearStatic::default()), &["clamp"], &["end"]), None).unwrap(); assert!(unsupported(&m).contains("plastic material `Steel` in a deck of 2 general steps")); } + +/// A buckling step is `*BUCKLE` with its loads as the perturbation loads; after a general step +/// (Abaqus's base state would be that step's) and with a shift it is refused by name. +#[test] +fn buckling_steps() { + let mut m = cantilever_mesh(); + m.add_load("end", Load::Force { target: Target::Set("tip".into()), vector: [-100.0, 0.0, 0.0], distribution: Distribution::Total }).unwrap(); + m.add_step(step("Buckle", Procedure::Buckle { modes: 4, shift: None }, &["clamp"], &["end"]), None).unwrap(); + let deck = abaqus::write(&m).unwrap(); + assert!(deck.contains("*BUCKLE, EIGENSOLVER=LANCZOS\n4\n") && deck.contains("*CLOAD"), "{deck}"); + let mut shifted = m.clone(); + shifted.replace_step(step("Buckle", Procedure::Buckle { modes: 4, shift: Some(2.0) }, &["clamp"], &["end"])).unwrap(); + assert!(abaqus::write(&shifted).unwrap_err().to_string().contains("shift")); + m.add_step(step("Bend", Procedure::Static {}, &["clamp"], &["end"]), Some("Buckle")).unwrap(); + assert!(abaqus::write(&m).unwrap_err().to_string().contains("after a general step")); + assert!(fieldcraft_decks::lsdyna::write(&m).is_err() && fieldcraft_decks::nastran::write(&m).is_err()); +} diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index 3515eab..e37cc8a 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -470,8 +470,8 @@ pub static CATALOGUE: &[Capability] = &[ "Static", Solve, "Implicit", - "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"], + "Linear static (or coupled temperature–displacement) solve of a step with FieldCraft's own solver; nonlinear static, buckling and heat transfer in the drop-down", + &["solve.nonlinear", "solve.buckle", "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 736d731..2d104e9 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -109,6 +109,7 @@ pub static ALL: &[CommandSpec] = &[ solve::STATIC, solve::NONLINEAR, solve::FREQUENCY, + solve::BUCKLE, solve::HEAT, solve::EXPLICIT, solve::HISTORY, diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 8903881..3ada6dc 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -7,8 +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_nonlinear, - solve_static, + FrequencyOptions, Limits, Location, MassKind, StaticResult, StepResult as ImplicitResult, solve_buckle, solve_frequency, solve_heat, + solve_nonlinear, solve_static, }; use serde_json::{Value, json}; @@ -60,6 +60,14 @@ pub const FREQUENCY: CommandSpec = CommandSpec { run: run_frequency, }; +pub const BUCKLE: CommandSpec = CommandSpec { + id: "solve.buckle", + label: "Solve Buckling Step", + params: r#"{"step": "Buckle", "max_dofs": 6000000, "max_factor_entries": 400000000, "max_iterations": 3000, "threads": 8} — linear buckling of a buckling step (the only one when `step` is left out): the step's loads as a linear preload, then the lowest load factors λ with (K + λ K_σ) φ = 0 and their mode shapes (largest component 1), or those above the step's shift (whose LDLᵀ factor counts the load factors below it); negative load factors found are reported; the result stays in the session and its modes display like frequency modes; `threads` as for solve.static"#, + kind: CommandKind::Model, + run: run_buckle, +}; + pub const EXPLICIT: CommandSpec = CommandSpec { id: "solve.explicit", label: "Solve Explicit Dynamic Step", @@ -159,6 +167,7 @@ fn run_all(s: &mut Session, _: &Value) -> Result { Procedure::Static {} | Procedure::CoupledTemperatureDisplacement {} => run_static, Procedure::HeatTransfer { .. } => run_heat, Procedure::Frequency { .. } => run_frequency, + Procedure::Buckle { .. } => run_buckle, Procedure::ExplicitDynamic(_) => run_explicit, Procedure::Fluid { .. } => run_fluid, Procedure::NonlinearStatic(_) => run_nonlinear, @@ -292,6 +301,24 @@ fn run_frequency(s: &mut Session, p: &Value) -> Result { Ok(summary) } +fn run_buckle(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_buckle(model, step.as_deref(), &limits))?; + let summary = json!({ + "step": r.step, + "dofs": r.stats.dofs, + "equations": r.stats.equations, + "lanczos_vectors": r.stats.lanczos_vectors, + "modes": r.modes.iter().enumerate().map(|(i, m)| json!({"mode": i + 1, "load_factor": m.eigenvalue})).collect::>(), + "buckling": r.buckling, + "warnings": r.warnings, + }); + s.store_result(StepResult::Implicit(ImplicitResult::Frequency(r))); + Ok(summary) +} + fn run_explicit(s: &mut Session, p: &Value) -> Result { use fieldcraft_explicit::Limits as ExplicitLimits; let limits = ExplicitLimits { diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs index 15ad54b..effdaff 100644 --- a/crates/model/src/check.rs +++ b/crates/model/src/check.rs @@ -322,6 +322,7 @@ impl Model { Procedure::ExplicitDynamic(_) => "is dynamic", Procedure::HeatTransfer { transient: Some(_) } => "is transient heat transfer", Procedure::Static {} + | Procedure::Buckle { .. } | Procedure::NonlinearStatic(_) | Procedure::Fluid { .. } | Procedure::HeatTransfer { transient: None } @@ -486,7 +487,7 @@ impl Model { } } match &s.procedure { - Procedure::Static {} | Procedure::NonlinearStatic(_) => { + Procedure::Static {} | Procedure::NonlinearStatic(_) | Procedure::Buckle { .. } => { if s.loads.is_empty() && !prescribes_motion { d.warning(&subject, "nothing loads the model in this step"); } diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index 67e808e..73ed648 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -55,6 +55,14 @@ pub enum Procedure { #[serde(default, skip_serializing_if = "Option::is_none")] max_hz: Option, }, + /// Linear buckling: the load factors λ at which the step's loads (and prescribed motions), + /// applied as a linear preload, make K + λ K_σ singular, and their mode shapes; the lowest + /// `modes` positive ones, or those nearest above `shift` when one is given. + Buckle { + modes: u32, + #[serde(default, skip_serializing_if = "Option::is_none")] + shift: Option, + }, /// Transient dynamics by explicit time integration up to `end_time`. ExplicitDynamic(ExplicitDynamic), /// Time-dependent incompressible flow in a fluid domain, from rest (or the domain's initial @@ -546,6 +554,7 @@ impl Procedure { match self { Procedure::Static {} => "static", Procedure::Frequency { .. } => "frequency", + Procedure::Buckle { .. } => "buckle", Procedure::ExplicitDynamic(_) => "explicit_dynamic", Procedure::HeatTransfer { .. } => "heat_transfer", Procedure::CoupledTemperatureDisplacement {} => "coupled_temperature_displacement", @@ -618,6 +627,14 @@ impl Procedure { if let Procedure::HeatTransfer { transient: Some(t) } = self { return t.validate(); } + if let Procedure::Buckle { modes, shift } = *self { + if !(1..=MAX_MODES).contains(&modes) { + return Err(Error::Invalid(format!("a buckling step extracts 1 to {MAX_MODES} modes, got {modes}"))); + } + if let Some(s) = shift { + finite("buckling shift", s)?; + } + } if let Procedure::Frequency { modes, min_hz, max_hz } = *self { if !(1..=MAX_MODES).contains(&modes) { return Err(Error::Invalid(format!("a frequency step extracts 1 to {MAX_MODES} modes, got {modes}"))); @@ -640,7 +657,7 @@ impl Procedure { pub fn default_output(&self) -> BTreeSet { match self { Procedure::Static {} | Procedure::NonlinearStatic(_) => OutputVariable::ALL.into_iter().collect(), - Procedure::Frequency { .. } => [OutputVariable::Displacement].into_iter().collect(), + Procedure::Frequency { .. } | Procedure::Buckle { .. } => [OutputVariable::Displacement].into_iter().collect(), Procedure::ExplicitDynamic(_) => [OutputVariable::Displacement, OutputVariable::Stress].into_iter().collect(), // Velocity, pressure and wall forces, always. Procedure::Fluid { .. } => BTreeSet::new(), diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 78cbfa8..ec57169 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -147,6 +147,9 @@ pub struct FrequencyResult { pub total_mass: f64, pub stats: Stats, pub warnings: Vec, + /// For a buckling step: its modes' eigenvalues are load factors (`frequency`, participation + /// and total mass are zero) and this says what else the eigenproblem found. + pub buckling: Option, } /// Mass matrix choice for a frequency analysis. @@ -168,8 +171,8 @@ pub struct FrequencyOptions { } /// The assembled matrices and the model's mass. -struct Assembled { - k: Csr, +pub(crate) struct Assembled { + pub(crate) k: Csr, m: Option, } @@ -180,7 +183,7 @@ pub(crate) const ROUND: usize = 4096; 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 { +pub(crate) fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option) -> Result { let layout = &mesh.layout; let mut k = graph.zero_matrix(layout); let mut m = mass.map(|_| graph.zero_matrix(layout)); @@ -540,10 +543,17 @@ impl Operator for ShiftInvert<'_> { } /// Block size of the Lanczos iteration: enough for the six rigid-body modes of a free body. -const BLOCK: usize = 6; +pub(crate) const BLOCK: usize = 6; + +/// Most Lanczos vectors for `neq` equations: the basis (vectors and their mass products) costs +/// 2 n × vectors, capped with the factor's entry limit. +pub(crate) fn lanczos_basis(limits: &Limits, neq: usize) -> usize { + let by_memory = limits.max_factor_entries / (2 * neq).max(1); + limits.max_lanczos_vectors.min(by_memory).max(BLOCK.min(neq)) +} /// Relative residual of an accepted Ritz pair. -const TOLERANCE: f64 = 1e-9; +pub(crate) const TOLERANCE: f64 = 1e-9; /// Solve a frequency step (the step named `step`, or the model's only frequency step): the /// lowest natural modes within the step's frequency range. @@ -613,10 +623,7 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt min: min_hz.map_or(f64::NEG_INFINITY, |f| (two_pi * f).powi(2)), max: max_hz.map_or(f64::INFINITY, |f| (two_pi * f).powi(2)), }; - // The basis (vectors and their mass products) costs 2 n × vectors; cap it with the factor's - // entry limit. - let by_memory = limits.max_factor_entries / (2 * neq).max(1); - let max_basis = limits.max_lanczos_vectors.min(by_memory).max(BLOCK.min(neq)); + let max_basis = lanczos_basis(limits, neq); let op = ShiftInvert { m: &mff, factor: &chol }; let pairs = lanczos(&op, shift, wanted, BLOCK, max_basis, TOLERANCE)?; // Participation: Γ_d = φᵀ M r_d with r_d the unit translation along d on the free equations @@ -655,6 +662,7 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt total_mass, stats: Stats { dofs: mesh.dofs(), equations: neq, factor_entries: chol.entries(), lanczos_vectors: pairs.basis }, warnings, + buckling: None, }) } diff --git a/crates/solver/src/buckle.rs b/crates/solver/src/buckle.rs new file mode 100644 index 0000000..d5d21dd --- /dev/null +++ b/crates/solver/src/buckle.rs @@ -0,0 +1,213 @@ +//! Linear buckling: the load factors λ at which a step's loads, applied as a linear preload, +//! make the structure lose stiffness, (K + λ K_σ) φ = 0, and the mode shapes φ (K.-J. Bathe, +//! *Finite Element Procedures*, 2nd ed., 2014, §6.8.2). +//! +//! **Preload.** The step is solved as a linear static step (its loads and prescribed motions); +//! the stresses at every element's stiffness points give the initial-stress stiffness +//! K_σ = ∫ (∇N_a · σ ∇N_b) I dV ([`fieldcraft_fem::finite`]). The eigenproblem holds every +//! boundary condition of the step at zero (the buckling motion is a perturbation of the +//! preloaded state). Pressures act as dead loads: the load stiffness of a follower pressure is +//! not part of K_σ. +//! +//! **Eigenproblem** by the block Lanczos of [`crate::eigen`], in the K inner product (K is +//! positive definite for a supported structure, while K_σ is indefinite), after R. G. Grimes, +//! J. G. Lewis, H. D. Simon, *A shifted block Lanczos algorithm for solving sparse symmetric +//! generalized eigenproblems*, SIAM J. Matrix Anal. Appl. 15(1):228–272, 1994, §3 (the buckling +//! transformations): +//! - **without a shift**, the operator K⁻¹(−K_σ) has eigenvalues θ = 1/λ: the lowest positive +//! load factors are its largest positive θ, well separated, and K is factored by Cholesky. +//! Converged negative θ are the load factors of the reversed load (reported, not returned as +//! modes); +//! - **with a shift σ ≠ 0**, the buckling spectral transformation (K + σ K_σ)⁻¹ K has +//! θ = λ/(λ − σ), largest for the load factors nearest above σ; K + σ K_σ is indefinite once σ +//! passes the first load factor and is factored by LDLᵀ ([`crate::ldlt`]), whose inertia counts +//! the load factors between 0 and σ (Sylvester's law), reported so none is missed silently. +//! +//! Modes are scaled to a largest component of 1 (positive). Solids and 2D elements only: +//! trusses, shells and beams have no initial-stress stiffness yet and are refused by name. + +use std::collections::BTreeMap; + +use fieldcraft_model::{Model, Procedure}; +use serde::Serialize; + +use crate::analysis::{ + BLOCK, FrequencyResult, Mode, Setup, Stats, TOLERANCE, assemble, budget, factor_error, factor_symmetric, lanczos_basis, pattern, reduction, + solve_static, split, +}; +use crate::eigen::{Operator, Wanted, lanczos}; +use crate::ldlt::{Kind, Symmetric}; +use crate::sparse::{Csr, add_block}; +use crate::system::Mesh; +use crate::{Error, Limits, Result}; + +/// What a buckling eigenproblem found besides its modes. +#[derive(Debug, Clone, Default, PartialEq, Serialize)] +pub struct Buckling { + /// Converged negative load factors (the load reversed buckles at |λ|), nearest to zero + /// first; found without a shift only. + pub negative: Vec, + /// The shift, and the number of load factors between 0 and it (the inertia of K + σ K_σ). + #[serde(skip_serializing_if = "Option::is_none")] + pub shift: Option, + #[serde(skip_serializing_if = "Option::is_none")] + pub below_shift: Option, +} + +/// K⁻¹(−K_σ) in the K inner product: θ = 1/λ. +struct Plain<'a> { + k: &'a Csr, + /// −K_σ. + g: &'a Csr, + factor: &'a Symmetric, +} + +impl Operator for Plain<'_> { + fn n(&self) -> usize { + self.k.n + } + fn mass(&self, x: &[f64], y: &mut [f64]) { + self.k.mul(x, y); + } + fn solve(&self, x: &mut [f64]) -> Result<()> { + self.factor.solve(x) + } + fn operand(&self, x: &[f64], _kx: &[f64]) -> Vec { + let mut y = vec![0.0; self.g.n]; + self.g.mul(x, &mut y); + y + } + fn eigenvalue(&self, _shift: f64, theta: f64) -> f64 { + 1.0 / theta + } +} + +/// (K + σ K_σ)⁻¹ K in the K inner product: θ = λ/(λ − σ). +struct Shifted<'a> { + k: &'a Csr, + /// K + σ K_σ and its factor, for iterative refinement. + shifted: &'a Csr, + factor: &'a Symmetric, +} + +impl Operator for Shifted<'_> { + fn n(&self) -> usize { + self.k.n + } + fn mass(&self, x: &[f64], y: &mut [f64]) { + self.k.mul(x, y); + } + fn solve(&self, x: &mut [f64]) -> Result<()> { + self.factor.solve_refined(self.shifted, x, 2) + } + fn eigenvalue(&self, shift: f64, theta: f64) -> f64 { + shift * theta / (theta - 1.0) + } +} + +/// Solve a buckling step (the step named `step`, or the model's only buckling step). +pub fn solve_buckle(model: &Model, step: Option<&str>, limits: &Limits) -> Result { + let step = model.select_step(step, |p| matches!(p, Procedure::Buckle { .. }), "buckling")?; + let Procedure::Buckle { modes, shift } = step.procedure else { + return Err(Error::Internal("buckling step".into())); + }; + step.procedure.validate()?; + let resolved = model.resolve()?; + let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; + if let Some(t) = mesh.trusses.first() { + return Err(Error::NotSupported(format!( + "trusses in buckling steps (element {}): their initial-stress stiffness is not implemented", + t.label + ))); + } + if let Some(s) = mesh.structural.first() { + return Err(Error::NotSupported(format!( + "shells and beams in buckling steps (element {}): their initial-stress stiffness is not implemented", + s.label + ))); + } + // The preload: the step solved as a linear static step. + let mut preload_model = model.clone(); + let mut preload_step = step.clone(); + preload_step.procedure = Procedure::Static {}; + preload_model.replace_step(preload_step)?; + let preload = solve_static(&preload_model, Some(&step.name), limits)?; + let mut warnings: Vec = preload.warnings.clone(); + let stresses: BTreeMap> = preload.elements.iter().map(|e| (e.label, e.points.iter().map(|p| p.stress).collect())).collect(); + + let (graph, order) = pattern(&mesh); + let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; + let red = reduction(&setup, false, &order, budget(limits, &mesh, &graph, 2), &mut Vec::new())?; + let neq = red.neq(); + if neq == 0 { + return Err(Error::Invalid("every degree of freedom is held; there is nothing to buckle".into())); + } + let k = assemble(&mesh, &graph, None)?.k; + let mut kg = graph.zero_matrix(&mesh.layout); + let matrices = fieldcraft_parallel::map_grained(mesh.elements.len(), crate::analysis::GRAIN, |i| -> Result<_> { + let e = mesh.elements.get(i).ok_or_else(|| Error::Internal("element index".into()))?; + let kernel = e.kernel()?; + let stress = stresses.get(&e.label).ok_or_else(|| Error::Internal(format!("element {} has no preload stress", e.label)))?; + Ok((kernel.dofs().node, kernel.geometric_stiffness(stress).map_err(|err| Error::element(e.label, err))?.data)) + }); + for (e, m) in mesh.elements.iter().zip(matrices) { + let (kinds, data) = m?; + add_block(&mut kg, &mesh.layout, &e.nodes, kinds, &data)?; + } + let [kff, gff] = red.reduce([&k, &kg])?; + drop((k, kg)); + if gff.val.iter().all(|v| *v == 0.0) { + return Err(Error::Invalid(format!("step `{}` puts no stress into the model, so nothing can buckle (give it loads)", step.name))); + } + let max_basis = lanczos_basis(limits, neq); + let count = usize::try_from(modes).unwrap_or(usize::MAX).min(neq); + let sigma = shift.unwrap_or(0.0); + let (pairs, entries, buckling) = if sigma == 0.0 { + let factor = factor_symmetric(&kff, limits, Kind::Definite).map_err(|e| factor_error(e, &mesh, &red))?; + let minus = Csr { val: gff.val.iter().map(|v| -v).collect(), ..gff.clone() }; + let op = Plain { k: &kff, g: &minus, factor: &factor }; + let pairs = lanczos(&op, 0.0, Wanted { count, min: 0.0, max: f64::INFINITY }, BLOCK, max_basis, TOLERANCE)?; + let buckling = Buckling { negative: pairs.negative.clone(), shift: None, below_shift: None }; + (pairs, factor.entries(), buckling) + } else { + let shifted = kff.add_scaled(sigma, &gff)?; + let factor = factor_symmetric(&shifted, limits, Kind::Indefinite).map_err(|e| match e { + crate::cholesky::FactorError::Singular { .. } => { + Error::Singular(format!("the shift {sigma} is a buckling load factor (K + σ K_σ is singular); move it a little")) + } + e => factor_error(e, &mesh, &red), + })?; + let below = factor.inertia().negative; + if below > 0 { + warnings.push(format!("{below} buckling load factor(s) lie between 0 and the shift {sigma}; the modes returned are those above it")); + } + let op = Shifted { k: &kff, shifted: &shifted, factor: &factor }; + let (min, max) = if sigma > 0.0 { (sigma, f64::INFINITY) } else { (f64::NEG_INFINITY, f64::INFINITY) }; + let pairs = lanczos(&op, sigma, Wanted { count, min, max }, BLOCK, max_basis, TOLERANCE)?; + (pairs, factor.entries(), Buckling { negative: Vec::new(), shift: Some(sigma), below_shift: Some(below) }) + }; + if let Some(first) = buckling.negative.first() { + warnings.push(format!("the reversed load buckles too: negative load factor {first:.6e} (and {} more found)", buckling.negative.len() - 1)); + } + let mut out = Vec::with_capacity(pairs.values.len()); + for (lambda, phi) in pairs.values.iter().zip(&pairs.vectors) { + // Scale to a largest component of 1 (Lanczos already made it positive). + let big = phi.iter().fold(0.0f64, |m, v| m.max(v.abs())); + let scaled: Vec = phi.iter().map(|v| if big > 0.0 { v / big } else { *v }).collect(); + let (shape, rotation, point_shape) = split(&red.expand(&scaled, false), &mesh); + out.push(Mode { eigenvalue: *lambda, frequency: 0.0, shape, rotation, participation: [0.0; 3], effective_mass: [0.0; 3], point_shape }); + } + if out.len() < count { + warnings.push(format!("found {} of the {modes} buckling modes asked for", out.len())); + } + Ok(FrequencyResult { + step: step.name.clone(), + nodes: mesh.labels.clone(), + reference_points: mesh.points.iter().map(|p| (p.name.clone(), p.label)).collect(), + modes: out, + total_mass: 0.0, + stats: Stats { dofs: mesh.dofs(), equations: neq, factor_entries: entries, lanczos_vectors: pairs.basis }, + warnings, + buckling: Some(buckling), + }) +} diff --git a/crates/solver/src/eigen.rs b/crates/solver/src/eigen.rs index b25dc93..54e8ecf 100644 --- a/crates/solver/src/eigen.rs +++ b/crates/solver/src/eigen.rs @@ -190,6 +190,17 @@ pub trait Operator: Sync { fn mass(&self, x: &[f64], y: &mut [f64]); /// x ← (K − σM)⁻¹ x. fn solve(&self, x: &mut [f64]) -> Result<()>; + /// The right-hand side the operator solves for, given a basis vector `x` and its M product + /// `mx`: M x for the spectral transformation (the default); another matrix's product for + /// other pencils (buckling, [`crate::buckle`]). + fn operand(&self, _x: &[f64], mx: &[f64]) -> Vec { + mx.to_vec() + } + /// The eigenvalue of a Ritz value θ of the operator: σ + 1/θ for the spectral + /// transformation. + fn eigenvalue(&self, shift: f64, theta: f64) -> f64 { + shift + 1.0 / theta + } } /// Converged eigenpairs: λ ascending, vectors M-orthonormal. @@ -197,6 +208,9 @@ pub trait Operator: Sync { pub struct Eigenpairs { pub values: Vec, pub vectors: Vec>, + /// Converged eigenvalues of negative Ritz values θ (outside what the transformation looks + /// for: negative buckling load factors), largest |θ| first. + pub negative: Vec, /// Lanczos vectors used. pub basis: usize, } @@ -287,7 +301,7 @@ pub fn lanczos(op: &impl Operator, shift: f64, wanted: Wanted, block: usize, max let mut w: Vec> = Vec::with_capacity(b); // The block's solves are independent: one task each. for y in fieldcraft_parallel::map(mcur.len(), |i| { - let mut y = mcur.get(i).cloned().unwrap_or_default(); + let mut y = op.operand(cur.get(i).map(Vec::as_slice).unwrap_or_default(), mcur.get(i).map(Vec::as_slice).unwrap_or_default()); op.solve(&mut y).map(|()| y) }) { w.push(y?); @@ -359,7 +373,7 @@ pub fn lanczos(op: &impl Operator, shift: f64, wanted: Wanted, block: usize, max if res2.sqrt() > tolerance * th { break; } - let lambda = shift + 1.0 / th; + let lambda = op.eigenvalue(shift, th); if lambda > wanted.max { done = true; break; @@ -384,10 +398,24 @@ pub fn lanczos(op: &impl Operator, shift: f64, wanted: Wanted, block: usize, max } } if done { + // Negative Ritz values that have converged, largest magnitude first. + let mut negative = Vec::new(); + for &i in idx.iter().rev() { + let th = theta[i]; + if !(th < 0.0) { + break; + } + let res2: f64 = + if full { 0.0 } else { (0..b).map(|r| (0..b).map(|c| bnext[r][c] * s[(m - b + c) * m + i]).sum::().powi(2)).sum() }; + if res2.sqrt() > tolerance * th.abs() { + break; + } + negative.push(op.eigenvalue(shift, th)); + } let mut values = Vec::new(); let mut vectors = Vec::new(); for &i in &accepted { - values.push(shift + 1.0 / theta[i]); + values.push(op.eigenvalue(shift, theta[i])); let mut x = vec![0.0; n]; for (k, qk) in q.iter().enumerate() { let c = s[k * m + i]; @@ -405,7 +433,7 @@ pub fn lanczos(op: &impl Operator, shift: f64, wanted: Wanted, block: usize, max } vectors.push(x); } - return Ok(Eigenpairs { values, vectors, basis: m }); + return Ok(Eigenpairs { values, vectors, negative, basis: m }); } } if full { diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index 36d35a9..3a777a3 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -28,6 +28,7 @@ #![allow(clippy::needless_range_loop)] pub mod analysis; +pub mod buckle; pub mod cholesky; pub mod constraints; pub mod discrete; @@ -43,6 +44,7 @@ pub use analysis::{ ElementStress, FrequencyOptions, FrequencyResult, Location, MassKind, Mode, PointResult, StaticResult, Stats, StepResult, solve_frequency, solve_static, }; +pub use buckle::{Buckling, solve_buckle}; pub use heat::{HeatResult, solve_heat}; pub use nonlinear::solve_nonlinear; diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 9ca9391..b7b2df7 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -155,6 +155,14 @@ impl Kernel<'_> { } } + /// The initial-stress stiffness of stresses at the stiffness points. + pub fn geometric_stiffness(&self, stress: &[[f64; 6]]) -> fieldcraft_fem::Result { + match self { + Kernel::Solid(s) => s.geometric_stiffness(stress), + Kernel::Plane(p) => p.geometric_stiffness(stress), + } + } + pub fn rule(&self) -> fieldcraft_fem::quadrature::Rule { match self { Kernel::Solid(s) => s.rule(), diff --git a/crates/solver/tests/buckle.rs b/crates/solver/tests/buckle.rs new file mode 100644 index 0000000..fe94921 --- /dev/null +++ b/crates/solver/tests/buckle.rs @@ -0,0 +1,93 @@ +//! Linear buckling against closed forms: Euler columns of hex20 solids (pinned–pinned by a +//! symmetry plane at midspan, and clamped–free), with Engesser's shear correction +//! P = P_e / (1 + P_e/(κ G A)) (S. P. Timoshenko, J. M. Gere, *Theory of Elastic Stability*, +//! 2nd ed., McGraw-Hill 1961, §2.1, §2.17), the shifted (LDLᵀ) form against the plain one, and +//! what buckling steps refuse. + +mod common; + +use common::*; +use fieldcraft_model::*; +use fieldcraft_solver::{Limits, solve_buckle}; + +const E: f64 = 1000.0; +const B: f64 = 1.0; +const L: f64 = 20.0; + +fn buckle(modes: u32, shift: Option) -> Procedure { + Procedure::Buckle { modes, shift } +} + +/// Engesser's critical load of a square column of side B and effective length `le`, ν = 0. +fn engesser(le: f64) -> f64 { + let (a, i, g) = (B * B, B.powi(4) / 12.0, E / 2.0); + let pe = std::f64::consts::PI.powi(2) * E * i / (le * le); + pe / (1.0 + pe / (5.0 / 6.0 * g * a)) +} + +/// A column of `length` along z, 2 × 2 × `cells` hex20, pressed by a unit pressure on its top. +fn column(length: f64, cells: usize, shift: Option) -> Grid { + let mut g = grid([2, 2, cells], ElementKind::Hex20, |x| [B * x[0], B * x[1], length * x[2]]); + g.define(iso(E, 0.0), None, Formulation::default(), buckle(3, shift)); + g.face_set("top", |x| near(x[2], 1.0)); + g.load("press", Load::Pressure { set: "top".into(), magnitude: 1.0 }); + g +} + +#[test] +fn euler_columns() { + // Pinned–pinned: half the column, the pinned end (z = 0 here is the loaded top) held + // laterally over its face (which leaves its rotation free), midspan a symmetry plane. + let mut g = grid([2, 2, 10], ElementKind::Hex20, |x| [B * x[0], B * x[1], L / 2.0 * x[2]]); + g.define(iso(E, 0.0), None, Formulation::default(), buckle(2, None)); + g.face_set("end", |x| near(x[2], 0.0)); + g.load("press", Load::Pressure { set: "end".into(), magnitude: 1.0 }); + g.node_set("pin", |x| near(x[2], 0.0)); + g.bc("pin", BoundaryCondition::Displacement { target: Target::Set("pin".into()), values: [Some(0.0), Some(0.0), None, None, None, None] }); + g.symmetry("mid", Symmetry::Z, |x| near(x[2], 1.0)); + let r = solve_buckle(&g.model, Some("Step"), &Limits::default()).unwrap(); + let want = engesser(L) / (B * B); + let got = r.modes[0].eigenvalue; + eprintln!("pinned: λ = {got:.6} (Engesser {want:.6}, {:+.3} %), second {:.6}", 100.0 * (got / want - 1.0), r.modes[1].eigenvalue); + assert!((got / want - 1.0).abs() < 5e-3); + assert!((r.modes[1].eigenvalue / got - 1.0).abs() < 1e-6, "the square section's two planes"); + assert!(r.buckling.as_ref().unwrap().negative.is_empty() || r.buckling.as_ref().unwrap().negative[0] < -got); + + // Clamped–free: the base held, the top free: le = 2L. + let mut g = column(L / 2.0, 10, None); + g.fixed("base", |x| near(x[2], 0.0)); + let r = solve_buckle(&g.model, Some("Step"), &Limits::default()).unwrap(); + let want = engesser(L) / (B * B); + let got = r.modes[0].eigenvalue; + eprintln!("clamped–free: λ = {got:.6} (Engesser {want:.6}, {:+.3} %)", 100.0 * (got / want - 1.0)); + assert!((got / want - 1.0).abs() < 5e-3); + // Scaled to a largest component of 1, which is lateral at the free top. + let big = r.modes[0].shape.iter().flatten().fold(0.0f64, |m, v| m.max(v.abs())); + assert!((big - 1.0).abs() < 1e-12); + + // Shifted between the first and the higher modes (LDLᵀ of K + σ K_σ): its inertia counts + // the two load factors below the shift, and the modes above it come out. + let mut s = column(L / 2.0, 10, Some(2.0 * got)); + s.fixed("base", |x| near(x[2], 0.0)); + let shifted = solve_buckle(&s.model, Some("Step"), &Limits::default()).unwrap(); + let info = shifted.buckling.as_ref().unwrap(); + assert_eq!(info.below_shift, Some(2)); + let third = r.modes[2].eigenvalue; + assert!((shifted.modes[0].eigenvalue / third - 1.0).abs() < 1e-7, "{} vs {third}", shifted.modes[0].eigenvalue); +} + +#[test] +fn buckling_steps_refuse_what_they_cannot_do() { + // No load: nothing to buckle. + let mut g = column(L, 4, None); + g.model.set_load_active("Step", "press", false).unwrap(); + g.fixed("base", |x| near(x[2], 0.0)); + assert!(solve_buckle(&g.model, Some("Step"), &Limits::default()).is_err()); + // Unsupported: a located singular stiffness. + let g = column(L, 4, None); + let e = solve_buckle(&g.model, Some("Step"), &Limits::default()).unwrap_err(); + assert!(e.to_string().contains("singular"), "{e}"); + // A shift that is not finite, and a step of another kind. + assert!(Procedure::Buckle { modes: 0, shift: None }.validate().is_err()); + assert!(Procedure::Buckle { modes: 1, shift: Some(f64::NAN) }.validate().is_err()); +} diff --git a/crates/ui-egui/src/results.rs b/crates/ui-egui/src/results.rs index b47ed96..30bbc41 100644 --- a/crates/ui-egui/src/results.rs +++ b/crates/ui-egui/src/results.rs @@ -64,11 +64,18 @@ pub struct Prepared { pub contour: Option>, } +/// A mode's value as the results show it: its frequency, or a buckling mode's load factor. +pub fn mode_value(buckling: bool, eigenvalue: f64, frequency: f64) -> String { + if buckling { format!("λ = {}", number(eigenvalue)) } else { format!("{} Hz", number(frequency)) } +} + impl Plot { - /// A frequency step's mode shown, and its frequency. - pub fn mode(&self) -> Option<(usize, f64)> { + /// A frequency or buckling step's mode shown, and its frequency (Hz) or load factor. + pub fn mode(&self) -> Option<(usize, String)> { match self.result.kind() { - Kind::Frequency(f) => f.modes.get(self.pick.mode.checked_sub(1)?).map(|m| (self.pick.mode, m.frequency)), + Kind::Frequency(f) => { + f.modes.get(self.pick.mode.checked_sub(1)?).map(|m| (self.pick.mode, mode_value(f.buckling.is_some(), m.eigenvalue, m.frequency))) + } Kind::Static(_) | Kind::Explicit(_) | Kind::Fluid(_) | Kind::Heat(_) => None, } } diff --git a/crates/ui-egui/src/results_panel.rs b/crates/ui-egui/src/results_panel.rs index 9ee918a..0de3138 100644 --- a/crates/ui-egui/src/results_panel.rs +++ b/crates/ui-egui/src/results_panel.rs @@ -35,8 +35,8 @@ fn legend(plot: &Plot, painter: &egui::Painter, rect: Rect, t: &Tokens) { (None, None) => "Deformed shape".into(), }; let mut step = format!("Step {}", plot.view.step); - if let Some((mode, hz)) = plot.mode() { - step.push_str(&format!(", mode {mode}: {} Hz", number(hz))); + if let Some((mode, value)) = plot.mode() { + step.push_str(&format!(", mode {mode}: {value}")); } if let Some((k, time, n)) = plot.frame() { step.push_str(&format!(", frame {k} of {}: t = {}", n.saturating_sub(1), number(time))); @@ -179,7 +179,13 @@ pub fn controls(app: &mut ModelApp, plot: &Plot, ui: &mut egui::Ui, t: &Tokens) egui::ScrollArea::vertical().id_salt("fc_results_modes").max_height(150.0).show(ui, |ui| { for (i, m) in f.modes.iter().enumerate() { let mode = i + 1; - if ui.selectable_label(plot.view.mode == mode, format!("{mode} {} Hz", number(m.frequency))).clicked() { + if ui + .selectable_label( + plot.view.mode == mode, + format!("{mode} {}", crate::results::mode_value(f.buckling.is_some(), m.eigenvalue, m.frequency)), + ) + .clicked() + { run.push(("results.modes", json!({"mode": mode}))); } } diff --git a/docs/features.md b/docs/features.md index 071152d..8d7bb7f 100644 --- a/docs/features.md +++ b/docs/features.md @@ -123,7 +123,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; nonlinear static and heat transfer in the drop-down | ready | `solve.static`, `solve.nonlinear`, `solve.heat` | +| Implicit | Static — Linear static (or coupled temperature–displacement) solve of a step with FieldCraft's own solver; nonlinear static, buckling and heat transfer in the drop-down | ready | `solve.static`, `solve.nonlinear`, `solve.buckle`, `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) | | diff --git a/docs/fem.md b/docs/fem.md index 655438a..f6d09d8 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -173,6 +173,26 @@ On a hex8 stiffness pattern of 140 000 equations (1.5 · 10⁸ entries) the LDL times the Cholesky factor, and 1.2 times on a shift with 10 600 negative eigenvalues and 7 700 2 × 2 pivots (`ldlt::tests::speed_against_cholesky`). +## Linear buckling + +A step `{"type": "buckle", "modes": 4}` (`Procedure::Buckle`, solved by `solve.buckle`, +`fieldcraft_solver::buckle`) solves its loads and prescribed motions as a linear static preload, +forms the initial-stress stiffness K_σ = ∫ (∇N_a · σ ∇N_b) I dV of every solid and 2D element +from the preload's stresses at the stiffness points (`Solid::geometric_stiffness`, +`Plane::geometric_stiffness`, with the axisymmetric hoop term), and finds the load factors λ +with (K + λ K_σ) φ = 0 by block Lanczos in the K inner product (Grimes, Lewis & Simon 1994): +without a shift on K⁻¹(−K_σ) (θ = 1/λ, Cholesky of K), the lowest positive λ; with `"shift": σ` +on the buckling transformation (K + σ K_σ)⁻¹ K (θ = λ/(λ − σ), LDLᵀ), the λ above σ, and the +LDLᵀ inertia counts the load factors between 0 and σ, reported. Negative load factors (the +reversed load buckles) are reported, not returned. Modes are scaled to a largest component of 1 +and display and probe like frequency modes (their value is λ; `FrequencyResult::buckling`). +Boundary conditions hold their directions at zero in the eigenproblem; pressures are dead loads +(no load stiffness); trusses, shells and beams are refused (no K_σ yet). Decks: Abaqus +`*BUCKLE, EIGENSOLVER=LANCZOS` (refused after a general step, whose loaded state Abaqus would +buckle from, and with a shift); Nastran SOL 105 and LS-DYNA refuse it by name. Verified by Euler +columns of hex20 (pinned and clamped–free, 0.18 % from Engesser's shear-corrected load) and the +shifted form against the plain one. + ## Nonlinear static analysis (finite strain) A step `{"type": "nonlinear_static", …}` (`Procedure::NonlinearStatic`, solved by diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index 5ff2ecb..e59ea68 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -1246,6 +1246,11 @@ update are written from (cited again in the module docs): 34(1):82–109, 1992. - Iterative refinement: N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM 2002, §12.2. +- Linear buckling and the initial-stress stiffness: K.-J. Bathe, *Finite Element Procedures*, + 2nd ed., 2014, §6.8.2; the buckling spectral transformations in the K inner product: R. G. + Grimes, J. G. Lewis, H. D. Simon, SIAM J. Matrix Anal. Appl. 15(1):228–272, 1994, §3; Euler's + columns with Engesser's shear correction: S. P. Timoshenko, J. M. Gere, *Theory of Elastic + Stability*, 2nd ed., McGraw-Hill 1961, §2.1, §2.17. ## Part Q: Contact and friction (explicit), as implemented