diff --git a/AGENTS.md b/AGENTS.md index d054069..525448b 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -75,9 +75,9 @@ People trust ModelCraft with their models; a crash loses their work. **This outr JSON `run`) + tests. The UI, scripts, CLI and control channel all reach the same commands. - **Programmatic calls never open dialogs.** - **Layering** is enforced by `cargo xtask layers` (`docs/architecture.md`): L0 `model` → L1 - `geometry`, `mesh`, `fem`, `decks`, `render` → L2 solvers (`solver`, `explicit`, `lbm`) → L3 - `fsi` → then `doc`, `io`, `engine`, UI and apps above (the layer table is renumbered when the - first solver crate lands). Nothing below L5 depends on egui/eframe/winit/wgpu/rfd. **Only + `geometry`, `mesh`, `fem`, `decks`, `render` → L2 solvers (`solver`; `explicit`, `lbm` later) → + L3 `fsi` (later) → L4 `doc` → L5 `io` → L6 `engine` → L7 `ui-egui`, `mcp` → apps. Nothing below + L7 depends on egui/eframe/winit/wgpu/rfd. **Only `crates/geometry` may use `solidcraft-*`.** - **Determinism.** The same input gives the same mesh, deck and CPU solver result, bit for bit, on every platform: no hash-order iteration into output, no thread-count-dependent results. GPU diff --git a/Cargo.lock b/Cargo.lock index 17865ae..b5f52d9 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1860,12 +1860,21 @@ dependencies = [ "modelcraft-io", "modelcraft-mesh", "modelcraft-model", + "modelcraft-solver", "serde", "serde_json", "thiserror 2.0.21", "zip", ] +[[package]] +name = "modelcraft-fem" +version = "0.1.0" +dependencies = [ + "modelcraft-model", + "thiserror 2.0.21", +] + [[package]] name = "modelcraft-geometry" version = "0.1.0" @@ -1926,6 +1935,16 @@ dependencies = [ "serde_json", ] +[[package]] +name = "modelcraft-solver" +version = "0.1.0" +dependencies = [ + "modelcraft-fem", + "modelcraft-model", + "serde", + "thiserror 2.0.21", +] + [[package]] name = "modelcraft-testkit" version = "0.1.0" diff --git a/Cargo.toml b/Cargo.toml index 8c208e0..a3e232f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,6 +17,8 @@ modelcraft-geometry = { path = "crates/geometry" } modelcraft-mesh = { path = "crates/mesh" } modelcraft-decks = { path = "crates/decks" } modelcraft-io = { path = "crates/io" } +modelcraft-fem = { path = "crates/fem" } +modelcraft-solver = { path = "crates/solver" } modelcraft-engine = { path = "crates/engine" } modelcraft-testkit = { path = "crates/testkit" } modelcraft-render = { path = "crates/render" } diff --git a/ROADMAP.md b/ROADMAP.md index 690cc53..40316c2 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -22,7 +22,7 @@ hex meshes today and it gives the best elements where stresses are read. | M1 Geometry in and surface preparation | Open `.solidcraft` designs with persistent face names, assemblies placed (done); face-tagged 3MF (done); test corpus, `cargo xtask corpus` (done); exact predicates (done); quality metrics and certified validity gate (done); weld, closed/manifold/orientation and self-intersection checks, feature edges (`modelcraft-mesh::surface`, `mesh.surface.check`, done) | in progress | | M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `modelcraft_mesh::size::SizeField` trait and edge-division helpers are in); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress | | M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.modelcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress | -| M4 Own solver and first results | Linear static and modal solver (tet4/10, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions; deformed shape, contours, mode shapes, probes, reaction sums | planned | +| M4 Own solver and first results | Linear static and modal solver (tet4/10, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests); deformed shape, contours, mode shapes, probes, reaction sums | started | | M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | | M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, springs, point masses; amplitudes, multiple steps; nonlinear static steps in decks | planned | diff --git a/crates/engine/Cargo.toml b/crates/engine/Cargo.toml index ad03d96..2c4bfb5 100644 --- a/crates/engine/Cargo.toml +++ b/crates/engine/Cargo.toml @@ -14,6 +14,7 @@ modelcraft-geometry = { workspace = true } modelcraft-mesh = { workspace = true } modelcraft-io = { workspace = true } modelcraft-decks = { workspace = true } +modelcraft-solver = { workspace = true } serde = { workspace = true } serde_json = { workspace = true } thiserror = { workspace = true } diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index 3ec2311..b32489e 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -7,6 +7,7 @@ pub(crate) mod file; mod geometry; mod mesh; mod model; +mod solve; mod tags; use modelcraft_model::MAX_LABEL; @@ -68,6 +69,9 @@ pub static ALL: &[CommandSpec] = &[ edit::UNDO, edit::REDO, deck::WRITE, + solve::STATIC, + solve::FREQUENCY, + solve::PROBE, ]; /// A required label parameter (1..=u32::MAX; the model checks its own range). diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs new file mode 100644 index 0000000..c0ddf0c --- /dev/null +++ b/crates/engine/src/cmd/solve.rs @@ -0,0 +1,268 @@ +//! `solve.*`: run the model's steps with the implicit solver, keep the results in the session, +//! and read them back. + +use modelcraft_solver::{FrequencyOptions, Limits, MassKind, StepResult, solve_frequency, solve_static}; +use serde_json::{Value, json}; + +use super::{labels, opt}; +use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; + +pub const STATIC: CommandSpec = CommandSpec { + id: "solve.static", + label: "Solve Static Step", + params: r#"{"step": "Load", "max_dofs": 6000000, "max_factor_entries": 400000000} — solves a linear static step (the only static step when `step` is left out) and keeps the result in the session; the caps may only be lowered"#, + kind: CommandKind::Model, + run: run_static, +}; + +pub const FREQUENCY: CommandSpec = CommandSpec { + id: "solve.frequency", + label: "Solve Frequency Step", + params: r#"{"step": "Modes", "mass": "consistent|lumped", "shift": 0.0, "max_dofs": 6000000, "max_factor_entries": 400000000, "max_iterations": 3000} — natural frequencies and mode shapes of a frequency step (its modes and range); `max_iterations` caps the Lanczos vectors; the result stays in the session"#, + kind: CommandKind::Model, + run: run_frequency, +}; + +pub const PROBE: CommandSpec = CommandSpec { + id: "solve.probe", + label: "Probe Results", + params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1} — the stored result at nodes: displacement, reaction, stress and von Mises (static), or the shape of mode `mode` (frequency, from 1)"#, + kind: CommandKind::Model, + run: probe, +}; + +/// Most nodes one probe lists. +const MAX_PROBE: usize = 10_000; + +/// An optional cap that may only lower the default. +fn cap(p: &Value, key: &str, default: usize) -> Result { + match p.get(key) { + None => Ok(default), + Some(v) => v + .as_u64() + .and_then(|n| usize::try_from(n).ok()) + .filter(|n| (1..=default).contains(n)) + .ok_or_else(|| EngineError::Params(format!("`{key}` must be a whole number from 1 to {default}"))), + } +} + +fn limits(p: &Value) -> Result { + Ok(Limits { + 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)?, + }) +} + +fn largest(values: impl Iterator) -> Value { + match values.fold(None, |best: Option<(u32, f64)>, (l, v)| if best.is_none_or(|(_, b)| v > b) { Some((l, v)) } else { best }) { + Some((node, value)) => json!({"node": node, "value": value}), + None => Value::Null, + } +} + +fn run_static(s: &mut Session, p: &Value) -> Result { + let limits = limits(p)?; + let step = opt::(p, "step")?; + let r = solve_static(&s.model, step.as_deref(), &limits)?; + let magnitude = |u: &[f64; 3]| (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]).sqrt(); + let summary = json!({ + "step": r.step, + "dofs": r.stats.dofs, + "equations": r.stats.equations, + "factor_entries": r.stats.factor_entries, + "max_displacement": largest(r.nodes.iter().zip(&r.displacement).map(|(l, u)| (*l, 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, + "warnings": r.warnings, + }); + s.store_result(StepResult::Static(r)); + Ok(summary) +} + +fn run_frequency(s: &mut Session, p: &Value) -> Result { + let limits = limits(p)?; + let step = opt::(p, "step")?; + let mass = opt::(p, "mass")?.unwrap_or_default(); + let shift = opt::(p, "shift")?; + let r = solve_frequency(&s.model, step.as_deref(), &FrequencyOptions { mass, shift }, &limits)?; + let summary = json!({ + "step": r.step, + "dofs": r.stats.dofs, + "equations": r.stats.equations, + "lanczos_vectors": r.stats.lanczos_vectors, + "total_mass": r.total_mass, + "modes": r.modes.iter().enumerate().map(|(i, m)| json!({ + "mode": i + 1, + "frequency": m.frequency, + "eigenvalue": m.eigenvalue, + "effective_mass": m.effective_mass, + })).collect::>(), + "warnings": r.warnings, + }); + s.store_result(StepResult::Frequency(r)); + Ok(summary) +} + +fn probe(s: &mut Session, p: &Value) -> Result { + let step = opt::(p, "step")?.ok_or_else(|| EngineError::Params("`step` is required".into()))?; + let result = s.result(&step).ok_or_else(|| EngineError::Params(format!("step `{step}` has no result; solve it first")))?; + let nodes = labels(p, "nodes", MAX_PROBE)?; + let index = |label: u32| { + result.nodes().binary_search(&label).map_err(|_| EngineError::Params(format!("node {label} has no result (it belongs to no solved element)"))) + }; + let out: Vec = match result.as_ref() { + StepResult::Static(r) => nodes + .iter() + .map(|&l| { + let i = index(l)?; + Ok(json!({ + "node": l, + "displacement": r.displacement.get(i), + "reaction": r.reaction.get(i), + "stress": r.stress.get(i), + "von_mises": r.von_mises.get(i), + })) + }) + .collect::>()?, + StepResult::Frequency(r) => { + let mode = opt::(p, "mode")?.unwrap_or(1); + let m = mode + .checked_sub(1) + .and_then(|k| r.modes.get(k)) + .ok_or_else(|| EngineError::Params(format!("`mode` must be from 1 to {}", r.modes.len())))?; + nodes.iter().map(|&l| Ok(json!({"node": l, "shape": m.shape.get(index(l)?)}))).collect::>()? + } + }; + Ok(json!({"step": step, "nodes": out})) +} + +#[cfg(test)] +pub(crate) mod tests { + use super::*; + + /// A cantilever 10 × 1 × 1 of ten hex8 elements with incompatible modes, clamped at x = 0, + /// a total force of 1 down on the free end; a static and a frequency step. + pub(crate) fn cantilever(s: &mut Session) { + let mut script = Vec::new(); + let label = |i: usize, j: usize, k: usize| 1 + i + 11 * (j + 2 * k); + for k in 0..2 { + for j in 0..2 { + for i in 0..11 { + script.push(json!({"command": "mesh.node.add", "params": {"label": label(i, j, k), "xyz": [i, j, k]}})); + } + } + } + for i in 0..10 { + let nodes = [ + label(i, 0, 0), + label(i + 1, 0, 0), + label(i + 1, 1, 0), + label(i, 1, 0), + label(i, 0, 1), + label(i + 1, 0, 1), + label(i + 1, 1, 1), + label(i, 1, 1), + ]; + script.push(json!({"command": "mesh.element.add", "params": {"label": i + 1, "kind": "hex8", "nodes": nodes}})); + } + let root: Vec = (0..4).map(|n| label(0, n % 2, n / 2)).collect(); + let tip: Vec = (0..4).map(|n| label(10, n % 2, n / 2)).collect(); + script.extend([ + json!({"command": "model.set.create", "params": {"name": "all", "elements": (1..=10).collect::>()}}), + json!({"command": "model.set.create", "params": {"name": "root", "nodes": root}}), + json!({"command": "model.set.create", "params": {"name": "tip", "nodes": tip}}), + json!({"command": "model.material.create", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 1000, "poisson": 0}, "density": 1}}), + json!({"command": "model.section.create", "params": {"name": "bar", "set": "all", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"incompatible_modes": true}}}), + json!({"command": "model.boundary_condition.create", "params": {"name": "clamp", "type": "fixed", "target": {"set": "root"}}}), + json!({"command": "model.load.create", "params": {"name": "push", "type": "force", "target": {"set": "tip"}, "vector": [0, 0, -1]}}), + json!({"command": "model.step.create", "params": {"name": "Bend", "procedure": {"type": "static"}, "boundary_conditions": ["clamp"], "loads": ["push"]}}), + json!({"command": "model.step.create", "params": {"name": "Modes", "procedure": {"type": "frequency", "modes": 4}, "boundary_conditions": ["clamp"]}}), + ]); + s.run_script(&Value::Array(script)).unwrap(); + } + + #[test] + fn solve_keep_and_probe_results() { + let mut s = Session::default(); + cantilever(&mut s); + let journal = s.journal.len(); + let r = s.execute("solve.static", &json!({})).unwrap(); + assert_eq!(r["step"], "Bend"); + // Euler–Bernoulli P L³ / (3 E I) = 4, plus shear; the clamp and one element through the + // depth keep it within a few percent. + let d = r["max_displacement"]["value"].as_f64().unwrap(); + assert!((d - 4.024).abs() / 4.024 < 0.03, "{r}"); + assert!((r["reaction_total"][2].as_f64().unwrap() - 1.0).abs() < 1e-9, "{r}"); + assert_eq!(s.journal.len(), journal, "solving does not change the model"); + let probe = s.execute("solve.probe", &json!({"step": "Bend", "nodes": [11, 1]})).unwrap(); + assert!(probe["nodes"][0]["displacement"][2].as_f64().unwrap() < -3.8, "{probe}"); + assert!(probe["nodes"][1]["reaction"][2].as_f64().unwrap() != 0.0, "{probe}"); + + let f = s.execute("solve.frequency", &json!({"mass": "lumped"})).unwrap(); + assert_eq!(f["modes"].as_array().unwrap().len(), 4, "{f}"); + let shape = s.execute("solve.probe", &json!({"step": "Modes", "nodes": [44], "mode": 2})).unwrap(); + assert!(shape["nodes"][0]["shape"].is_array(), "{shape}"); + assert_eq!(s.solved_steps().count(), 2); + + // Results describe the model they were solved for: an edit clears them, a failed command + // keeps them. + assert!(s.execute("model.material.edit", &json!({"name": "Mat", "elastic": {"type": "isotropic", "young": -1}})).is_err()); + assert_eq!(s.solved_steps().count(), 2); + s.execute("model.material.edit", &json!({"name": "Mat", "elastic": {"type": "isotropic", "young": 2000, "poisson": 0}, "density": 1})) + .unwrap(); + assert_eq!(s.solved_steps().count(), 0); + // Undo returns to the solved state, and its results are current again. + s.execute("edit.undo", &json!({})).unwrap(); + assert_eq!(s.solved_steps().count(), 2); + s.execute("edit.redo", &json!({})).unwrap(); + assert!(s.execute("solve.probe", &json!({"step": "Bend", "nodes": [1]})).is_err()); + } + + #[test] + fn solve_hostile_params_never_panic() { + let mut s = Session::default(); + cantilever(&mut s); + s.execute("solve.static", &json!({})).unwrap(); + s.execute("solve.frequency", &json!({})).unwrap(); + let hostile = [ + json!(null), + json!({"step": 7}), + json!({"step": "Nope"}), + json!({"step": "Modes"}), + json!({"step": "Bend", "max_dofs": 0}), + json!({"step": "Bend", "max_dofs": 1e300}), + json!({"step": "Bend", "max_dofs": 3}), + json!({"step": "Bend", "max_factor_entries": -1}), + json!({"step": "Bend", "max_factor_entries": 10}), + json!({"step": "Modes", "max_iterations": 1}), + json!({"step": "Modes", "max_iterations": 7}), + json!({"step": "Modes", "mass": "heavy"}), + json!({"step": "Modes", "shift": 1e308}), + json!({"step": "Modes", "shift": -1e308}), + json!({"step": "Modes", "shift": 5000.0}), + json!({"step": "Modes", "shift": "x"}), + json!({"step": "Modes", "nodes": [1], "mode": 0}), + json!({"step": "Modes", "nodes": [1], "mode": 1e300}), + json!({"step": "Bend", "nodes": [0, 99999999, 4294967296u64]}), + json!({"step": "Bend", "nodes": vec![1; 20_000]}), + ]; + for id in ["solve.static", "solve.frequency", "solve.probe"] { + for p in &hostile { + let mut t = s.clone(); + let before = t.model.clone(); + let results = t.solved_steps().count(); + if t.execute(id, p).is_err() { + assert_eq!(t.model, before); + assert_eq!(t.solved_steps().count(), results, "{id} {p}"); + } + } + } + // An unsupported model (a coupling, a reference point load) is an error, not a panic. + let mut defined = Session::default(); + crate::cmd::define::tests::defined_bar(&mut defined); + let e = defined.execute("solve.static", &json!({"step": "Bend"})).unwrap_err().to_string(); + assert!(e.contains("not supported yet"), "{e}"); + } +} diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 96c919c..1dab7e5 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -17,11 +17,13 @@ pub mod cmd; pub mod history; pub mod sample; +use std::collections::BTreeMap; use std::sync::Arc; use modelcraft_geometry::Design; use modelcraft_geometry::tagged_mesh::TaggedSurface; use modelcraft_model::Model; +use modelcraft_solver::StepResult; use serde_json::Value; use crate::history::History; @@ -42,6 +44,8 @@ pub enum EngineError { File(#[from] modelcraft_io::Error), #[error(transparent)] Deck(#[from] modelcraft_decks::Error), + #[error(transparent)] + Solver(#[from] modelcraft_solver::Error), /// A command panicked; the model is unchanged. #[error("internal error in {command}: {message}")] Internal { command: String, message: String }, @@ -110,6 +114,9 @@ pub struct Session { /// The revision last saved to (or opened from) `path`; `None` when the content was never /// saved, where only the empty start (revision 0) is clean. saved_revision: Option, + /// The latest result of each solved step, by step name, with the revision it was solved at; + /// [`Session::result`] hands out only those that describe the current state. + results: BTreeMap)>, } impl Session { @@ -163,6 +170,22 @@ impl Session { self.saved_revision.unwrap_or(0) != self.revision } + /// The stored result of a step, if it was solved in the current state (an edit makes it stale; + /// undoing back to the solved state makes it current again). + pub fn result(&self, step: &str) -> Option<&Arc> { + self.results.get(step).filter(|(revision, _)| *revision == self.revision).map(|(_, r)| r) + } + + /// The steps with a current result. + pub fn solved_steps(&self) -> impl Iterator { + self.results.iter().filter(|(_, (revision, _))| *revision == self.revision).map(|(name, _)| name.as_str()) + } + + /// Keep a step's result for the current state. + pub(crate) fn store_result(&mut self, result: StepResult) { + self.results.insert(result.step().to_string(), (self.revision, Arc::new(result))); + } + /// The current state is what `path` now holds. pub(crate) fn mark_saved(&mut self, path: String) { self.path = Some(path); diff --git a/crates/fem/Cargo.toml b/crates/fem/Cargo.toml new file mode 100644 index 0000000..61a5d0f --- /dev/null +++ b/crates/fem/Cargo.toml @@ -0,0 +1,13 @@ +[package] +name = "modelcraft-fem" +description = "ModelCraft's finite element kernels: shape functions, quadrature, strain-displacement matrices, element stiffness and mass, elastic constitutive matrices" +version.workspace = true +edition.workspace = true +license.workspace = true + +[lints] +workspace = true + +[dependencies] +modelcraft-model = { workspace = true } +thiserror = { workspace = true } diff --git a/crates/fem/src/dense.rs b/crates/fem/src/dense.rs new file mode 100644 index 0000000..307161c --- /dev/null +++ b/crates/fem/src/dense.rs @@ -0,0 +1,138 @@ +//! Small dense matrices: element matrices and the little systems inside an element (static +//! condensation, extrapolation fits). + +use crate::{Error, Result}; + +/// A square matrix, row-major. +#[derive(Debug, Clone, PartialEq)] +pub struct Matrix { + pub n: usize, + pub data: Vec, +} + +impl Matrix { + pub fn zeros(n: usize) -> Matrix { + Matrix { n, data: vec![0.0; n * n] } + } + + pub fn get(&self, i: usize, j: usize) -> f64 { + self.data.get(i * self.n + j).copied().unwrap_or(0.0) + } + + pub fn add(&mut self, i: usize, j: usize, v: f64) { + if let Some(x) = self.data.get_mut(i * self.n + j) { + *x += v; + } + } + + /// y = A x (x shorter than n is padded with zeros). + pub fn mul_vec(&self, x: &[f64]) -> Vec { + self.data.chunks_exact(self.n.max(1)).map(|row| row.iter().zip(x).map(|(a, b)| a * b).sum()).collect() + } +} + +/// The Cholesky factor of a small symmetric positive definite matrix (lower triangle, row-major). +#[derive(Debug, Clone)] +pub struct Cholesky { + n: usize, + l: Vec, +} + +impl Cholesky { + /// Factor `a` (only its lower triangle is read). Fails when a pivot is not positive relative + /// to the matrix's diagonal. + pub fn new(a: &Matrix) -> Result { + let n = a.n; + let mut l = vec![0.0; n * n]; + let scale = (0..n).map(|i| a.get(i, i).abs()).fold(0.0, f64::max); + for j in 0..n { + let mut d = a.get(j, j); + for k in 0..j { + let v = l[j * n + k]; + d -= v * v; + } + if !(d > 1e-13 * scale) { + return Err(Error::Singular(format!("a {n} × {n} element matrix is not positive definite (pivot {j})"))); + } + let djj = d.sqrt(); + l[j * n + j] = djj; + for i in j + 1..n { + let mut s = a.get(i, j); + for k in 0..j { + s -= l[i * n + k] * l[j * n + k]; + } + l[i * n + j] = s / djj; + } + } + Ok(Cholesky { n, l }) + } + + /// Solve A x = b in place. + pub fn solve(&self, b: &mut [f64]) { + let n = self.n.min(b.len()); + for i in 0..n { + let mut s = b[i]; + for k in 0..i { + s -= self.l[i * self.n + k] * b[k]; + } + b[i] = s / self.l[i * self.n + i]; + } + for i in (0..n).rev() { + let mut s = b[i]; + for k in i + 1..n { + s -= self.l[k * self.n + i] * b[k]; + } + b[i] = s / self.l[i * self.n + i]; + } + } +} + +/// The inverse and determinant of a 3 × 3 matrix; `None` when it is singular or not finite. +pub fn inverse3(m: &[[f64; 3]; 3]) -> Option<([[f64; 3]; 3], f64)> { + let det = det3(m); + if !(det.abs() > 0.0 && det.is_finite()) { + return None; + } + let mut inv = [[0.0; 3]; 3]; + for (i, row) in inv.iter_mut().enumerate() { + for (j, v) in row.iter_mut().enumerate() { + let (j1, j2) = ((j + 1) % 3, (j + 2) % 3); + let (i1, i2) = ((i + 1) % 3, (i + 2) % 3); + // inv[i][j] = cofactor[j][i] / det. + *v = (m[j1][i1] * m[j2][i2] - m[j1][i2] * m[j2][i1]) / det; + } + } + Some((inv, det)) +} + +pub fn det3(m: &[[f64; 3]; 3]) -> f64 { + m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn cholesky_solves_and_rejects_indefinite() { + let mut a = Matrix::zeros(3); + a.data = vec![4.0, 2.0, 0.6, 2.0, 5.0, 1.0, 0.6, 1.0, 3.0]; + let x = [1.0, -2.0, 0.5]; + let mut b = a.mul_vec(&x); + Cholesky::new(&a).unwrap().solve(&mut b); + assert!(b.iter().zip(x).all(|(p, q)| (p - q).abs() < 1e-14), "{b:?}"); + a.data[8] = -1.0; + assert!(Cholesky::new(&a).is_err()); + let m = [[2.0, 0.5, 0.0], [0.1, 3.0, 0.2], [0.0, 0.4, 1.0]]; + let (inv, det) = inverse3(&m).unwrap(); + assert!((det - det3(&m)).abs() < 1e-15); + for i in 0..3 { + for j in 0..3 { + let p: f64 = (0..3).map(|k| m[i][k] * inv[k][j]).sum(); + assert!((p - if i == j { 1.0 } else { 0.0 }).abs() < 1e-14); + } + } + assert!(inverse3(&[[1.0, 2.0, 3.0], [2.0, 4.0, 6.0], [0.0, 0.0, 1.0]]).is_none()); + } +} diff --git a/crates/fem/src/element.rs b/crates/fem/src/element.rs new file mode 100644 index 0000000..eefe788 --- /dev/null +++ b/crates/fem/src/element.rs @@ -0,0 +1,610 @@ +//! Solid continuum elements: the element table, stiffness, mass, body forces, stresses at the +//! integration points and their extrapolation to the nodes. +//! +//! The isoparametric displacement formulation (Zienkiewicz, Taylor & Zhu 2013, ch. 6; Hughes +//! 2000, ch. 3): K = ∫ Bᵀ D B dV and M = ∫ ρ Nᵀ N dV, integrated by Gauss quadrature on the +//! reference element through the Jacobian of the map from natural to physical coordinates. +//! +//! Formulations: +//! - **full integration**: the rule that integrates the stiffness of an undistorted element +//! exactly (1 point for tet4, 4 for tet10, 2 × 2 × 2 for hex8, 3 × 3 × 3 for hex20); +//! - **reduced integration** for hex20 (2 × 2 × 2), which avoids locking at the cost of one +//! spurious mode that a mesh of more than one element restrains. Reduced hex8 needs hourglass +//! control and is not offered here yet (it arrives with the explicit solver, Flanagan & +//! Belytschko 1981); +//! - **incompatible modes** for hex8 (E. L. Wilson, R. L. Taylor, W. P. Doherty, J. Ghaboussi, +//! *Incompatible displacement models*, in Numerical and Computer Methods in Structural +//! Mechanics, Academic Press 1973, 43–57): nine internal displacement modes (1 − ξ²), +//! (1 − η²), (1 − ζ²) per direction, condensed out of the element. Their derivatives use the +//! Jacobian at the element centre scaled by det J₀ / det J, the correction of R. L. Taylor, +//! P. J. Beresford, E. L. Wilson, *A non-conforming element for stress analysis*, IJNME +//! 10(6):1211–1219, 1976, so the element passes the patch test when distorted. It bends +//! correctly with one element through the thickness. +//! +//! Lumped mass is the diagonal scaling of Hinton, Rock and Zienkiewicz (E. Hinton, T. Rock, +//! O. C. Zienkiewicz, *A note on mass lumping and related processes in the finite element +//! method*, Earthquake Eng. Struct. Dyn. 4(3):245–249, 1976): the consistent diagonal scaled to +//! the element's mass, positive for quadratic elements too. +//! +//! Nodal stresses extrapolate the integration-point values with the element's linear corner +//! functions fitted by least squares (E. Hinton, J. S. Campbell, *Local and global smoothing of +//! discontinuous finite element functions using a least squares method*, IJNME 8(3):461–480, +//! 1974); mid-side nodes take the mean of their edge's corners. + +use modelcraft_model::ElementKind; + +use crate::dense::{Cholesky, Matrix, inverse3}; +use crate::material::{Elasticity, stress}; +use crate::quadrature::{Point, Rule}; +use crate::shape::{self, MAX_NODES, Shape}; +use crate::{Error, Result}; + +/// The reference shape of an element type. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Family { + Tetrahedron, + Hexahedron, +} + +/// One row of the element table: everything the kernels need to know about an element kind. +/// Adding an element (wedge6, shells, beams) is adding a row and its shape functions. +#[derive(Debug)] +pub struct ElementType { + pub kind: ElementKind, + pub family: Family, + pub nodes: usize, + /// The corner nodes come first in the node list. + pub corners: usize, + pub shape: fn([f64; 3]) -> Shape, + /// Full integration of the stiffness. + pub full: Rule, + /// Reduced integration, where the element is stable with it. + pub reduced: Option, + /// Exact for the consistent mass of an undistorted element. + pub mass: Rule, + /// Whether the element has an incompatible-modes formulation. + pub incompatible_modes: bool, +} + +/// The element table. +pub static ELEMENT_TYPES: [ElementType; 4] = [ + ElementType { + kind: ElementKind::Tet4, + family: Family::Tetrahedron, + nodes: 4, + corners: 4, + shape: shape::tet4, + full: Rule::TetCentroid, + reduced: None, + mass: Rule::Tet4, + incompatible_modes: false, + }, + ElementType { + kind: ElementKind::Tet10, + family: Family::Tetrahedron, + nodes: 10, + corners: 4, + shape: shape::tet10, + full: Rule::Tet4, + reduced: None, + mass: Rule::TetConical(4), + incompatible_modes: false, + }, + ElementType { + kind: ElementKind::Hex8, + family: Family::Hexahedron, + nodes: 8, + corners: 8, + shape: shape::hex8, + full: Rule::Hex(2), + reduced: None, + mass: Rule::Hex(2), + incompatible_modes: true, + }, + ElementType { + kind: ElementKind::Hex20, + family: Family::Hexahedron, + nodes: 20, + corners: 8, + shape: shape::hex20, + full: Rule::Hex(3), + reduced: Some(Rule::Hex(2)), + mass: Rule::Hex(3), + incompatible_modes: false, + }, +]; + +/// The table row of an element kind. +pub fn element_type(kind: ElementKind) -> Result<&'static ElementType> { + ELEMENT_TYPES + .iter() + .find(|t| t.kind == kind) + .ok_or_else(|| Error::NotSupported(format!("{} elements (the solver has tet4, tet10, hex8 and hex20 so far)", kind.name()))) +} + +/// How a solid element is formulated. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Formulation { + Full, + Reduced, + IncompatibleModes, +} + +impl Formulation { + /// The formulation a model section asks for, if this element type has it. + pub fn from_section(kind: ElementKind, f: &modelcraft_model::Formulation) -> Result { + let ty = element_type(kind)?; + if f.hybrid { + return Err(Error::NotSupported("hybrid (mixed displacement–pressure) elements".into())); + } + if f.incompatible_modes { + return if ty.incompatible_modes { + Ok(Formulation::IncompatibleModes) + } else { + Err(Error::Invalid(format!("incompatible modes apply to hex8 elements, not {}", kind.name()))) + }; + } + match f.integration { + modelcraft_model::Integration::Full => Ok(Formulation::Full), + modelcraft_model::Integration::Reduced if ty.reduced.is_some() => Ok(Formulation::Reduced), + modelcraft_model::Integration::Reduced => Err(Error::NotSupported(format!( + "reduced integration of {} elements (hex8 needs hourglass control; use full integration or incompatible modes)", + kind.name() + ))), + } + } +} + +/// Strain, stress and position at one integration point. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct PointStress { + pub xyz: [f64; 3], + pub strain: [f64; 6], + pub stress: [f64; 6], +} + +/// Shape function gradients in physical coordinates at one point. +struct Gradients { + shape: Shape, + dx: [[f64; 3]; MAX_NODES], + det: f64, +} + +/// A solid element: its type, its node coordinates and its formulation. +#[derive(Debug, Clone, Copy)] +pub struct Solid<'a> { + pub ty: &'static ElementType, + xyz: &'a [[f64; 3]], + formulation: Formulation, +} + +/// The stiffness blocks of the incompatible-modes element. +struct IncompatibleParts { + kuu: Matrix, + kua: Vec, + kaa: Matrix, + j0: ([[f64; 3]; 3], f64), +} + +/// The incompatible modes' count of internal degrees of freedom. +const INCOMPATIBLE_DOFS: usize = 9; + +impl<'a> Solid<'a> { + pub fn new(kind: ElementKind, xyz: &'a [[f64; 3]], formulation: Formulation) -> Result> { + let ty = element_type(kind)?; + if xyz.len() != ty.nodes { + return Err(Error::Invalid(format!("a {} element has {} nodes, got {}", kind.name(), ty.nodes, xyz.len()))); + } + modelcraft_model::check_coordinates(xyz.as_flattened()).map_err(|e| Error::Invalid(e.to_string()))?; + let ok = match formulation { + Formulation::Full => true, + Formulation::Reduced => ty.reduced.is_some(), + Formulation::IncompatibleModes => ty.incompatible_modes, + }; + if !ok { + return Err(Error::NotSupported(format!("{formulation:?} formulation of {} elements", kind.name()))); + } + Ok(Solid { ty, xyz, formulation }) + } + + /// The stiffness integration rule. + pub fn rule(&self) -> Rule { + match (self.formulation, self.ty.reduced) { + (Formulation::Reduced, Some(r)) => r, + _ => self.ty.full, + } + } + + /// The stiffness integration points (where stresses are reported). + pub fn integration_points(&self) -> Vec { + self.rule().points() + } + + /// The Jacobian matrix J[a][b] = ∂x_b/∂ξ_a. + fn jacobian(&self, s: &Shape) -> [[f64; 3]; 3] { + let mut j = [[0.0; 3]; 3]; + for (d, x) in s.derivatives().iter().zip(self.xyz) { + for a in 0..3 { + for b in 0..3 { + j[a][b] += d[a] * x[b]; + } + } + } + j + } + + fn gradients(&self, xi: [f64; 3]) -> Result { + let shape = (self.ty.shape)(xi); + let j = self.jacobian(&shape); + let (inv, det) = inverse3(&j).ok_or_else(|| self.degenerate())?; + let norms: f64 = j.iter().map(|r| (r[0] * r[0] + r[1] * r[1] + r[2] * r[2]).sqrt()).product(); + if !(det > 1e-12 * norms) { + return Err(self.degenerate()); + } + let mut dx = [[0.0; 3]; MAX_NODES]; + for (g, d) in dx.iter_mut().zip(shape.derivatives()) { + for b in 0..3 { + g[b] = inv[b][0] * d[0] + inv[b][1] * d[1] + inv[b][2] * d[2]; + } + } + Ok(Gradients { shape, dx, det }) + } + + fn degenerate(&self) -> Error { + Error::Degenerate(format!("a {} element is inverted or degenerate (its Jacobian is not positive)", self.ty.kind.name())) + } + + /// Position of a natural point. + pub fn position(&self, xi: [f64; 3]) -> [f64; 3] { + let s = (self.ty.shape)(xi); + let mut p = [0.0; 3]; + for (n, x) in s.values().iter().zip(self.xyz) { + for k in 0..3 { + p[k] += n * x[k]; + } + } + p + } + + /// The element's volume. + pub fn volume(&self) -> Result { + let mut v = 0.0; + for p in self.ty.mass.points() { + v += self.gradients(p.xi)?.det * p.weight; + } + Ok(v) + } + + /// The incompatible modes' gradients at a point: det J₀/det J · J₀⁻¹ ∂P/∂ξ. + fn incompatible_gradients(&self, xi: [f64; 3], det: f64, j0: &([[f64; 3]; 3], f64)) -> [[f64; 3]; 3] { + let (inv0, det0) = j0; + let mut g = [[0.0; 3]; 3]; + for (k, gk) in g.iter_mut().enumerate() { + // ∂(1 − ξ_k²)/∂ξ_k = −2 ξ_k, other derivatives zero. + let dp = -2.0 * xi[k]; + for b in 0..3 { + gk[b] = det0 / det * inv0[b][k] * dp; + } + } + g + } + + fn centre_jacobian(&self) -> Result<([[f64; 3]; 3], f64)> { + let s = (self.ty.shape)([0.0; 3]); + inverse3(&self.jacobian(&s)).filter(|(_, d)| *d > 0.0).ok_or_else(|| self.degenerate()) + } + + /// The parts of the stiffness with the incompatible modes: K_uu (when `with_uu`), K_ua + /// (3n × 9), K_aa, and the centre Jacobian they were built with. + fn incompatible_parts(&self, d: &Elasticity, with_uu: bool) -> Result { + let n = 3 * self.ty.nodes; + let j0 = self.centre_jacobian()?; + let mut kuu = Matrix::zeros(if with_uu { n } else { 0 }); + let mut kua = vec![0.0; n * INCOMPATIBLE_DOFS]; + let mut kaa = Matrix::zeros(INCOMPATIBLE_DOFS); + for p in self.rule().points() { + let g = self.gradients(p.xi)?; + let w = g.det * p.weight; + let ga = self.incompatible_gradients(p.xi, g.det, &j0); + let gu = g.dx.get(..self.ty.nodes).unwrap_or_default(); + if with_uu { + accumulate_symmetric(&mut kuu.data, gu, d, w); + } + accumulate(&mut kua, INCOMPATIBLE_DOFS, gu, &ga, d, w); + accumulate_symmetric(&mut kaa.data, &ga, d, w); + } + Ok(IncompatibleParts { kuu, kua, kaa, j0 }) + } + + /// The stiffness matrix, 3n × 3n, degrees of freedom ordered node by node (x, y, z). + pub fn stiffness(&self, d: &Elasticity) -> Result { + let n = 3 * self.ty.nodes; + if self.formulation == Formulation::IncompatibleModes { + let IncompatibleParts { kuu: mut k, kua, kaa, .. } = self.incompatible_parts(d, true)?; + let chol = Cholesky::new(&kaa)?; + // K = K_uu − K_ua K_aa⁻¹ K_au, column by column of K_au. + for i in 0..n { + let mut col: Vec = kua.get(i * INCOMPATIBLE_DOFS..(i + 1) * INCOMPATIBLE_DOFS).unwrap_or_default().to_vec(); + chol.solve(&mut col); + for j in 0..n { + let row = kua.get(j * INCOMPATIBLE_DOFS..(j + 1) * INCOMPATIBLE_DOFS).unwrap_or_default(); + let v: f64 = row.iter().zip(&col).map(|(a, b)| a * b).sum(); + k.add(j, i, -v); + } + } + symmetrize(&mut k); + return Ok(k); + } + let mut k = Matrix::zeros(n); + for p in self.rule().points() { + let g = self.gradients(p.xi)?; + let gu = g.dx.get(..self.ty.nodes).unwrap_or_default(); + accumulate_symmetric(&mut k.data, gu, d, g.det * p.weight); + } + Ok(k) + } + + /// The consistent mass matrix per displacement component, n × n: m_ij = ∫ ρ Nᵢ Nⱼ dV (the + /// 3n × 3n matrix has it on each component's block). + pub fn mass(&self, density: f64) -> Result { + let n = self.ty.nodes; + let mut m = Matrix::zeros(n); + for p in self.ty.mass.points() { + let g = self.gradients(p.xi)?; + let w = density * g.det * p.weight; + let v = g.shape.values(); + for (i, a) in v.iter().enumerate() { + for (j, b) in v.iter().enumerate() { + m.add(i, j, w * a * b); + } + } + } + Ok(m) + } + + /// The lumped (diagonal) mass per node: the consistent diagonal scaled to the element's mass. + pub fn lumped_mass(&self, density: f64) -> Result> { + let m = self.mass(density)?; + // The consistent entries sum to the element's mass (the shape functions sum to one). + let total: f64 = m.data.iter().sum(); + let diag: Vec = (0..m.n).map(|i| m.get(i, i)).collect(); + let sum: f64 = diag.iter().sum(); + if !(sum > 0.0) { + return Err(self.degenerate()); + } + Ok(diag.iter().map(|d| d * total / sum).collect()) + } + + /// Consistent nodal forces of a uniform force per volume: fᵢ = ∫ Nᵢ b dV. + pub fn body_force(&self, b: [f64; 3]) -> Result> { + let mut f = vec![[0.0; 3]; self.ty.nodes]; + for p in self.ty.mass.points() { + let g = self.gradients(p.xi)?; + for (fi, n) in f.iter_mut().zip(g.shape.values()) { + for k in 0..3 { + fi[k] += n * b[k] * g.det * p.weight; + } + } + } + Ok(f) + } + + /// Strain and stress at every integration point for the nodal displacements `u`. + pub fn stresses(&self, d: &Elasticity, u: &[[f64; 3]]) -> Result> { + if u.len() != self.ty.nodes { + return Err(Error::Invalid(format!("{} displacements for a {}-node element", u.len(), self.ty.nodes))); + } + // Incompatible modes: recover the internal amplitudes α = −K_aa⁻¹ K_au u. + let alpha = if self.formulation == Formulation::IncompatibleModes { + let IncompatibleParts { kua, kaa, j0, .. } = self.incompatible_parts(d, false)?; + let flat: Vec = u.iter().flatten().copied().collect(); + let mut rhs = [0.0; INCOMPATIBLE_DOFS]; + for (i, ui) in flat.iter().enumerate() { + for (a, r) in rhs.iter_mut().enumerate() { + *r -= kua.get(i * INCOMPATIBLE_DOFS + a).copied().unwrap_or(0.0) * ui; + } + } + Cholesky::new(&kaa)?.solve(&mut rhs); + Some((rhs, j0)) + } else { + None + }; + let mut out = Vec::new(); + for p in self.rule().points() { + let g = self.gradients(p.xi)?; + let mut grad = [[0.0; 3]; 3]; // ∂u_c/∂x_b + for (dn, ui) in g.dx.iter().zip(u) { + for c in 0..3 { + for b in 0..3 { + grad[c][b] += ui[c] * dn[b]; + } + } + } + if let Some((a, j0)) = &alpha { + let ga = self.incompatible_gradients(p.xi, g.det, j0); + for (k, gk) in ga.iter().enumerate() { + for c in 0..3 { + for b in 0..3 { + grad[c][b] += a[3 * k + c] * gk[b]; + } + } + } + } + let strain = [grad[0][0], grad[1][1], grad[2][2], grad[0][1] + grad[1][0], grad[0][2] + grad[2][0], grad[1][2] + grad[2][1]]; + let mut xyz = [0.0; 3]; + for (n, x) in g.shape.values().iter().zip(self.xyz) { + for k in 0..3 { + xyz[k] += n * x[k]; + } + } + out.push(PointStress { xyz, strain, stress: stress(d, &strain) }); + } + Ok(out) + } + + /// Extrapolate values at the integration points to the element's nodes. + pub fn extrapolate(&self, values: &[[f64; 6]]) -> Result> { + let e = extrapolation(self.ty, self.rule())?; + if values.len() != e.points { + return Err(Error::Invalid(format!("{} values for {} integration points", values.len(), e.points))); + } + Ok(e.apply(values)) + } +} + +/// Make a matrix exactly symmetric (the mean of it and its transpose), so assembly can rely on it. +fn symmetrize(k: &mut Matrix) { + let n = k.n; + for i in 0..n { + for j in i + 1..n { + let v = 0.5 * (k.get(i, j) + k.get(j, i)); + if let Some(x) = k.data.get_mut(i * n + j) { + *x = v; + } + if let Some(x) = k.data.get_mut(j * n + i) { + *x = v; + } + } + } +} + +/// The strain-displacement block of one gradient g: ε = B u for a single node's (x, y, z). +fn b_block(g: &[f64; 3]) -> [[f64; 3]; 6] { + [[g[0], 0.0, 0.0], [0.0, g[1], 0.0], [0.0, 0.0, g[2]], [g[1], g[0], 0.0], [g[2], 0.0, g[0]], [0.0, g[2], g[1]]] +} + +/// out[3i + r][3j + c] += w (Bᵢᵀ D Bⱼ)[r][c] for gradients `ga` (rows) and `gb` (columns), with +/// `cols` columns in the row-major `out`. +fn accumulate(out: &mut [f64], cols: usize, ga: &[[f64; 3]], gb: &[[f64; 3]], d: &Elasticity, w: f64) { + for (j, gj) in gb.iter().enumerate() { + let bj = b_block(gj); + // D Bⱼ: 6 × 3. + let mut db = [[0.0; 3]; 6]; + for s in 0..6 { + for c in 0..3 { + db[s][c] = (0..6).map(|t| d[s][t] * bj[t][c]).sum(); + } + } + for (i, gi) in ga.iter().enumerate() { + let bi = b_block(gi); + for r in 0..3 { + for c in 0..3 { + let v: f64 = (0..6).map(|s| bi[s][r] * db[s][c]).sum(); + if let Some(x) = out.get_mut((3 * i + r) * cols + 3 * j + c) { + *x += w * v; + } + } + } + } + } +} + +/// The symmetric case of [`accumulate`] (rows and columns from the same gradients, `out` square): +/// computes the blocks with j ≥ i and mirrors them, so the result is exactly symmetric. +fn accumulate_symmetric(out: &mut [f64], g: &[[f64; 3]], d: &Elasticity, w: f64) { + let cols = 3 * g.len(); + for (j, gj) in g.iter().enumerate() { + let bj = b_block(gj); + let mut db = [[0.0; 3]; 6]; + for s in 0..6 { + for c in 0..3 { + db[s][c] = (0..6).map(|t| d[s][t] * bj[t][c]).sum(); + } + } + for (i, gi) in g.iter().enumerate().take(j + 1) { + let bi = b_block(gi); + for r in 0..3 { + // On a diagonal block only the upper triangle is computed, then mirrored. + for c in if i == j { r } else { 0 }..3 { + let v = w * (0..6).map(|s| bi[s][r] * db[s][c]).sum::(); + let (row, col) = (3 * i + r, 3 * j + c); + if let Some(x) = out.get_mut(row * cols + col) { + *x += v; + } + if row != col + && let Some(x) = out.get_mut(col * cols + row) + { + *x += v; + } + } + } + } + } +} + +/// The map from integration-point values to nodal values of an element type and rule. +#[derive(Debug, Clone)] +pub struct Extrapolation { + pub points: usize, + /// nodes × points, row-major. + pub weights: Vec, +} + +impl Extrapolation { + pub fn apply(&self, values: &[[f64; 6]]) -> Vec<[f64; 6]> { + self.weights + .chunks_exact(self.points.max(1)) + .map(|row| { + let mut v = [0.0; 6]; + for (w, x) in row.iter().zip(values) { + for k in 0..6 { + v[k] += w * x[k]; + } + } + v + }) + .collect() + } +} + +/// The least-squares extrapolation of a rule's points to the nodes with the linear corner +/// functions (Hinton & Campbell 1974); the mean when the rule has fewer points than corners. +pub fn extrapolation(ty: &ElementType, rule: Rule) -> Result { + let pts = rule.points(); + let g = pts.len(); + let c = ty.corners; + let corner_basis = |xi: [f64; 3]| match ty.family { + Family::Tetrahedron => shape::tet4(xi), + Family::Hexahedron => shape::hex8(xi), + }; + let mut corner_rows = vec![0.0; c * g]; + if g < c { + corner_rows.iter_mut().for_each(|w| *w = 1.0 / g.max(1) as f64); + } else { + // (AᵀA)⁻¹ Aᵀ with A[p][j] = corner function j at point p. + let a: Vec> = pts.iter().map(|p| corner_basis(p.xi).values().to_vec()).collect(); + let mut ata = Matrix::zeros(c); + for row in &a { + for i in 0..c { + for j in 0..c { + ata.add(i, j, row.get(i).copied().unwrap_or(0.0) * row.get(j).copied().unwrap_or(0.0)); + } + } + } + let chol = Cholesky::new(&ata)?; + for (p, row) in a.iter().enumerate() { + let mut col = row.clone(); + chol.solve(&mut col); + for (i, v) in col.iter().enumerate() { + if let Some(x) = corner_rows.get_mut(i * g + p) { + *x = *v; + } + } + } + } + let mut weights = vec![0.0; ty.nodes * g]; + for (w, v) in weights.iter_mut().zip(&corner_rows) { + *w = *v; + } + for edge in ty.kind.edges() { + if let [a, b, m] = **edge { + for p in 0..g { + let v = 0.5 * (corner_rows.get(a * g + p).copied().unwrap_or(0.0) + corner_rows.get(b * g + p).copied().unwrap_or(0.0)); + if let Some(x) = weights.get_mut(m * g + p) { + *x = v; + } + } + } + } + Ok(Extrapolation { points: g, weights }) +} diff --git a/crates/fem/src/face.rs b/crates/fem/src/face.rs new file mode 100644 index 0000000..a8f1b5b --- /dev/null +++ b/crates/fem/src/face.rs @@ -0,0 +1,76 @@ +//! Element faces: consistent nodal forces of surface loads. +//! +//! A pressure p on a face gives node i the force fᵢ = −p ∫ Nᵢ n dA, where n dA = +//! (∂x/∂ξ × ∂x/∂η) dξ dη on the face's own, possibly curved, geometry; faces are numbered so this +//! normal points out of the element ([`modelcraft_model::ElementKind::faces`]), and a positive +//! pressure pushes in. This is the consistent load vector of the virtual work of the surface load +//! (Zienkiewicz, Taylor & Zhu 2013, §6.6), integrated with the model's face rules (exact for a +//! uniform pressure on a flat face). + +use crate::{Error, Result}; + +/// ∫ Nᵢ n dA for every node of a face (node coordinates in the face's node order): the model's +/// face integration ([`modelcraft_model::nodal::face_normal_weights`]), so pressures and the +/// tractions the model distributes use the same rule. +pub fn normal_weights(xyz: &[[f64; 3]]) -> Result> { + modelcraft_model::nodal::face_normal_weights(xyz).map_err(|e| Error::Invalid(e.to_string())) +} + +/// Consistent nodal forces of a uniform pressure (positive pushes into the element). +pub fn pressure(xyz: &[[f64; 3]], magnitude: f64) -> Result> { + Ok(normal_weights(xyz)?.into_iter().map(|w| w.map(|c| -magnitude * c)).collect()) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn pressure_on_flat_faces_is_the_textbook_split() { + // A 2 × 3 quadrilateral in z = 0, normal +z: quad4 shares p A equally, quad8 gives the + // corners −A/12 and the mid-sides A/3. + let quad = [[0., 0., 0.], [2., 0., 0.], [2., 3., 0.], [0., 3., 0.]]; + let f = pressure(&quad, 1.0).unwrap(); + assert!(f.iter().all(|v| (v[2] + 1.5).abs() < 1e-14 && v[0] == 0.0 && v[1] == 0.0), "{f:?}"); + let mut q8 = quad.to_vec(); + for (a, b) in [(0, 1), (1, 2), (2, 3), (3, 0)] { + q8.push([0, 1, 2].map(|k| (quad[a][k] + quad[b][k]) / 2.0)); + } + let f = pressure(&q8, 1.0).unwrap(); + for (i, v) in f.iter().enumerate() { + let want = if i < 4 { 0.5 } else { -2.0 }; + assert!((v[2] - want).abs() < 1e-13, "{i}: {v:?}"); + } + // Triangles: tri3 A/3 each; tri6 corners 0, mid-sides A/3. + let tri = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]]; + let f = pressure(&tri, 6.0).unwrap(); + assert!(f.iter().all(|v| (v[2] + 1.0).abs() < 1e-14), "{f:?}"); + let mut t6 = tri.to_vec(); + for (a, b) in [(0, 1), (1, 2), (2, 0)] { + t6.push([0, 1, 2].map(|k| (tri[a][k] + tri[b][k]) / 2.0)); + } + let f = pressure(&t6, 6.0).unwrap(); + for (i, v) in f.iter().enumerate() { + let want = if i < 3 { 0.0 } else { -1.0 }; + assert!((v[2] - want).abs() < 1e-13, "{i}: {v:?}"); + } + assert!(pressure(&tri[..2], 1.0).is_err()); + } + + #[test] + fn pressure_on_a_curved_face_sums_to_the_projected_area() { + // A quadratic face on a quarter cylinder r = 2, z in 0..1, normal outward: the total + // pressure force is −p times the projected areas (2 in x and 2 in y). + let r: f64 = 2.0; + let at = |t: f64, z: f64| [r * t.cos(), r * t.sin(), z]; + let h = std::f64::consts::FRAC_PI_2; + // Corners counter-clockwise seen from outside: (0,0), (0,1)... with the normal outward. + let xyz = vec![at(0.0, 0.0), at(h, 0.0), at(h, 1.0), at(0.0, 1.0), at(h / 2.0, 0.0), at(h, 0.5), at(h / 2.0, 1.0), at(0.0, 0.5)]; + let w = normal_weights(&xyz).unwrap(); + let total = w.iter().fold([0.0; 3], |a, v| [a[0] + v[0], a[1] + v[1], a[2] + v[2]]); + // A quadratic arc interpolates the circle closely; the projected area is exact for the + // interpolated surface (its end points are exact). + assert!((total[0].abs() - 2.0).abs() < 1e-12 && (total[1].abs() - 2.0).abs() < 1e-12 && total[2].abs() < 1e-12, "{total:?}"); + assert!(total[0] > 0.0 && total[1] > 0.0, "this ordering points outwards: {total:?}"); + } +} diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs new file mode 100644 index 0000000..fc33091 --- /dev/null +++ b/crates/fem/src/lib.rs @@ -0,0 +1,48 @@ +//! Finite element kernels shared by ModelCraft's structural solvers. +//! +//! Pure functions over element geometry and material constants, with no solver loop and no +//! global state, so the implicit solver (`modelcraft-solver`) and the explicit solver that follows +//! it use the same element technology: +//! - [`quadrature`]: Gauss–Legendre and tetrahedron rules; +//! - [`shape`]: shape functions of tet4, tet10, hex8, hex20 and of their faces; +//! - [`material`]: isotropic and orthotropic elasticity matrices, von Mises stress; +//! - [`element`]: the element table, stiffness (full, reduced, incompatible modes), consistent +//! and lumped mass, body forces, integration-point stresses and their extrapolation to nodes; +//! - [`face`]: consistent nodal forces of pressure on element faces; +//! - [`dense`]: the small dense linear algebra inside an element. +//! +//! Conventions: three translational degrees of freedom per node, ordered node by node; stress +//! and strain 6-vectors in the order xx, yy, zz, xy, xz, yz with engineering shear strains. +#![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] +#![forbid(unsafe_code)] +// Matrix kernels over fixed small dimensions read as the textbook formulas with index loops. +#![allow(clippy::needless_range_loop)] + +pub mod dense; +pub mod element; +pub mod face; +pub mod material; +pub mod quadrature; +pub mod shape; + +pub use element::{ELEMENT_TYPES, ElementType, Family, Formulation, PointStress, Solid, element_type}; +pub use material::{Elasticity, elasticity, von_mises}; + +#[derive(Debug, Clone, PartialEq, thiserror::Error)] +pub enum Error { + #[error("invalid input: {0}")] + Invalid(String), + #[error("not supported yet: {0}")] + NotSupported(String), + /// An element's geometry is inverted or collapsed. + #[error("{0}")] + Degenerate(String), + /// A matrix that must be positive definite is not. + #[error("{0}")] + Singular(String), +} + +pub type Result = std::result::Result; + +#[cfg(test)] +mod tests; diff --git a/crates/fem/src/material.rs b/crates/fem/src/material.rs new file mode 100644 index 0000000..9540c43 --- /dev/null +++ b/crates/fem/src/material.rs @@ -0,0 +1,113 @@ +//! Linear elastic constitutive matrices. +//! +//! Stress and strain are 6-vectors in the order xx, yy, zz, xy, xz, yz, with engineering shear +//! strains (γ = 2ε), so σ = D ε. The order puts the orthotropic shear moduli G₁₂, G₁₃, G₂₃ on +//! the diagonal in the order the model lists them. Material axes are the global axes until the +//! model has material orientations. +//! +//! Orthotropic stiffness is the inverse of the compliance written with engineering constants +//! (R. M. Jones, *Mechanics of Composite Materials*, 2nd ed., Taylor & Francis 1999, §2.4): +//! S₁₁ = 1/E₁, S₁₂ = −ν₁₂/E₁, S₁₃ = −ν₁₃/E₁, S₂₃ = −ν₂₃/E₂ (ν_ji/E_j = ν_ij/E_i), and +//! 1/G on the shear diagonal. + +use modelcraft_model::Elastic; + +use crate::dense::inverse3; +use crate::{Error, Result}; + +/// A 6 × 6 elasticity matrix, row-major. +pub type Elasticity = [[f64; 6]; 6]; + +/// The isotropic elasticity matrix. +pub fn isotropic(young: f64, poisson: f64) -> Result { + if !(young > 0.0 && young.is_finite() && poisson > -1.0 && poisson < 0.5) { + return Err(Error::Invalid(format!("isotropic constants E = {young}, ν = {poisson} are not a stable material"))); + } + let lambda = young * poisson / ((1.0 + poisson) * (1.0 - 2.0 * poisson)); + let mu = young / (2.0 * (1.0 + poisson)); + let mut d = [[0.0; 6]; 6]; + for i in 0..3 { + for j in 0..3 { + d[i][j] = lambda; + } + d[i][i] = lambda + 2.0 * mu; + d[i + 3][i + 3] = mu; + } + Ok(d) +} + +/// The orthotropic elasticity matrix from engineering constants. +pub fn orthotropic(young: [f64; 3], poisson: [f64; 3], shear: [f64; 3]) -> Result { + let [e1, e2, e3] = young; + let [n12, n13, n23] = poisson; + let s = [[1.0 / e1, -n12 / e1, -n13 / e1], [-n12 / e1, 1.0 / e2, -n23 / e2], [-n13 / e1, -n23 / e2, 1.0 / e3]]; + let scale = (1.0 / e1) * (1.0 / e2) * (1.0 / e3); + let inverse = inverse3(&s).filter(|(_, det)| *det > 1e-12 * scale); + let Some((normal, _)) = inverse.filter(|_| shear.iter().all(|g| *g > 0.0 && g.is_finite())) else { + return Err(Error::Invalid("the orthotropic constants do not describe a stable material".into())); + }; + let mut d = [[0.0; 6]; 6]; + for i in 0..3 { + d[i][..3].copy_from_slice(&normal[i]); + } + for (k, g) in shear.iter().enumerate() { + d[3 + k][3 + k] = *g; + } + Ok(d) +} + +/// The elasticity matrix of a model material's elastic behaviour. +pub fn elasticity(e: &Elastic) -> Result { + match *e { + Elastic::Isotropic { young, poisson } => isotropic(young, poisson), + Elastic::EngineeringConstants { young, poisson, shear } => orthotropic(young, poisson, shear), + } +} + +/// The von Mises equivalent stress of a stress 6-vector. +pub fn von_mises(s: &[f64; 6]) -> f64 { + let [xx, yy, zz, xy, xz, yz] = *s; + let normal = (xx - yy) * (xx - yy) + (yy - zz) * (yy - zz) + (zz - xx) * (zz - xx); + (0.5 * normal + 3.0 * (xy * xy + xz * xz + yz * yz)).sqrt() +} + +/// σ = D ε. +pub fn stress(d: &Elasticity, strain: &[f64; 6]) -> [f64; 6] { + let mut s = [0.0; 6]; + for (si, row) in s.iter_mut().zip(d) { + *si = row.iter().zip(strain).map(|(a, b)| a * b).sum(); + } + s +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn orthotropic_with_isotropic_constants_is_isotropic() { + let (e, nu) = (210e3, 0.3); + let g = e / (2.0 * (1.0 + nu)); + let a = isotropic(e, nu).unwrap(); + let b = orthotropic([e; 3], [nu; 3], [g; 3]).unwrap(); + for i in 0..6 { + for j in 0..6 { + assert!((a[i][j] - b[i][j]).abs() < 1e-9 * e, "{i} {j}: {} vs {}", a[i][j], b[i][j]); + } + } + assert!(isotropic(1.0, 0.5).is_err()); + assert!(orthotropic([1.0; 3], [0.9; 3], [1.0; 3]).is_err()); + } + + #[test] + fn orthotropic_uniaxial_response() { + // Stress along 1 alone gives strains σ/E1, −ν12 σ/E1, −ν13 σ/E1. + let (young, poisson) = ([140e3, 10e3, 12e3], [0.3, 0.25, 0.4]); + let d = orthotropic(young, poisson, [5e3, 4e3, 3.8e3]).unwrap(); + let strain = [1.0 / young[0], -poisson[0] / young[0], -poisson[1] / young[0], 0.0, 0.0, 0.0]; + let s = stress(&d, &strain); + assert!((s[0] - 1.0).abs() < 1e-12 && s[1].abs() < 1e-12 && s[2].abs() < 1e-12, "{s:?}"); + assert!((von_mises(&[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]) - 1.0).abs() < 1e-15); + assert!((von_mises(&[0.0, 0.0, 0.0, 1.0, 0.0, 0.0]) - 3f64.sqrt()).abs() < 1e-15); + } +} diff --git a/crates/fem/src/quadrature.rs b/crates/fem/src/quadrature.rs new file mode 100644 index 0000000..22ebb00 --- /dev/null +++ b/crates/fem/src/quadrature.rs @@ -0,0 +1,161 @@ +//! Numerical integration rules on the reference elements. +//! +//! - Gauss–Legendre rules on [−1, 1] with 1 to 4 points (exact for polynomials of degree +//! 2n − 1), and their tensor products on the reference hexahedron [−1, 1]³. +//! - On the reference tetrahedron {ξ, η, ζ ≥ 0, ξ + η + ζ ≤ 1} (volume 1/6): the centroid rule +//! (degree 1), the symmetric four-point rule (degree 2), and the conical product rule of Stroud +//! (A. H. Stroud, *Approximate Calculation of Multiple Integrals*, Prentice-Hall 1971, §2.7): +//! the collapsed map ξ = u, η = v(1 − u), ζ = w(1 − u)(1 − v) from the unit cube, with +//! Jacobian (1 − u)²(1 − v), integrated by Gauss–Legendre in u, v and w. With n points per +//! direction it is exact for polynomials of degree 2n − 3 on the tetrahedron, with positive +//! weights only. +//! - On the reference triangle {ξ, η ≥ 0, ξ + η ≤ 1}: the same collapsed product in two +//! dimensions (Jacobian 1 − u), exact for degree 2n − 2. +//! - On the reference square [−1, 1]²: tensor Gauss–Legendre. + +/// A quadrature point: natural coordinates and weight. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Point { + pub xi: [f64; 3], + pub weight: f64, +} + +/// Gauss–Legendre abscissae and weights on [−1, 1], `n` in 1..=4 (clamped). +pub fn gauss_legendre(n: usize) -> &'static [(f64, f64)] { + const G1: [(f64, f64); 1] = [(0.0, 2.0)]; + const G2: [(f64, f64); 2] = [(-0.577_350_269_189_625_8, 1.0), (0.577_350_269_189_625_8, 1.0)]; + const G3: [(f64, f64); 3] = [(-0.774_596_669_241_483_4, 5.0 / 9.0), (0.0, 8.0 / 9.0), (0.774_596_669_241_483_4, 5.0 / 9.0)]; + const G4: [(f64, f64); 4] = [ + (-0.861_136_311_594_052_6, 0.347_854_845_137_453_9), + (-0.339_981_043_584_856_3, 0.652_145_154_862_546_1), + (0.339_981_043_584_856_3, 0.652_145_154_862_546_1), + (0.861_136_311_594_052_6, 0.347_854_845_137_453_9), + ]; + match n { + 0 | 1 => &G1, + 2 => &G2, + 3 => &G3, + _ => &G4, + } +} + +/// An integration rule on a reference element. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub enum Rule { + /// n × n × n Gauss–Legendre on the reference hexahedron. + Hex(u8), + /// The tetrahedron's centroid, weight 1/6 (degree 1). + TetCentroid, + /// Four symmetric points on the tetrahedron, weight 1/24 each (degree 2). + Tet4, + /// Stroud's conical product on the tetrahedron, n points per direction. + TetConical(u8), + /// n × n Gauss–Legendre on the reference square (z = 0). + Quad(u8), + /// The collapsed product on the reference triangle, n points per direction (z = 0). + TriConical(u8), +} + +impl Rule { + pub fn points(self) -> Vec { + match self { + Rule::Hex(n) => { + let g = gauss_legendre(usize::from(n)); + let mut out = Vec::with_capacity(g.len().pow(3)); + // ξ fastest, then η, then ζ: the corner order of the reference hexahedron for 2 × 2 × 2. + for &(z, wz) in g { + for &(y, wy) in g { + for &(x, wx) in g { + out.push(Point { xi: [x, y, z], weight: wx * wy * wz }); + } + } + } + out + } + Rule::TetCentroid => vec![Point { xi: [0.25; 3], weight: 1.0 / 6.0 }], + Rule::Tet4 => { + // (5 ∓ √5)/20: the barycentric coordinates of the degree-2 rule. + let a = 0.585_410_196_624_968_5; + let b = 0.138_196_601_125_010_5; + // Point k sits near corner k (barycentric a there, b elsewhere). + [[b, b, b], [a, b, b], [b, a, b], [b, b, a]].into_iter().map(|xi| Point { xi, weight: 1.0 / 24.0 }).collect() + } + Rule::TetConical(n) => { + let g: Vec<(f64, f64)> = gauss_legendre(usize::from(n)).iter().map(|&(x, w)| ((x + 1.0) / 2.0, w / 2.0)).collect(); + let mut out = Vec::with_capacity(g.len().pow(3)); + for &(u, wu) in &g { + for &(v, wv) in &g { + for &(w, ww) in &g { + let xi = [u, v * (1.0 - u), w * (1.0 - u) * (1.0 - v)]; + out.push(Point { xi, weight: wu * wv * ww * (1.0 - u) * (1.0 - u) * (1.0 - v) }); + } + } + } + out + } + Rule::Quad(n) => { + let g = gauss_legendre(usize::from(n)); + g.iter().flat_map(|&(y, wy)| g.iter().map(move |&(x, wx)| Point { xi: [x, y, 0.0], weight: wx * wy })).collect() + } + Rule::TriConical(n) => { + let g: Vec<(f64, f64)> = gauss_legendre(usize::from(n)).iter().map(|&(x, w)| ((x + 1.0) / 2.0, w / 2.0)).collect(); + g.iter() + .flat_map(|&(u, wu)| g.iter().map(move |&(v, wv)| Point { xi: [u, v * (1.0 - u), 0.0], weight: wu * wv * (1.0 - u) })) + .collect() + } + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// ∫ ξ^a η^b ζ^c over the reference tetrahedron = a! b! c! / (a + b + c + 3)!. + fn tet_monomial(a: u32, b: u32, c: u32) -> f64 { + let f = |n: u32| (1..=n).map(f64::from).product::(); + f(a) * f(b) * f(c) / f(a + b + c + 3) + } + + fn integrate(rule: Rule, a: i32, b: i32, c: i32) -> f64 { + rule.points().iter().map(|p| p.weight * p.xi[0].powi(a) * p.xi[1].powi(b) * p.xi[2].powi(c)).sum() + } + + #[test] + fn tetrahedron_rules_are_exact_to_their_degree() { + for (rule, degree) in [(Rule::TetCentroid, 1), (Rule::Tet4, 2), (Rule::TetConical(3), 3), (Rule::TetConical(4), 5)] { + for a in 0..=degree { + for b in 0..=degree - a { + for c in 0..=degree - a - b { + let exact = tet_monomial(a, b, c); + let got = integrate(rule, a as i32, b as i32, c as i32); + assert!((got - exact).abs() < 1e-14, "{rule:?} ξ^{a} η^{b} ζ^{c}: {got} vs {exact}"); + } + } + } + } + } + + #[test] + fn hexahedron_and_face_rules() { + for n in 1..=4u8 { + let d = 2 * i32::from(n) - 1; + for a in 0..=d { + // ∫_{-1}^{1} x^a dx. + let exact_1d = if a % 2 == 1 { 0.0 } else { 2.0 / f64::from(a + 1) }; + let got = integrate(Rule::Hex(n), a, 0, 0); + assert!((got - exact_1d * 4.0).abs() < 1e-14, "hex {n} x^{a}"); + let got = integrate(Rule::Quad(n), 0, a, 0); + assert!((got - exact_1d * 2.0).abs() < 1e-14, "quad {n} y^{a}"); + } + } + // Triangle: ∫ ξ^a η^b = a! b! / (a + b + 2)!. + let f = |n: i32| (1..=n).map(f64::from).product::(); + for a in 0..=4 { + for b in 0..=4 - a { + let got = integrate(Rule::TriConical(3), a, b, 0); + assert!((got - f(a) * f(b) / f(a + b + 2)).abs() < 1e-14, "tri ξ^{a} η^{b}"); + } + } + } +} diff --git a/crates/fem/src/shape.rs b/crates/fem/src/shape.rs new file mode 100644 index 0000000..3036ec2 --- /dev/null +++ b/crates/fem/src/shape.rs @@ -0,0 +1,192 @@ +//! Shape functions and their derivatives in natural coordinates. +//! +//! Node order is the model's (Abaqus/VTK, [`modelcraft_model::ElementKind`]): +//! - tetrahedra on the reference element with corners (0,0,0), (1,0,0), (0,1,0), (0,0,1) and the +//! volume coordinates L₁ = 1 − ξ − η − ζ, L₂ = ξ, L₃ = η, L₄ = ζ; tet10 mid-side nodes on edges +//! 1-2, 2-3, 3-1, 1-4, 2-4, 3-4 with N = 4 Lₐ L_b, corners N = L(2L − 1); +//! - hexahedra on [−1, 1]³, corners 1–4 at ζ = −1 counter-clockwise from (−1,−1), 5–8 above +//! them; hex20 is the serendipity element (corners +//! N = ⅛(1 + ξξᵢ)(1 + ηηᵢ)(1 + ζζᵢ)(ξξᵢ + ηηᵢ + ζζᵢ − 2), a mid-side node with ξᵢ = 0 +//! N = ¼(1 − ξ²)(1 + ηηᵢ)(1 + ζζᵢ), and likewise), mid-side nodes 9–12 on the bottom edges, +//! 13–16 on the top edges, 17–20 on the vertical edges (mid-side positions come from +//! [`modelcraft_model::ElementKind::edges`], the one place the node order is written down). +//! +//! The functions are the textbook ones (O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, *The Finite +//! Element Method: Its Basis and Fundamentals*, 7th ed., Butterworth-Heinemann 2013, ch. 6; +//! T. J. R. Hughes, *The Finite Element Method*, Dover 2000, ch. 3). + +use modelcraft_model::ElementKind; + +/// Most nodes an element in this crate has (hex20). +pub const MAX_NODES: usize = 20; + +/// Shape function values and natural derivatives at one point. +#[derive(Debug, Clone, Copy)] +pub struct Shape { + pub count: usize, + pub n: [f64; MAX_NODES], + /// ∂Nᵢ/∂ξ, ∂Nᵢ/∂η, ∂Nᵢ/∂ζ. + pub dn: [[f64; 3]; MAX_NODES], +} + +impl Shape { + fn new(count: usize) -> Shape { + Shape { count, n: [0.0; MAX_NODES], dn: [[0.0; 3]; MAX_NODES] } + } + + pub fn values(&self) -> &[f64] { + self.n.get(..self.count).unwrap_or_default() + } + + pub fn derivatives(&self) -> &[[f64; 3]] { + self.dn.get(..self.count).unwrap_or_default() + } +} + +/// Natural coordinates of the hexahedron's corners. +pub const HEX_CORNERS: [[f64; 3]; 8] = + [[-1., -1., -1.], [1., -1., -1.], [1., 1., -1.], [-1., 1., -1.], [-1., -1., 1.], [1., -1., 1.], [1., 1., 1.], [-1., 1., 1.]]; + +/// Natural coordinates of the tetrahedron's corners. +pub const TET_CORNERS: [[f64; 3]; 4] = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]; + +/// Four-node tetrahedron. +pub fn tet4(xi: [f64; 3]) -> Shape { + let mut s = Shape::new(4); + s.n[0] = 1.0 - xi[0] - xi[1] - xi[2]; + s.n[1] = xi[0]; + s.n[2] = xi[1]; + s.n[3] = xi[2]; + s.dn[0] = [-1.0, -1.0, -1.0]; + s.dn[1] = [1.0, 0.0, 0.0]; + s.dn[2] = [0.0, 1.0, 0.0]; + s.dn[3] = [0.0, 0.0, 1.0]; + s +} + +/// Ten-node tetrahedron. +pub fn tet10(xi: [f64; 3]) -> Shape { + let lin = tet4(xi); + let l = [lin.n[0], lin.n[1], lin.n[2], lin.n[3]]; + let dl = [lin.dn[0], lin.dn[1], lin.dn[2], lin.dn[3]]; + let mut s = Shape::new(10); + for i in 0..4 { + s.n[i] = l[i] * (2.0 * l[i] - 1.0); + s.dn[i] = dl[i].map(|d| (4.0 * l[i] - 1.0) * d); + } + // Mid-side nodes on the model's tet10 edges. + for edge in ElementKind::Tet10.edges() { + if let [a, b, m] = **edge { + s.n[m] = 4.0 * l[a] * l[b]; + s.dn[m] = [0, 1, 2].map(|c| 4.0 * (dl[a][c] * l[b] + l[a] * dl[b][c])); + } + } + s +} + +/// Eight-node hexahedron (trilinear). +pub fn hex8(xi: [f64; 3]) -> Shape { + let mut s = Shape::new(8); + for (i, c) in HEX_CORNERS.iter().enumerate() { + let f = [0, 1, 2].map(|k| 1.0 + xi[k] * c[k]); + s.n[i] = f[0] * f[1] * f[2] / 8.0; + s.dn[i] = [c[0] * f[1] * f[2] / 8.0, c[1] * f[0] * f[2] / 8.0, c[2] * f[0] * f[1] / 8.0]; + } + s +} + +/// Twenty-node serendipity hexahedron. +pub fn hex20(xi: [f64; 3]) -> Shape { + let mut s = Shape::new(20); + let [x, y, z] = xi; + for (i, c) in HEX_CORNERS.iter().enumerate() { + let f = [0, 1, 2].map(|k| 1.0 + xi[k] * c[k]); + let g = x * c[0] + y * c[1] + z * c[2] - 2.0; + s.n[i] = f[0] * f[1] * f[2] * g / 8.0; + // ∂/∂ξ_k: c_k f_l f_m g / 8 + f_0 f_1 f_2 c_k / 8. + s.dn[i] = [c[0] * f[1] * f[2] * (g + f[0]) / 8.0, c[1] * f[0] * f[2] * (g + f[1]) / 8.0, c[2] * f[0] * f[1] * (g + f[2]) / 8.0]; + } + // Mid-side nodes on the model's hex20 edges, at the midpoints of their corners. + for edge in ElementKind::Hex20.edges() { + let [a, b, i] = **edge else { continue }; + let c = [0, 1, 2].map(|k| (HEX_CORNERS[a][k] + HEX_CORNERS[b][k]) / 2.0); + // The direction along the edge has cᵢ = 0: N = ¼(1 − t²)(1 + a aᵢ)(1 + b bᵢ). + let along = if c[0] == 0.0 { + 0 + } else if c[1] == 0.0 { + 1 + } else { + 2 + }; + let (p, q) = ((along + 1) % 3, (along + 2) % 3); + let t = xi[along]; + let fp = 1.0 + xi[p] * c[p]; + let fq = 1.0 + xi[q] * c[q]; + s.n[i] = (1.0 - t * t) * fp * fq / 4.0; + let mut d = [0.0; 3]; + d[along] = -2.0 * t * fp * fq / 4.0; + d[p] = (1.0 - t * t) * c[p] * fq / 4.0; + d[q] = (1.0 - t * t) * fp * c[q] / 4.0; + s.dn[i] = d; + } + s +} + +#[cfg(test)] +mod tests { + use super::*; + + fn nodes_of(kind: &str) -> Vec<[f64; 3]> { + match kind { + "tet4" => TET_CORNERS.to_vec(), + "tet10" => { + let mut v = TET_CORNERS.to_vec(); + for e in ElementKind::Tet10.edges() { + v.push([0, 1, 2].map(|k| (TET_CORNERS[e[0]][k] + TET_CORNERS[e[1]][k]) / 2.0)); + } + v + } + "hex8" => HEX_CORNERS.to_vec(), + _ => { + let mut v = HEX_CORNERS.to_vec(); + for e in ElementKind::Hex20.edges() { + v.push([0, 1, 2].map(|k| (HEX_CORNERS[e[0]][k] + HEX_CORNERS[e[1]][k]) / 2.0)); + } + v + } + } + } + + #[test] + fn kronecker_delta_partition_of_unity_and_derivatives() { + let kinds: [(&str, fn([f64; 3]) -> Shape); 4] = [("tet4", tet4), ("tet10", tet10), ("hex8", hex8), ("hex20", hex20)]; + for (name, f) in kinds { + let nodes = nodes_of(name); + for (i, p) in nodes.iter().enumerate() { + let s = f(*p); + for (j, v) in s.values().iter().enumerate() { + assert!((v - if i == j { 1.0 } else { 0.0 }).abs() < 1e-14, "{name} N{j} at node {i}: {v}"); + } + } + // At an arbitrary point: ΣN = 1, ΣdN = 0, Σ N x = x (linear completeness), and the + // derivatives match central differences. + let p = [0.21, 0.17, 0.33]; + let s = f(p); + assert!((s.values().iter().sum::() - 1.0).abs() < 1e-14); + for k in 0..3 { + assert!(s.derivatives().iter().map(|d| d[k]).sum::().abs() < 1e-13); + let x: f64 = s.values().iter().zip(&nodes).map(|(n, x)| n * x[k]).sum(); + assert!((x - p[k]).abs() < 1e-14, "{name}"); + let h = 1e-6; + let (mut a, mut b) = (p, p); + a[k] += h; + b[k] -= h; + let (sa, sb) = (f(a), f(b)); + for j in 0..s.count { + let fd = (sa.n[j] - sb.n[j]) / (2.0 * h); + assert!((fd - s.dn[j][k]).abs() < 1e-8, "{name} dN{j}/dξ{k}: {fd} vs {}", s.dn[j][k]); + } + } + } + } +} diff --git a/crates/fem/src/tests.rs b/crates/fem/src/tests.rs new file mode 100644 index 0000000..ddc0cd8 --- /dev/null +++ b/crates/fem/src/tests.rs @@ -0,0 +1,206 @@ +//! Element-level checks: rigid-body modes, constant-strain reproduction (the single-element form +//! of the patch test), mass and load totals, extrapolation. + +use modelcraft_model::ElementKind; + +use crate::dense::{Cholesky, Matrix}; +use crate::element::{Formulation, Solid}; +use crate::material::{isotropic, orthotropic, stress}; + +/// A distorted but valid element of each kind, mid-side nodes off the straight edges for the +/// quadratic kinds when `curved`. +fn distorted(kind: ElementKind, curved: bool) -> Vec<[f64; 3]> { + let corners: Vec<[f64; 3]> = match kind { + ElementKind::Tet4 | ElementKind::Tet10 => vec![[0.1, -0.05, 0.0], [1.2, 0.1, 0.05], [0.2, 0.9, -0.1], [0.05, 0.15, 1.1]], + _ => vec![ + [0.0, 0.0, 0.0], + [2.1, 0.1, -0.1], + [2.3, 1.2, 0.05], + [-0.1, 0.9, 0.1], + [0.1, -0.1, 0.8], + [2.0, 0.05, 1.1], + [2.2, 1.1, 0.9], + [0.05, 1.05, 1.2], + ], + }; + let mut xyz = corners.clone(); + for (i, e) in kind.edges().iter().enumerate() { + if let [a, b, _] = **e { + let bump = if curved { 0.03 * ((i % 3) as f64 - 1.0) } else { 0.0 }; + xyz.push([0, 1, 2].map(|k| (corners[a][k] + corners[b][k]) / 2.0 + bump)); + } + } + xyz +} + +const CASES: [(ElementKind, Formulation); 6] = [ + (ElementKind::Tet4, Formulation::Full), + (ElementKind::Tet10, Formulation::Full), + (ElementKind::Hex8, Formulation::Full), + (ElementKind::Hex8, Formulation::IncompatibleModes), + (ElementKind::Hex20, Formulation::Full), + (ElementKind::Hex20, Formulation::Reduced), +]; + +fn rigid_modes(xyz: &[[f64; 3]]) -> Vec> { + let mut modes = Vec::new(); + for k in 0..3 { + modes.push(xyz.iter().flat_map(|_| [0, 1, 2].map(|c| if c == k { 1.0 } else { 0.0 })).collect()); + } + for axis in 0..3 { + // ω × x for ω along `axis`. + let w = [0, 1, 2].map(|c| if c == axis { 1.0 } else { 0.0 }); + modes.push(xyz.iter().flat_map(|x| [w[1] * x[2] - w[2] * x[1], w[2] * x[0] - w[0] * x[2], w[0] * x[1] - w[1] * x[0]]).collect()); + } + modes +} + +#[test] +fn stiffness_is_symmetric_with_exactly_six_rigid_body_modes() { + let d = isotropic(1000.0, 0.3).unwrap(); + for (kind, f) in CASES { + let xyz = distorted(kind, true); + let k = Solid::new(kind, &xyz, f).unwrap().stiffness(&d).unwrap(); + let scale = k.data.iter().fold(0.0f64, |m, v| m.max(v.abs())); + for i in 0..k.n { + for j in 0..k.n { + assert_eq!(k.get(i, j), k.get(j, i), "{kind:?} {f:?} symmetric"); + } + } + for m in rigid_modes(&xyz) { + let r = k.mul_vec(&m); + assert!(r.iter().all(|v| v.abs() < 1e-10 * scale), "{kind:?} {f:?}: rigid motion strains the element"); + } + // Holding six degrees of freedom (node 1 fully, node 2 in y and z, node 3 in z) leaves a + // positive definite matrix: no other zero-energy mode. Reduced hex20 has spurious modes + // on a single element, so it is skipped. + if f != Formulation::Reduced { + let mut held = k.clone(); + for dof in [0, 1, 2, 4, 5, 8] { + held.add(dof, dof, 1e3 * scale); + } + assert!(Cholesky::new(&held).is_ok(), "{kind:?} {f:?} has a spurious zero-energy mode"); + } + } +} + +#[test] +fn constant_strain_is_reproduced_exactly() { + // u = A x for a general A: every integration point sees ε(A), the stiffness's internal + // forces balance (K u is the boundary load of a uniform stress), and the incompatible modes + // stay idle. + let a = [[1e-3, 2e-4, -3e-4], [5e-4, -2e-3, 1e-4], [-1e-4, 3e-4, 7e-4]]; + let strain = [a[0][0], a[1][1], a[2][2], a[0][1] + a[1][0], a[0][2] + a[2][0], a[1][2] + a[2][1]]; + for d in [isotropic(210e3, 0.3).unwrap(), orthotropic([140e3, 10e3, 12e3], [0.3, 0.25, 0.4], [5e3, 4e3, 3.8e3]).unwrap()] { + let want = stress(&d, &strain); + for (kind, f) in CASES { + // Curved mid-sides keep the field in the isoparametric span only for straight edges. + let xyz = distorted(kind, false); + let u: Vec<[f64; 3]> = xyz.iter().map(|x| [0, 1, 2].map(|i| (0..3).map(|j| a[i][j] * x[j]).sum())).collect(); + let e = Solid::new(kind, &xyz, f).unwrap(); + for p in e.stresses(&d, &u).unwrap() { + for k in 0..6 { + assert!((p.strain[k] - strain[k]).abs() < 1e-14, "{kind:?} {f:?} strain {k}"); + assert!((p.stress[k] - want[k]).abs() < 1e-9 * want[0].abs().max(1.0), "{kind:?} {f:?} stress {k}"); + } + } + let nodal = e.extrapolate(&e.stresses(&d, &u).unwrap().iter().map(|p| p.stress).collect::>()).unwrap(); + assert_eq!(nodal.len(), xyz.len()); + assert!(nodal.iter().all(|s| (0..6).all(|k| (s[k] - want[k]).abs() < 1e-9 * want[0].abs())), "{kind:?} {f:?}"); + } + } +} + +#[test] +fn mass_volume_and_body_force_totals() { + let rho = 7.85e-9; + for (kind, f) in CASES { + let xyz = distorted(kind, true); + let e = Solid::new(kind, &xyz, f).unwrap(); + let v = e.volume().unwrap(); + assert!(v > 0.0); + let m = e.mass(rho).unwrap(); + let total: f64 = m.data.iter().sum(); + assert!((total - rho * v).abs() < 1e-12 * rho * v, "{kind:?}: {total} vs {}", rho * v); + let lumped = e.lumped_mass(rho).unwrap(); + assert!(lumped.iter().all(|x| *x > 0.0), "{kind:?}: lumped masses are positive"); + assert!((lumped.iter().sum::() - rho * v).abs() < 1e-12 * rho * v); + let b = e.body_force([0.0, 0.0, -9.81]).unwrap(); + assert!((b.iter().map(|f| f[2]).sum::() + 9.81 * v).abs() < 1e-12 * v); + } + // A tetrahedron's volume is a sixth of its edge triple product. + let xyz = distorted(ElementKind::Tet4, false); + let e = |i: usize| [0, 1, 2].map(|k| xyz[i][k] - xyz[0][k]); + let (a, b, c) = (e(1), e(2), e(3)); + let exact = (a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0])) / 6.0; + let v = Solid::new(ElementKind::Tet4, &xyz, Formulation::Full).unwrap().volume().unwrap(); + assert!((v - exact).abs() < 1e-15); +} + +#[test] +fn extrapolation_reproduces_linear_fields_on_affine_elements() { + // Values linear in x at the integration points come back exactly at the nodes (tet4 has + // one point and gives the mean). + let f = |x: [f64; 3]| 1.0 + 2.0 * x[0] - 0.5 * x[1] + 3.0 * x[2]; + let affine = |kind: ElementKind| -> Vec<[f64; 3]> { + let mut xyz: Vec<[f64; 3]> = match kind { + ElementKind::Tet10 => vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]], + _ => vec![[0., 0., 0.], [2., 0., 0.], [2.5, 1., 0.], [0.5, 1., 0.], [0., 0., 1.], [2., 0., 1.], [2.5, 1., 1.], [0.5, 1., 1.]], + }; + let c = xyz.clone(); + for e in kind.edges() { + if let [a, b, _] = **e { + xyz.push([0, 1, 2].map(|k| (c[a][k] + c[b][k]) / 2.0)); + } + } + xyz + }; + for (kind, form) in CASES.into_iter().filter(|(k, _)| *k != ElementKind::Tet4) { + let xyz = affine(kind); + let e = Solid::new(kind, &xyz, form).unwrap(); + let values: Vec<[f64; 6]> = e.integration_points().iter().map(|p| [f(e.position(p.xi)); 6]).collect(); + let nodal = e.extrapolate(&values).unwrap(); + for (x, v) in xyz.iter().zip(&nodal) { + assert!((v[0] - f(*x)).abs() < 1e-12, "{kind:?} {form:?} at {x:?}: {} vs {}", v[0], f(*x)); + } + } +} + +#[test] +fn bad_input_is_an_error() { + let xyz = distorted(ElementKind::Hex8, false); + assert!(Solid::new(ElementKind::Hex8, &xyz[..7], Formulation::Full).is_err()); + assert!(Solid::new(ElementKind::Hex8, &xyz, Formulation::Reduced).is_err()); + assert!(Solid::new(ElementKind::Tet4, &xyz[..4], Formulation::IncompatibleModes).is_err()); + assert!(Solid::new(ElementKind::Wedge6, &xyz[..6], Formulation::Full).is_err()); + let mut nan = xyz.clone(); + nan[3][1] = f64::NAN; + assert!(Solid::new(ElementKind::Hex8, &nan, Formulation::Full).is_err()); + // An inverted element (top and bottom swapped) and a flat one. + let mut inverted = xyz.clone(); + inverted.rotate_left(4); + let d = isotropic(1.0, 0.3).unwrap(); + assert!(Solid::new(ElementKind::Hex8, &inverted, Formulation::Full).unwrap().stiffness(&d).is_err()); + let flat: Vec<[f64; 3]> = xyz.iter().map(|p| [p[0], p[1], 0.0]).collect(); + assert!(Solid::new(ElementKind::Hex8, &flat, Formulation::IncompatibleModes).unwrap().stiffness(&d).is_err()); + let e = Solid::new(ElementKind::Hex8, &xyz, Formulation::Full).unwrap(); + assert!(e.stresses(&d, &[[0.0; 3]; 3]).is_err()); + assert!(e.extrapolate(&[[0.0; 6]; 3]).is_err()); + let _ = Matrix::zeros(0).mul_vec(&[]); +} + +#[test] +fn section_formulations() { + use modelcraft_model::{Formulation as F, Integration}; + let full = F::default(); + let reduced = F { integration: Integration::Reduced, ..F::default() }; + let im = F { incompatible_modes: true, ..F::default() }; + let hybrid = F { hybrid: true, ..F::default() }; + assert_eq!(Formulation::from_section(ElementKind::Hex8, &full), Ok(Formulation::Full)); + assert_eq!(Formulation::from_section(ElementKind::Hex8, &im), Ok(Formulation::IncompatibleModes)); + assert_eq!(Formulation::from_section(ElementKind::Hex20, &reduced), Ok(Formulation::Reduced)); + assert!(Formulation::from_section(ElementKind::Hex8, &reduced).unwrap_err().to_string().contains("not supported yet")); + assert!(Formulation::from_section(ElementKind::Tet10, &im).is_err()); + assert!(Formulation::from_section(ElementKind::Tet4, &hybrid).is_err()); + assert!(Formulation::from_section(ElementKind::Quad4, &full).is_err()); +} diff --git a/crates/model/src/nodal.rs b/crates/model/src/nodal.rs index 613c6e8..1a22dc0 100644 --- a/crates/model/src/nodal.rs +++ b/crates/model/src/nodal.rs @@ -95,6 +95,48 @@ fn rule(corners: usize) -> Vec<(f64, f64, f64)> { } } +/// Visit the quadrature points of a face (node coordinates in the face's node order): the shape +/// function values and the weighted normal (∂x/∂ξ × ∂x/∂η) w, whose length is the point's share of +/// the area and whose direction is the face's right-hand normal. +fn integrate_face(xyz: &[[f64; 3]], mut visit: impl FnMut(&[f64], [f64; 3])) -> Result<()> { + if xyz.iter().flatten().any(|c| !c.is_finite()) { + return Err(Error::Invalid("face coordinates must be finite".into())); + } + let (corners, shape): (usize, fn(usize, f64, f64) -> Shape) = match xyz.len() { + 3 | 6 => (3, triangle), + 4 | 8 => (4, quadrilateral), + n => return Err(Error::Invalid(format!("a face has 3, 4, 6 or 8 nodes, not {n}"))), + }; + for (xi, eta, w) in rule(corners) { + let s = shape(xyz.len(), xi, eta); + let mut tx = [0.0; 3]; + let mut te = [0.0; 3]; + for ((p, dx), de) in xyz.iter().zip(&s.dxi).zip(&s.deta) { + for k in 0..3 { + tx[k] += dx * p[k]; + te[k] += de * p[k]; + } + } + let c = [tx[1] * te[2] - tx[2] * te[1], tx[2] * te[0] - tx[0] * te[2], tx[0] * te[1] - tx[1] * te[0]]; + visit(&s.n, c.map(|v| v * w)); + } + Ok(()) +} + +/// ∫ Nᵢ n dA for every node of a face: the consistent nodal forces of a unit pressure pulling +/// along the face's right-hand (outward, for solids) normal. +pub fn face_normal_weights(xyz: &[[f64; 3]]) -> Result> { + let mut out = vec![[0.0; 3]; xyz.len()]; + integrate_face(xyz, |n, nda| { + for (o, ni) in out.iter_mut().zip(n) { + for k in 0..3 { + o[k] += ni * nda[k]; + } + } + })?; + Ok(out) +} + impl Model { /// The nodal weights `∫ N_i dA` of a set of element faces, summed per node. Their sum is the /// faces' area. @@ -104,27 +146,16 @@ impl Model { let labels = self.face_nodes(f).ok_or_else(|| Error::NotFound(format!("element {} face {}", f.element, f.face)))?; let xyz: Vec<[f64; 3]> = labels.iter().map(|n| self.node(*n).map(|p| p.xyz).ok_or_else(|| Error::NotFound(format!("node {n}")))).collect::>()?; - let (corners, shape): (usize, fn(usize, f64, f64) -> Shape) = match labels.len() { - 3 | 6 => (3, triangle), - 4 | 8 => (4, quadrilateral), - n => return Err(Error::Invalid(format!("element {} face {} has {n} nodes, not a surface face", f.element, f.face))), - }; - for (xi, eta, w) in rule(corners) { - let s = shape(labels.len(), xi, eta); - let mut tx = [0.0; 3]; - let mut te = [0.0; 3]; - for ((p, dx), de) in xyz.iter().zip(&s.dxi).zip(&s.deta) { - for k in 0..3 { - tx[k] += dx * p[k]; - te[k] += de * p[k]; - } + integrate_face(&xyz, |n, nda| { + let da = (nda[0] * nda[0] + nda[1] * nda[1] + nda[2] * nda[2]).sqrt(); + for (label, ni) in labels.iter().zip(n) { + *out.entry(*label).or_default() += ni * da; } - let c = [tx[1] * te[2] - tx[2] * te[1], tx[2] * te[0] - tx[0] * te[2], tx[0] * te[1] - tx[1] * te[0]]; - let da = (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]).sqrt() * w; - for (label, n) in labels.iter().zip(&s.n) { - *out.entry(*label).or_default() += n * da; - } - } + }) + .map_err(|e| match e { + Error::Invalid(m) => Error::Invalid(format!("element {} face {}: {m}", f.element, f.face)), + e => e, + })?; } Ok(out) } diff --git a/crates/solver/Cargo.toml b/crates/solver/Cargo.toml new file mode 100644 index 0000000..84c99b2 --- /dev/null +++ b/crates/solver/Cargo.toml @@ -0,0 +1,15 @@ +[package] +name = "modelcraft-solver" +description = "ModelCraft's implicit solver: assembly, boundary conditions, sparse Cholesky, linear static and frequency analysis" +version.workspace = true +edition.workspace = true +license.workspace = true + +[lints] +workspace = true + +[dependencies] +modelcraft-model = { workspace = true } +modelcraft-fem = { workspace = true } +thiserror = { workspace = true } +serde = { workspace = true } diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs new file mode 100644 index 0000000..66e7b30 --- /dev/null +++ b/crates/solver/src/analysis.rs @@ -0,0 +1,431 @@ +//! Linear static and frequency analysis of a model's step, and their results. + +use modelcraft_fem::element::{Extrapolation, extrapolation}; +use modelcraft_fem::quadrature::Rule; +use modelcraft_fem::{PointStress, von_mises}; +use modelcraft_model::{ElementKind, Model, Procedure, Resolved, Step}; +use serde::Serialize; + +use crate::cholesky::{Factor, FactorError, UpperColumns}; +use crate::eigen::{Operator, Wanted, lanczos}; +use crate::ordering::nested_dissection; +use crate::sparse::{Csr, NodeGraph, add_block, add_diagonal, add_scalar}; +use crate::system::{Mesh, equations, loads, prescribed}; +use crate::{Error, Limits, Result}; + +/// Size and effort of a solution. +#[derive(Debug, Clone, Default, PartialEq, Serialize)] +pub struct Stats { + /// Degrees of freedom of the mesh (3 per node). + pub dofs: usize, + /// Free degrees of freedom solved for. + pub equations: usize, + /// Entries in the Cholesky factor. + pub factor_entries: usize, + /// Lanczos vectors used (frequency analysis). + pub lanczos_vectors: usize, +} + +/// Stresses of one element. +#[derive(Debug, Clone, PartialEq)] +pub struct ElementStress { + pub label: u32, + pub kind: ElementKind, + /// At the integration points. + pub points: Vec, + /// Extrapolated to the element's nodes, in its node order. + pub nodal: Vec<[f64; 6]>, +} + +/// The result of a static step. +#[derive(Debug, Clone, PartialEq)] +pub struct StaticResult { + pub step: String, + /// Node labels with degrees of freedom, sorted; the per-node vectors follow this order. + pub nodes: Vec, + pub displacement: Vec<[f64; 3]>, + /// Reaction forces (zero where nothing is prescribed). + pub reaction: Vec<[f64; 3]>, + /// Stress averaged at the nodes over the elements that share them (xx, yy, zz, xy, xz, yz). + pub stress: Vec<[f64; 6]>, + /// Von Mises stress of the averaged nodal stress. + pub von_mises: Vec, + pub elements: Vec, + /// Sum of the applied loads and of the reactions. + pub applied_total: [f64; 3], + pub reaction_total: [f64; 3], + pub stats: Stats, + pub warnings: Vec, +} + +/// One natural mode. +#[derive(Debug, Clone, PartialEq)] +pub struct Mode { + /// λ = ω², in (rad per time unit)². + pub eigenvalue: f64, + /// Cycles per time unit (0 for a rigid-body mode). + pub frequency: f64, + /// Mass-normalized (φᵀ M φ = 1), per node in [`FrequencyResult::nodes`] order. + pub shape: Vec<[f64; 3]>, + /// Participation factors Γ = φᵀ M r for unit translations along X, Y, Z. + pub participation: [f64; 3], + /// Effective mass Γ² along X, Y, Z. + pub effective_mass: [f64; 3], +} + +/// The result of a frequency step. +#[derive(Debug, Clone, PartialEq)] +pub struct FrequencyResult { + pub step: String, + pub nodes: Vec, + pub modes: Vec, + /// Mass of the model. + pub total_mass: f64, + pub stats: Stats, + pub warnings: Vec, +} + +/// Mass matrix choice for a frequency analysis. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, serde::Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum MassKind { + #[default] + Consistent, + Lumped, +} + +/// Options of a frequency analysis. +#[derive(Debug, Clone, Copy, Default)] +pub struct FrequencyOptions { + pub mass: MassKind, + /// The spectral shift σ (eigenvalue units); by default 0 for a supported structure and a + /// small negative value for a free one. + pub shift: Option, +} + +/// The step named `name`, or the only step of the procedure when no name is given; it must be +/// of that procedure. +fn select_step<'m>(model: &'m Model, name: Option<&str>, wanted: fn(&Procedure) -> bool, what: &str) -> Result<&'m Step> { + let step = match name { + Some(n) => model.step(n).ok_or_else(|| Error::Invalid(format!("step `{n}` does not exist")))?, + None => { + let matching: Vec<&Step> = model.steps().iter().filter(|s| wanted(&s.procedure)).collect(); + match matching.as_slice() { + [one] => *one, + [] => return Err(Error::Invalid(format!("the model has no {what} step"))), + many => { + let names: Vec<&str> = many.iter().map(|s| s.name.as_str()).collect(); + return Err(Error::Invalid(format!("name the step: the model has {what} steps {}", names.join(", ")))); + } + } + } + }; + if !wanted(&step.procedure) { + return Err(Error::Invalid(format!("step `{}` is a {} step, not a {what} one", step.name, step.procedure.name()))); + } + Ok(step) +} + +/// The assembled matrices and the model's mass. +struct Assembled { + k: Csr, + m: Option, +} + +/// 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 { + let mut k = graph.zero_matrix(); + let mut m = mass.map(|_| graph.zero_matrix()); + for e in &mesh.elements { + let solid = e.solid()?; + let at = |err| Error::element(e.label, err); + add_block(&mut k, &e.nodes, &solid.stiffness(&e.d).map_err(at)?.data)?; + if let (Some(mm), Some(kind)) = (m.as_mut(), mass) { + let rho = e.density.ok_or_else(|| Error::Invalid(format!("element {} has no density", e.label)))?; + match kind { + MassKind::Consistent => add_scalar(mm, &e.nodes, &solid.mass(rho).map_err(at)?.data)?, + MassKind::Lumped => add_diagonal(mm, &e.nodes, &solid.lumped_mass(rho).map_err(at)?)?, + } + } + } + Ok(Assembled { k, m }) +} + +/// Factor a free-dof matrix; a singular pivot is reported as the free-dof column it hit. +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 and direction of a singular pivot. +fn factor_error(e: FactorError, mesh: &Mesh, map: &[Option]) -> 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)")) + } + FactorError::NotPositiveDefinite { column } => { + let dof = map.iter().position(|m| *m == Some(column)).unwrap_or(0); + let label = mesh.labels.get(dof / 3).copied().unwrap_or(0); + let axis = ["X", "Y", "Z"].get(dof % 3).copied().unwrap_or("?"); + Error::Singular(format!( + "the stiffness is singular at node {label} in {axis}: the model can move without straining (missing or insufficient boundary conditions, or a part not connected to the rest)" + )) + } + } +} + +/// Spread a free-equation vector over every degree of freedom (`base` holds the prescribed +/// values) and group it per node. +fn per_node(x: &[f64], map: &[Option], base: Vec) -> Vec<[f64; 3]> { + let mut full = base; + for (dof, m) in map.iter().enumerate() { + if let (Some(eq), Some(v)) = (m, full.get_mut(dof)) { + *v = x.get(*eq).copied().unwrap_or(0.0); + } + } + full.as_chunks::<3>().0.to_vec() +} + +/// Solve a linear static step: the step named `step`, or the model's only static step. +pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Result { + let resolved = model.resolve()?; + let step = select_step(model, step, |p| matches!(p, Procedure::Static {}), "static")?; + let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; + let fixed = prescribed(model, &resolved, &mesh, step, true)?; + let f = loads(model, &resolved, &mesh, step)?; + let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); + let k = assemble(&mesh, &graph, None)?.k; + let (map, neq) = equations(&nested_dissection(&mesh.xyz, &graph.adjacency), &fixed); + if neq == 0 { + return Err(Error::Invalid("every degree of freedom is prescribed; there is nothing to solve".into())); + } + let chol = factor_raw(&k.restrict(&map, neq), limits).map_err(|e| factor_error(e, &mesh, &map))?; + // Right-hand side: f_free − K_fp u_p (K u_p only when something moves). + let up: Vec = fixed.iter().map(|v| v.unwrap_or(0.0)).collect(); + let mut kup = vec![0.0; mesh.dofs()]; + if up.iter().any(|v| *v != 0.0) { + k.mul(&up, &mut kup); + } + let mut rhs = vec![0.0; neq]; + for (dof, m) in map.iter().enumerate() { + if let (Some(eq), Some(fv), Some(kv)) = (m, f.get(dof), kup.get(dof)) + && let Some(r) = rhs.get_mut(*eq) + { + *r = fv - kv; + } + } + chol.solve(&mut rhs)?; + let displacement = per_node(&rhs, &map, up); + if displacement.iter().flatten().any(|x| !x.is_finite()) { + return Err(Error::Singular("the solution is not finite; the model is ill-conditioned".into())); + } + // Reactions R = K u − f at the prescribed degrees of freedom (only their rows of K). + let mut reaction = vec![[0.0; 3]; mesh.labels.len()]; + let (mut applied_total, mut reaction_total) = ([0.0; 3], [0.0; 3]); + for (dof, fv) in f.iter().enumerate() { + let (node, c) = (dof / 3, dof % 3); + applied_total[c] += fv; + if fixed.get(dof).is_some_and(Option::is_some) { + let (cols, vals) = k.row(dof); + let ku: f64 = cols.iter().zip(vals).map(|(j, v)| v * displacement.get(j / 3).map_or(0.0, |u| u[j % 3])).sum(); + if let Some(rn) = reaction.get_mut(node) { + rn[c] = ku - fv; + } + reaction_total[c] += ku - fv; + } + } + // Stresses: one extrapolation per element type and rule. + let mut extrapolations: Vec<(ElementKind, 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 e in &mesh.elements { + let solid = e.solid()?; + let at = |err| Error::element(e.label, err); + let ue: Vec<[f64; 3]> = e.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); + let points = solid.stresses(&e.d, &ue).map_err(at)?; + let values: Vec<[f64; 6]> = points.iter().map(|p| p.stress).collect(); + let rule = solid.rule(); + let cached = extrapolations.iter().position(|(k, r, _)| *k == e.kind && *r == rule); + let index = match cached { + Some(i) => i, + None => { + extrapolations.push((e.kind, rule, extrapolation(solid.ty, rule).map_err(at)?)); + extrapolations.len() - 1 + } + }; + let nodal = extrapolations.get(index).map(|(_, _, x)| x.apply(&values)).unwrap_or_default(); + for (n, s) in e.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: e.label, kind: e.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 mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); + let scale = applied_total.iter().chain(&reaction_total).fold(0.0f64, |m, v| m.max(v.abs())); + let imbalance = (0..3).map(|c| (applied_total[c] + reaction_total[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, + reaction, + stress, + von_mises, + elements, + applied_total, + reaction_total, + stats: Stats { dofs: mesh.dofs(), equations: neq, factor_entries: chol.entries(), lanczos_vectors: 0 }, + warnings, + }) +} + +struct ShiftInvert<'a> { + m: &'a Csr, + factor: &'a Factor, +} + +impl Operator for ShiftInvert<'_> { + fn n(&self) -> usize { + self.m.n + } + fn mass(&self, x: &[f64], y: &mut [f64]) { + self.m.mul(x, y); + } + fn solve(&self, x: &mut [f64]) -> Result<()> { + self.factor.solve(x) + } +} + +/// Block size of the Lanczos iteration: enough for the six rigid-body modes of a free body. +const BLOCK: usize = 6; + +/// Relative residual of an accepted Ritz pair. +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. +pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOptions, limits: &Limits) -> Result { + let resolved: Resolved = model.resolve()?; + let step = select_step(model, step, |p| matches!(p, Procedure::Frequency { .. }), "frequency")?; + let Procedure::Frequency { modes, min_hz, max_hz } = step.procedure else { + return Err(Error::Internal("frequency step".into())); + }; + if let Some(s) = options.shift + && !s.is_finite() + { + return Err(Error::Invalid(format!("the shift must be a finite number, got {s}"))); + } + let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; + let fixed = prescribed(model, &resolved, &mesh, step, false)?; + let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); + let Assembled { k, m } = assemble(&mesh, &graph, Some(options.mass))?; + let m = m.ok_or_else(|| Error::Internal("mass matrix".into()))?; + // Every element's mass entries sum to its mass, once per direction. + let total_mass = m.val.iter().sum::() / 3.0; + let (map, neq) = equations(&nested_dissection(&mesh.xyz, &graph.adjacency), &fixed); + if neq == 0 { + return Err(Error::Invalid("every degree of freedom is held; there is nothing to vibrate".into())); + } + let kff = k.restrict(&map, neq); + let mff = m.restrict(&map, neq); + drop((k, m)); + let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); + let shifted = |s: f64| -> Result { factor_raw(&kff.add_scaled(-s, &mff)?, limits).map_err(|e| factor_error(e, &mesh, &map)) }; + let (shift, chol) = match options.shift { + Some(s) => (s, shifted(s)?), + None => match factor_raw(&kff, limits) { + Ok(f) => (0.0, f), + Err(FactorError::NotPositiveDefinite { .. }) => { + // A free (or partly free) body: shift below zero so the rigid-body modes come out. + let (tk, tm): (f64, f64) = (kff.diagonal().iter().sum(), mff.diagonal().iter().sum()); + if !(tk > 0.0 && tm > 0.0) { + return Err(Error::Singular("the stiffness or mass has no positive diagonal".into())); + } + let s = -1e-6 * tk / tm; + warnings.push("the structure is not fully supported; rigid-body modes appear at (near) zero frequency".into()); + (s, shifted(s)?) + } + Err(e) => return Err(factor_error(e, &mesh, &map)), + }, + }; + drop(kff); + let two_pi = 2.0 * std::f64::consts::PI; + let wanted = Wanted { + count: usize::try_from(modes).unwrap_or(usize::MAX).min(neq), + 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 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 dofs. + let mr: Vec> = (0..3) + .map(|d| { + let mut r = vec![0.0; neq]; + for (dof, mm) in map.iter().enumerate() { + if let Some(x) = mm.filter(|_| dof % 3 == d).and_then(|eq| r.get_mut(eq)) { + *x = 1.0; + } + } + let mut y = vec![0.0; neq]; + mff.mul(&r, &mut y); + y + }) + .collect(); + let mut out = Vec::with_capacity(pairs.values.len()); + for (lambda, phi) in pairs.values.iter().zip(&pairs.vectors) { + let participation = [0, 1, 2].map(|d| mr.get(d).map_or(0.0, |y| y.iter().zip(phi).map(|(a, b)| a * b).sum())); + out.push(Mode { + eigenvalue: *lambda, + frequency: lambda.max(0.0).sqrt() / two_pi, + shape: per_node(phi, &map, vec![0.0; mesh.dofs()]), + participation, + effective_mass: participation.map(|g| g * g), + }); + } + if out.len() < wanted.count && max_hz.is_none() { + warnings.push(format!("found {} of the {modes} modes asked for", out.len())); + } + Ok(FrequencyResult { + step: step.name.clone(), + nodes: mesh.labels.clone(), + modes: out, + total_mass, + stats: Stats { dofs: mesh.dofs(), equations: neq, factor_entries: chol.entries(), lanczos_vectors: pairs.basis }, + warnings, + }) +} + +/// The result of one step, as a session keeps it. +#[derive(Debug, Clone, PartialEq)] +pub enum StepResult { + Static(StaticResult), + Frequency(FrequencyResult), +} + +impl StepResult { + pub fn step(&self) -> &str { + match self { + StepResult::Static(r) => &r.step, + StepResult::Frequency(r) => &r.step, + } + } + + pub fn nodes(&self) -> &[u32] { + match self { + StepResult::Static(r) => &r.nodes, + StepResult::Frequency(r) => &r.nodes, + } + } +} diff --git a/crates/solver/src/cholesky.rs b/crates/solver/src/cholesky.rs new file mode 100644 index 0000000..506ad85 --- /dev/null +++ b/crates/solver/src/cholesky.rs @@ -0,0 +1,326 @@ +//! Sparse Cholesky factorization A = L Lᵀ of a symmetric positive definite matrix. +//! +//! The up-looking method of T. A. Davis, *Direct Methods for Sparse Linear Systems*, SIAM 2006, +//! §4.1–4.8: the elimination tree (J. W. H. Liu, *The role of elimination trees in sparse +//! factorization*, SIAM J. Matrix Anal. Appl. 11(1):134–172, 1990) gives the pattern of each row +//! of L as a reach in the tree, the symbolic pass counts every column before any value is +//! computed (so the factor's size is checked against a cap before memory is committed), and each +//! row of L is a sparse triangular solve. Operations run in a fixed order, with no threads and no +//! data-dependent reduction order, so the factor is the same bit for bit on every platform. +//! +//! A pivot at or below `PIVOT_TOLERANCE` times the matrix's own diagonal entry means the matrix +//! is singular to working precision (an unconstrained rigid-body motion, a disconnected part, a +//! zero-stiffness direction); the factorization stops and names the column. + +use crate::{Error, Result}; + +/// Relative pivot size below which the matrix counts as singular. +pub const PIVOT_TOLERANCE: f64 = 1e-10; + +/// The upper triangle of a symmetric matrix by columns: column k holds A(i, k) for i ≤ k, rows +/// sorted. +#[derive(Debug, Clone)] +pub struct UpperColumns { + pub n: usize, + pub ptr: Vec, + pub idx: Vec, + pub val: Vec, +} + +impl UpperColumns { + /// From a symmetric matrix stored by full rows (row k's entries with column ≤ k are column k + /// of the upper triangle). + pub fn from_symmetric(a: &crate::sparse::Csr) -> UpperColumns { + let mut ptr = Vec::with_capacity(a.n + 1); + ptr.push(0); + let mut idx = Vec::new(); + let mut val = Vec::new(); + for k in 0..a.n { + let (cols, vals) = a.row(k); + for (c, v) in cols.iter().zip(vals) { + if *c <= k { + idx.push(*c); + val.push(*v); + } + } + ptr.push(idx.len()); + } + UpperColumns { n: a.n, ptr, idx, val } + } + + fn column(&self, k: usize) -> (&[usize], &[f64]) { + let (a, b) = (self.ptr.get(k).copied().unwrap_or(0), self.ptr.get(k + 1).copied().unwrap_or(0)); + (self.idx.get(a..b).unwrap_or_default(), self.val.get(a..b).unwrap_or_default()) + } +} + +const NONE: usize = usize::MAX; + +/// The elimination tree: parent of each column, `NONE` at roots. +fn etree(a: &UpperColumns) -> Vec { + let mut parent = vec![NONE; a.n]; + let mut ancestor = vec![NONE; a.n]; + for k in 0..a.n { + let (rows, _) = a.column(k); + for &i0 in rows { + let mut i = i0; + // Walk from i up to the root of its current subtree, compressing the path to k. + while i != NONE && i < k { + let next = ancestor.get(i).copied().unwrap_or(NONE); + if let Some(x) = ancestor.get_mut(i) { + *x = k; + } + if next == NONE + && let Some(p) = parent.get_mut(i) + { + *p = k; + } + i = next; + } + } + } + parent +} + +/// The pattern of row k of L (excluding the diagonal): the reach of column k's entries in the +/// elimination tree. Fills `stack[top..]` in topological order and returns `top`; `mark` holds +/// `k` for visited columns. +fn ereach(a: &UpperColumns, k: usize, parent: &[usize], stack: &mut [usize], mark: &mut [usize], path: &mut Vec) -> usize { + let n = a.n; + let mut top = n; + if let Some(m) = mark.get_mut(k) { + *m = k; + } + let (rows, _) = a.column(k); + for &i0 in rows { + let mut i = i0; + if i > k { + continue; + } + path.clear(); + while i != NONE && mark.get(i).copied() != Some(k) { + path.push(i); + if let Some(m) = mark.get_mut(i) { + *m = k; + } + i = parent.get(i).copied().unwrap_or(NONE); + } + while let Some(v) = path.pop() { + top -= 1; + if let Some(s) = stack.get_mut(top) { + *s = v; + } + } + } + top +} + +/// A Cholesky factor, L by columns with the diagonal first in each column. +#[derive(Debug, Clone)] +pub struct Factor { + pub n: usize, + ptr: Vec, + idx: Vec, + val: Vec, +} + +/// Why a factorization stopped. +#[derive(Debug, Clone, PartialEq)] +pub enum FactorError { + /// The pivot of this column is not positive (relative to the diagonal). + NotPositiveDefinite { column: usize }, + /// The factor would hold more entries than allowed. + TooLarge { entries: usize, limit: usize }, +} + +impl Factor { + /// The elimination tree and the entry count of each column of L, without computing values. + fn symbolic(a: &UpperColumns) -> (Vec, Vec) { + let parent = etree(a); + let mut counts = vec![1usize; a.n]; + let mut stack = vec![0; a.n]; + let mut mark = vec![NONE; a.n]; + let mut path = Vec::new(); + for k in 0..a.n { + let top = ereach(a, k, &parent, &mut stack, &mut mark, &mut path); + for &j in stack.get(top..).unwrap_or_default() { + if let Some(c) = counts.get_mut(j) { + *c += 1; + } + } + } + (parent, counts) + } + + pub fn new(a: &UpperColumns, max_entries: usize) -> std::result::Result { + let n = a.n; + let (parent, counts) = Factor::symbolic(a); + let total: usize = counts.iter().sum(); + if total > max_entries { + return Err(FactorError::TooLarge { entries: total, limit: max_entries }); + } + let mut ptr = Vec::with_capacity(n + 1); + ptr.push(0); + for c in &counts { + ptr.push(ptr.last().copied().unwrap_or(0) + c); + } + let mut idx = vec![0usize; total]; + let mut val = vec![0.0f64; total]; + // Next free slot in each column. + let mut next: Vec = ptr.get(..n).unwrap_or_default().to_vec(); + let mut x = vec![0.0f64; n]; + let mut stack = vec![0usize; n]; + let mut mark = vec![NONE; n]; + let mut path = Vec::new(); + for k in 0..n { + let top = ereach(a, k, &parent, &mut stack, &mut mark, &mut path); + let (rows, vals) = a.column(k); + let mut diag = 0.0; + for (i, v) in rows.iter().zip(vals) { + if *i < k { + if let Some(xi) = x.get_mut(*i) { + *xi = *v; + } + } else if *i == k { + diag = *v; + } + } + let mut d = diag; + for &i in stack.get(top..).unwrap_or_default() { + let start = ptr.get(i).copied().unwrap_or(0); + let lii = val.get(start).copied().unwrap_or(1.0); + let lki = x.get(i).copied().unwrap_or(0.0) / lii; + if let Some(xi) = x.get_mut(i) { + *xi = 0.0; + } + let end = next.get(i).copied().unwrap_or(start); + for p in start + 1..end { + let (r, l) = (idx.get(p).copied().unwrap_or(0), val.get(p).copied().unwrap_or(0.0)); + if let Some(xr) = x.get_mut(r) { + *xr -= l * lki; + } + } + d -= lki * lki; + if let Some(slot) = next.get_mut(i) { + if let (Some(ri), Some(vi)) = (idx.get_mut(*slot), val.get_mut(*slot)) { + *ri = k; + *vi = lki; + } + *slot += 1; + } + } + if !(d > PIVOT_TOLERANCE * diag.abs()) || !d.is_finite() { + return Err(FactorError::NotPositiveDefinite { column: k }); + } + if let Some(slot) = next.get_mut(k) { + if let (Some(ri), Some(vi)) = (idx.get_mut(*slot), val.get_mut(*slot)) { + *ri = k; + *vi = d.sqrt(); + } + *slot += 1; + } + } + Ok(Factor { n, ptr, idx, val }) + } + + /// Entries in L. + pub fn entries(&self) -> usize { + self.val.len() + } + + /// Solve A x = b in place. + pub fn solve(&self, b: &mut [f64]) -> Result<()> { + if b.len() != self.n { + return Err(Error::Internal("right-hand side length".into())); + } + // L y = b, by columns. + for j in 0..self.n { + let (start, end) = (self.ptr.get(j).copied().unwrap_or(0), self.ptr.get(j + 1).copied().unwrap_or(0)); + let yj = b.get(j).copied().unwrap_or(0.0) / self.val.get(start).copied().unwrap_or(1.0); + if let Some(x) = b.get_mut(j) { + *x = yj; + } + for p in start + 1..end { + let (r, l) = (self.idx.get(p).copied().unwrap_or(0), self.val.get(p).copied().unwrap_or(0.0)); + if let Some(x) = b.get_mut(r) { + *x -= l * yj; + } + } + } + // Lᵀ x = y. + for j in (0..self.n).rev() { + let (start, end) = (self.ptr.get(j).copied().unwrap_or(0), self.ptr.get(j + 1).copied().unwrap_or(0)); + let mut s = b.get(j).copied().unwrap_or(0.0); + for p in start + 1..end { + let (r, l) = (self.idx.get(p).copied().unwrap_or(0), self.val.get(p).copied().unwrap_or(0.0)); + s -= l * b.get(r).copied().unwrap_or(0.0); + } + if let Some(x) = b.get_mut(j) { + *x = s / self.val.get(start).copied().unwrap_or(1.0); + } + } + Ok(()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::sparse::Csr; + + fn dense_to_csr(a: &[Vec]) -> Csr { + let n = a.len(); + let mut ptr = vec![0]; + let mut idx = Vec::new(); + let mut val = Vec::new(); + for row in a { + for (j, v) in row.iter().enumerate() { + if *v != 0.0 { + idx.push(j); + val.push(*v); + } + } + ptr.push(idx.len()); + } + Csr { n, ptr, idx, val } + } + + #[test] + fn factor_and_solve_a_sparse_spd_matrix() { + // A 2D Laplacian-like matrix on a 7 × 5 grid plus a diagonal shift: sparse with fill. + let (w, h) = (7, 5); + let n = w * h; + let mut a = vec![vec![0.0; n]; n]; + for i in 0..w { + for j in 0..h { + let k = i * h + j; + a[k][k] = 4.5; + if i + 1 < w { + a[k][k + h] = -1.0; + a[k + h][k] = -1.0; + } + if j + 1 < h { + a[k][k + 1] = -1.0; + a[k + 1][k] = -1.0; + } + } + } + let csr = dense_to_csr(&a); + let f = Factor::new(&UpperColumns::from_symmetric(&csr), usize::MAX).unwrap(); + let x: Vec = (0..n).map(|i| (i as f64 * 0.37).sin()).collect(); + let mut b = vec![0.0; n]; + csr.mul(&x, &mut b); + f.solve(&mut b).unwrap(); + assert!(b.iter().zip(&x).all(|(p, q)| (p - q).abs() < 1e-13), "{b:?}"); + // The cap is checked before the numeric factorization. + assert!(matches!(Factor::new(&UpperColumns::from_symmetric(&csr), 10), Err(FactorError::TooLarge { .. }))); + // A singular matrix (a pure Laplacian has the constant vector in its null space). + let mut sing = a.clone(); + for (k, row) in sing.iter_mut().enumerate() { + let off: f64 = row.iter().enumerate().filter(|(j, _)| *j != k).map(|(_, v)| -v).sum(); + row[k] = off; + } + let e = Factor::new(&UpperColumns::from_symmetric(&dense_to_csr(&sing)), usize::MAX).unwrap_err(); + assert_eq!(e, FactorError::NotPositiveDefinite { column: n - 1 }); + } +} diff --git a/crates/solver/src/eigen.rs b/crates/solver/src/eigen.rs new file mode 100644 index 0000000..263d44f --- /dev/null +++ b/crates/solver/src/eigen.rs @@ -0,0 +1,521 @@ +//! Eigenvalues: the dense symmetric eigenproblem of the Lanczos projection, and the shifted +//! block Lanczos method for the lowest modes of K φ = λ M φ. +//! +//! **Dense.** Householder reduction to tridiagonal form, then the implicit symmetric QR +//! algorithm with Wilkinson shifts, accumulating the transformations (G. H. Golub, C. F. Van +//! Loan, *Matrix Computations*, 4th ed., Johns Hopkins 2013, §8.3, algorithms 8.3.1 and 8.3.3). +//! +//! **Sparse.** The spectral transformation (T. Ericsson, A. Ruhe, *The spectral transformation +//! Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue +//! problems*, Math. Comp. 35(152):1251–1268, 1980): with a shift σ below the wanted +//! eigenvalues, the operator (K − σM)⁻¹M is symmetric in the M inner product and has eigenvalues +//! θ = 1/(λ − σ), so the lowest λ become the largest, well separated θ. A block of vectors +//! (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) +//! finds repeated eigenvalues up to the block size, such as the six rigid-body modes of a free +//! body. Every new block is reorthogonalized against all earlier ones, twice (full +//! reorthogonalization, B. N. Parlett, *The Symmetric Eigenvalue Problem*, SIAM 1998, ch. 13), +//! so no spurious copies appear. A Ritz pair has converged when its residual bound +//! ‖B_{j+1} s‖ is below a relative tolerance. The start block is pseudo-random from a fixed seed, +//! so results are reproducible. + +use crate::{Error, Result}; + +/// Eigenvalues (ascending) and eigenvectors (column k of the row-major `vectors`) of a dense +/// symmetric matrix (row-major, n × n). +pub fn symmetric_eigen(a: &[f64], n: usize) -> Result<(Vec, Vec)> { + if a.len() != n * n { + return Err(Error::Internal("dense eigenproblem size".into())); + } + let mut t = a.to_vec(); + let mut q = vec![0.0; n * n]; + for i in 0..n { + q[i * n + i] = 1.0; + } + // Householder tridiagonalization: for each column k, reflect A[k+1.., k] onto e1. + let mut v = vec![0.0; n]; + let mut p = vec![0.0; n]; + for k in 0..n.saturating_sub(2) { + let m = k + 1; + let norm: f64 = (m..n).map(|i| t[i * n + k] * t[i * n + k]).sum::().sqrt(); + if norm == 0.0 { + continue; + } + let x0 = t[m * n + k]; + let alpha = if x0 >= 0.0 { -norm } else { norm }; + for i in m..n { + v[i] = t[i * n + k]; + } + v[m] -= alpha; + let vnorm: f64 = (m..n).map(|i| v[i] * v[i]).sum::().sqrt(); + if vnorm == 0.0 { + continue; + } + for i in m..n { + v[i] /= vnorm; + } + // p = A22 v, w = p − (vᵀp) v, A22 −= 2 (v wᵀ + w vᵀ). + for i in m..n { + p[i] = (m..n).map(|j| t[i * n + j] * v[j]).sum(); + } + let kk: f64 = (m..n).map(|i| v[i] * p[i]).sum(); + for i in m..n { + p[i] -= kk * v[i]; + } + for i in m..n { + for j in m..n { + t[i * n + j] -= 2.0 * (v[i] * p[j] + p[i] * v[j]); + } + } + t[m * n + k] = alpha; + t[k * n + m] = alpha; + for i in m + 1..n { + t[i * n + k] = 0.0; + t[k * n + i] = 0.0; + } + // Q = Q (I − 2 v vᵀ). + for r in 0..n { + let s: f64 = (m..n).map(|j| q[r * n + j] * v[j]).sum(); + for j in m..n { + q[r * n + j] -= 2.0 * s * v[j]; + } + } + } + let mut d: Vec = (0..n).map(|i| t[i * n + i]).collect(); + let mut e: Vec = (0..n.saturating_sub(1)).map(|i| t[(i + 1) * n + i]).collect(); + tridiagonal_qr(&mut d, &mut e, &mut q, n)?; + // Sort ascending, carrying the vectors. + let mut order: Vec = (0..n).collect(); + order.sort_by(|&a, &b| d[a].total_cmp(&d[b]).then(a.cmp(&b))); + let values = order.iter().map(|&i| d[i]).collect(); + let mut vectors = vec![0.0; n * n]; + for (new, &old) in order.iter().enumerate() { + for r in 0..n { + vectors[r * n + new] = q[r * n + old]; + } + } + Ok((values, vectors)) +} + +/// Diagonalize a symmetric tridiagonal matrix (diagonal d, sub-diagonal e) by implicit QR steps +/// with Wilkinson shifts, applying the rotations to the columns of q. +fn tridiagonal_qr(d: &mut [f64], e: &mut [f64], q: &mut [f64], n: usize) -> Result<()> { + let small = |e: f64, a: f64, b: f64| e.abs() <= f64::EPSILON * (a.abs() + b.abs()); + let mut h = n.saturating_sub(1); + let mut iterations = 0usize; + let cap = 30 * n.max(1) + 30; + while h > 0 { + if small(e[h - 1], d[h - 1], d[h]) { + e[h - 1] = 0.0; + h -= 1; + continue; + } + let mut l = h - 1; + while l > 0 && !small(e[l - 1], d[l - 1], d[l]) { + l -= 1; + } + if l > 0 { + e[l - 1] = 0.0; + } + iterations += 1; + if iterations > cap { + return Err(Error::NotConverged("the tridiagonal QR iteration did not converge".into())); + } + // Wilkinson shift from the trailing 2 × 2 block. + let dd = (d[h - 1] - d[h]) / 2.0; + let eh = e[h - 1]; + let sign = if dd >= 0.0 { 1.0 } else { -1.0 }; + let mu = d[h] - eh * eh / (dd + sign * dd.hypot(eh)); + let mut x = d[l] - mu; + let mut z = e[l]; + for k in l..h { + let r = x.hypot(z); + let (c, s) = if r == 0.0 { (1.0, 0.0) } else { (x / r, -z / r) }; + if k > l { + e[k - 1] = c * e[k - 1] - s * z; + } + let (a, b, dq) = (d[k], e[k], d[k + 1]); + d[k] = c * c * a - 2.0 * c * s * b + s * s * dq; + d[k + 1] = s * s * a + 2.0 * c * s * b + c * c * dq; + e[k] = c * s * (a - dq) + (c * c - s * s) * b; + if k + 1 < h { + z = -s * e[k + 1]; + e[k + 1] *= c; + x = e[k]; + } + for row in 0..n { + let (qp, qq) = (q[row * n + k], q[row * n + k + 1]); + q[row * n + k] = c * qp - s * qq; + q[row * n + k + 1] = s * qp + c * qq; + } + } + } + Ok(()) +} + +/// A small deterministic pseudo-random generator (xorshift64*, S. Vigna, *An experimental +/// exploration of Marsaglia's xorshift generators, scrambled*, ACM TOMS 42(4), 2016) for +/// Lanczos start vectors. +pub struct Random(u64); + +impl Random { + pub fn new(seed: u64) -> Random { + Random(seed.max(1)) + } + + /// Uniform in [−1, 1). + pub fn uniform(&mut self) -> f64 { + self.0 ^= self.0 >> 12; + self.0 ^= self.0 << 25; + self.0 ^= self.0 >> 27; + let v = self.0.wrapping_mul(0x2545_f491_4f6c_dd1d); + (v >> 11) as f64 / (1u64 << 52) as f64 - 1.0 + } +} + +/// What the Lanczos iteration looks for, in terms of λ = σ + 1/θ. +#[derive(Debug, Clone, Copy)] +pub struct Wanted { + /// Modes to return within the range. + pub count: usize, + /// Lower and upper bounds on λ (eigenvalues outside are skipped). + pub min: f64, + pub max: f64, +} + +/// The operators the Lanczos iteration needs. +pub trait Operator { + fn n(&self) -> usize; + /// y = M x. + fn mass(&self, x: &[f64], y: &mut [f64]); + /// x ← (K − σM)⁻¹ x. + fn solve(&self, x: &mut [f64]) -> Result<()>; +} + +/// Converged eigenpairs: λ ascending, vectors M-orthonormal. +#[derive(Debug, Clone)] +pub struct Eigenpairs { + pub values: Vec, + pub vectors: Vec>, + /// Lanczos vectors used. + pub basis: usize, +} + +fn dot(a: &[f64], b: &[f64]) -> f64 { + a.iter().zip(b).map(|(x, y)| x * y).sum() +} + +/// Orthonormalize a block against the basis and itself in the M inner product (classical +/// Gram–Schmidt, twice). A column that vanishes is replaced by a fresh random direction and its +/// diagonal coefficient stays zero. Returns the new vectors, their M products and the +/// upper-triangular coefficients `r[j][i]` of column i on new vector j. +fn orthonormalize( + op: &impl Operator, + w: Vec>, + q: &[Vec], + mq: &[Vec], + rng: &mut Random, + scratch: &mut [f64], +) -> (Vec>, Vec>, Vec>) { + let n = op.n(); + let k = w.len(); + let mut r = vec![vec![0.0; k]; k]; + let mut out: Vec> = Vec::with_capacity(k); + let mut mout: Vec> = Vec::with_capacity(k); + for (i, mut v) in w.into_iter().enumerate() { + let mut fresh = false; + for _attempt in 0..4 { + op.mass(&v, scratch); + let before = dot(&v, scratch).max(0.0).sqrt(); + for _ in 0..2 { + for (qj, mqj) in q.iter().zip(mq) { + let c = dot(mqj, &v); + v.iter_mut().zip(qj).for_each(|(x, y)| *x -= c * y); + } + for (j, (qj, mqj)) in out.iter().zip(&mout).enumerate() { + let c = dot(mqj, &v); + v.iter_mut().zip(qj).for_each(|(x, y)| *x -= c * y); + if !fresh && let Some(x) = r.get_mut(j).and_then(|row| row.get_mut(i)) { + *x += c; + } + } + } + op.mass(&v, scratch); + let norm = dot(&v, scratch).max(0.0).sqrt(); + if norm > 1e-8 * before && norm > 0.0 { + if !fresh && let Some(x) = r.get_mut(i).and_then(|row| row.get_mut(i)) { + *x = norm; + } + v.iter_mut().for_each(|x| *x /= norm); + mout.push(scratch.iter().map(|x| x / norm).collect()); + out.push(v); + break; + } + fresh = true; + v = (0..n).map(|_| rng.uniform()).collect(); + } + } + (out, mout, r) +} + +/// Shift-invert block Lanczos with full reorthogonalization. +pub fn lanczos(op: &impl Operator, shift: f64, wanted: Wanted, block: usize, max_basis: usize, tolerance: f64) -> Result { + let n = op.n(); + if n == 0 { + return Err(Error::Invalid("there are no free degrees of freedom".into())); + } + let b = block.clamp(1, n); + let max_basis = max_basis.min(n).max(b); + let mut rng = Random::new(0x9e37_79b9_7f4a_7c15); + let mut q: Vec> = Vec::new(); + let mut mq: Vec> = Vec::new(); + // The projected matrix T, grown block by block (row-major in a growing square). + let mut t: Vec> = Vec::new(); + let mut scratch = vec![0.0; n]; + + // Start block. + let start: Vec> = (0..b).map(|_| (0..n).map(|_| rng.uniform()).collect()).collect(); + let (mut cur, mut mcur, _) = orthonormalize(op, start, &q, &mq, &mut rng, &mut scratch); + if cur.len() < b { + return Err(Error::NotConverged("could not build a Lanczos start block".into())); + } + let mut prev_b: Option>> = None; + let mut next_check = (wanted.count + b).max(2 * b); + loop { + let j0 = q.len(); + // W = Op Q_j. + let mut w: Vec> = Vec::with_capacity(b); + for mx in &mcur { + let mut y = mx.clone(); + op.solve(&mut y)?; + w.push(y); + } + // A_j = Q_jᵀ M W (symmetrized), W −= Q_j A_j (the earlier blocks go in the full + // reorthogonalization below). + let mut a = vec![vec![0.0; b]; b]; + for (i, mqi) in mcur.iter().enumerate() { + for (j, wj) in w.iter().enumerate() { + a[i][j] = dot(mqi, wj); + } + } + for i in 0..b { + for j in i + 1..b { + let s = 0.5 * (a[i][j] + a[j][i]); + a[i][j] = s; + a[j][i] = s; + } + } + q.append(&mut cur); + mq.append(&mut mcur); + let m = q.len(); + for row in &mut t { + row.resize(m, 0.0); + } + t.resize(m, vec![0.0; m]); + for i in 0..b { + for j in 0..b { + t[j0 + i][j0 + j] = a[i][j]; + } + } + if let Some(bp) = &prev_b { + // T[block j][block j−1] = B_j, T[block j−1][block j] = B_jᵀ. + for i in 0..b { + for j in 0..b { + t[j0 + i][j0 - b + j] = bp[i][j]; + t[j0 - b + j][j0 + i] = bp[i][j]; + } + } + } + // Full reorthogonalization happens inside the QR against the whole basis. + let full = m >= max_basis || m + b > n; + let (next, mnext, bnext) = + if full { (Vec::new(), Vec::new(), vec![vec![0.0; b]; b]) } else { orthonormalize(op, w, &q, &mq, &mut rng, &mut scratch) }; + // Convergence check on a schedule (the dense eigenproblem costs m³). + if m >= next_check || full || next.len() < b { + next_check = m + (m / 4).max(b); + let dense: Vec = t.iter().flatten().copied().collect(); + let (theta, s) = symmetric_eigen(&dense, m)?; + // Order by λ ascending = θ descending (θ > 0 below the spectrum). + let mut idx: Vec = (0..m).collect(); + idx.sort_by(|&x, &y| theta[y].total_cmp(&theta[x]).then(x.cmp(&y))); + let mut done = false; + let mut accepted = Vec::new(); + let mut found = 0; + for &i in &idx { + let th = theta[i]; + if !(th > 0.0) { + break; + } + // Residual bound ‖B_{j+1} s_last‖ with s_last the last block of the Ritz vector. + let mut res2 = 0.0; + if !full { + for r in 0..b { + let v: f64 = (0..b).map(|c| bnext[r][c] * s[(m - b + c) * m + i]).sum(); + res2 += v * v; + } + } + if res2.sqrt() > tolerance * th { + break; + } + let lambda = shift + 1.0 / th; + if lambda > wanted.max { + done = true; + break; + } + if lambda >= wanted.min { + accepted.push(i); + found += 1; + if found == wanted.count { + done = true; + break; + } + } + } + if full && !done { + if m >= n { + done = true; + } else { + return Err(Error::NotConverged(format!( + "found {found} of {} modes after {m} Lanczos vectors (the cap); raise the iteration limit or ask for fewer modes", + wanted.count + ))); + } + } + if done { + let mut values = Vec::new(); + let mut vectors = Vec::new(); + for &i in &accepted { + values.push(shift + 1.0 / theta[i]); + let mut x = vec![0.0; n]; + for (k, qk) in q.iter().enumerate() { + let c = s[k * m + i]; + x.iter_mut().zip(qk).for_each(|(a, b)| *a += c * b); + } + op.mass(&x, &mut scratch); + let norm = dot(&x, &scratch).max(0.0).sqrt(); + if norm > 0.0 { + x.iter_mut().for_each(|a| *a /= norm); + } + // Deterministic sign: the largest component positive. + let big = x.iter().copied().fold(0.0f64, |m, v| if v.abs() > m.abs() { v } else { m }); + if big < 0.0 { + x.iter_mut().for_each(|a| *a = -*a); + } + vectors.push(x); + } + return Ok(Eigenpairs { values, vectors, basis: m }); + } + } + if full { + return Err(Error::NotConverged("the Lanczos basis is exhausted".into())); + } + cur = next; + mcur = mnext; + if cur.len() < b { + return Err(Error::NotConverged("the Lanczos block lost rank".into())); + } + prev_b = Some(bnext); + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn dense_symmetric_eigen() { + let n = 9; + let mut rng = Random::new(7); + let mut a = vec![0.0; n * n]; + for i in 0..n { + for j in 0..=i { + let v = rng.uniform(); + a[i * n + j] = v; + a[j * n + i] = v; + } + } + let (values, vectors) = symmetric_eigen(&a, n).unwrap(); + assert!(values.windows(2).all(|w| w[0] <= w[1])); + for k in 0..n { + for i in 0..n { + let av: f64 = (0..n).map(|j| a[i * n + j] * vectors[j * n + k]).sum(); + assert!((av - values[k] * vectors[i * n + k]).abs() < 1e-12, "pair {k}"); + } + let norm: f64 = (0..n).map(|i| vectors[i * n + k].powi(2)).sum(); + assert!((norm - 1.0).abs() < 1e-12); + } + assert_eq!(symmetric_eigen(&[3.0], 1).unwrap().0, vec![3.0]); + } + + /// K and M diagonal. + struct Diagonal { + k: Vec, + m: Vec, + shift: f64, + } + + impl Operator for Diagonal { + fn n(&self) -> usize { + self.k.len() + } + fn mass(&self, x: &[f64], y: &mut [f64]) { + for i in 0..x.len() { + y[i] = self.m[i] * x[i]; + } + } + fn solve(&self, x: &mut [f64]) -> Result<()> { + for i in 0..x.len() { + x[i] /= self.k[i] - self.shift * self.m[i]; + } + Ok(()) + } + } + + #[test] + fn lanczos_finds_repeated_and_zero_eigenvalues() { + // Six zero eigenvalues, then 1 twice, then distinct values; M not the identity. + let n = 200; + let lambdas: Vec = (0..n) + .map(|i| { + if i < 6 { + 0.0 + } else if i < 8 { + 1.0 + } else { + 1.0 + (i - 7) as f64 * 0.5 + } + }) + .collect(); + let m: Vec = (0..n).map(|i| 1.0 + (i % 7) as f64 * 0.1).collect(); + let k: Vec = lambdas.iter().zip(&m).map(|(l, m)| l * m).collect(); + let shift = -1e-3; + let op = Diagonal { k, m, shift }; + let wanted = Wanted { count: 10, min: f64::NEG_INFINITY, max: f64::INFINITY }; + let pairs = lanczos(&op, shift, wanted, 6, n, 1e-9).unwrap(); + assert_eq!(pairs.values.len(), 10); + let want = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.5, 2.0]; + for (got, w) in pairs.values.iter().zip(want) { + assert!((got - w).abs() < 1e-8, "{:?}", pairs.values); + } + // M-orthonormal vectors. + let mut y = vec![0.0; n]; + for a in &pairs.vectors { + for b in &pairs.vectors { + op.mass(b, &mut y); + let d = dot(a, &y); + assert!((d - if std::ptr::eq(a, b) { 1.0 } else { 0.0 }).abs() < 1e-8); + } + } + // A range: values in [1.2, 3.2] (λ), at most 100. + let ranged = lanczos(&op, shift, Wanted { count: 100, min: 1.2, max: 3.2 }, 6, n, 1e-9).unwrap(); + assert_eq!(ranged.values.len(), 4, "{:?}", ranged.values); + assert!(ranged.values.iter().zip([1.5, 2.0, 2.5, 3.0]).all(|(g, w)| (g - w).abs() < 1e-8), "{:?}", ranged.values); + // Caps: a basis too small to converge is an error, not a hang. + assert!(lanczos(&op, shift, Wanted { count: 50, min: f64::NEG_INFINITY, max: f64::INFINITY }, 6, 24, 1e-9).is_err()); + // Tiny problems: the whole space. + let small = Diagonal { k: vec![2.0, 1.0, 3.0], m: vec![1.0; 3], shift: 0.0 }; + let all = lanczos(&small, 0.0, Wanted { count: 5, min: f64::NEG_INFINITY, max: f64::INFINITY }, 6, 100, 1e-9).unwrap(); + assert_eq!(all.values.len(), 3); + assert!((all.values[0] - 1.0).abs() < 1e-12 && (all.values[2] - 3.0).abs() < 1e-12, "{:?}", all.values); + } +} diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs new file mode 100644 index 0000000..9e7e452 --- /dev/null +++ b/crates/solver/src/lib.rs @@ -0,0 +1,89 @@ +//! ModelCraft's implicit solver: linear static and frequency analysis of a model's steps. +//! +//! The pipeline, for one step of a [`modelcraft_model::Model`]: +//! 1. **System** ([`system`]): the model is resolved; every section's elements become solid +//! elements with their material's elasticity and the section's formulation +//! ([`modelcraft_fem`] kernels); the step's boundary conditions give the prescribed degrees of +//! freedom, its loads the load vector (forces and tractions as the model distributes them, +//! consistent pressure loads, gravity and body forces). +//! 2. **Assembly** ([`sparse`]) on the node graph, elements in label order. +//! 3. **Ordering** ([`ordering`]): nested dissection of the nodes; the free degrees of freedom +//! are numbered in that order and the prescribed ones eliminated (u_p moved to the right-hand +//! side), with no penalty numbers. +//! 4. **Solution** ([`cholesky`]): sparse Cholesky; reactions K u − f at the prescribed degrees +//! of freedom; stresses at the integration points, extrapolated and averaged at the nodes, and +//! von Mises. Frequency analysis ([`eigen`]) runs shift-invert block Lanczos on the same +//! kind of factorization. +//! +//! Everything is sequential and in a fixed order, so the same model gives the same numbers bit +//! for bit; sizes and iterations are capped by [`Limits`]. +#![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] +#![forbid(unsafe_code)] +// Matrix kernels read as the textbook formulas with index loops. +#![allow(clippy::needless_range_loop)] + +pub mod analysis; +pub mod cholesky; +pub mod eigen; +pub mod ordering; +pub mod sparse; +pub mod system; + +pub use analysis::{ + ElementStress, FrequencyOptions, FrequencyResult, MassKind, Mode, StaticResult, Stats, StepResult, solve_frequency, solve_static, +}; + +#[derive(Debug, Clone, PartialEq, thiserror::Error)] +pub enum Error { + #[error(transparent)] + Model(#[from] modelcraft_model::Error), + #[error("invalid input: {0}")] + Invalid(String), + #[error("not supported yet: {0}")] + NotSupported(String), + /// An element kernel failed (inverted geometry, a formulation it lacks, …). + #[error("element {label}: {source}")] + Element { label: u32, source: modelcraft_fem::Error }, + /// The stiffness (or shifted matrix) is singular. + #[error("{0}")] + Singular(String), + #[error("limit reached: {0}")] + Limit(String), + #[error("did not converge: {0}")] + NotConverged(String), + /// A broken internal invariant (a bug), reported instead of a panic. + #[error("internal error: {0}")] + Internal(String), +} + +pub type Result = std::result::Result; + +impl Error { + /// A kernel error, naming the element. + pub fn element(label: u32, source: modelcraft_fem::Error) -> Error { + Error::Element { label, source } + } +} + +/// Caps on the size of a solution. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct Limits { + /// Degrees of freedom of the mesh. + pub max_dofs: usize, + /// Entries of the Cholesky factor (8 bytes each); also bounds the Lanczos basis memory. + pub max_factor_entries: usize, + /// Lanczos vectors (iterations × block size). + pub max_lanczos_vectors: 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; +} + +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 } + } +} diff --git a/crates/solver/src/ordering.rs b/crates/solver/src/ordering.rs new file mode 100644 index 0000000..4e04340 --- /dev/null +++ b/crates/solver/src/ordering.rs @@ -0,0 +1,119 @@ +//! Fill-reducing ordering: nested dissection by coordinate bisection. +//! +//! Nested dissection (A. George, *Nested dissection of a regular finite element mesh*, SIAM J. +//! Numer. Anal. 10(2):345–363, 1973) orders a separator after the two parts it splits, so the +//! Cholesky factor of a 3D mesh with n nodes fills O(n^{4/3}) instead of the O(n^{5/3}) of a band +//! or envelope ordering. The separator here is geometric: the nodes are split at the median of +//! the longest axis of their bounding box, and the nodes of the smaller side that touch the other +//! side form the separator (the coordinate nested dissection of H. D. Simon, *Partitioning of +//! unstructured problems for parallel processing*, Computing Systems in Engineering 2:135–148, +//! 1991, with vertex separators). Ties are broken by node index, so the ordering is +//! deterministic. + +/// Parts this small are ordered as they come. +const LEAF: usize = 16; + +/// A permutation of the nodes: `order[k]` is the node eliminated k-th. +pub fn nested_dissection(xyz: &[[f64; 3]], adjacency: &[Vec]) -> Vec { + let n = xyz.len(); + let mut order = Vec::with_capacity(n); + // 0 = not in the current part, 1 = left, 2 = right. + let mut side = vec![0u8; n]; + // An explicit stack instead of recursion: (part, separator to emit after it). + enum Task { + Split(Vec), + Emit(Vec), + } + let mut stack = vec![Task::Split((0..n).collect())]; + while let Some(task) = stack.pop() { + let mut part = match task { + Task::Emit(sep) => { + order.extend(sep); + continue; + } + Task::Split(p) => p, + }; + if part.len() <= LEAF { + order.extend(part); + continue; + } + let mut lo = [f64::INFINITY; 3]; + let mut hi = [f64::NEG_INFINITY; 3]; + for &v in &part { + let p = xyz.get(v).copied().unwrap_or_default(); + for k in 0..3 { + lo[k] = lo[k].min(p[k]); + hi[k] = hi[k].max(p[k]); + } + } + let axis = (0..3).fold(0, |best, k| if hi[k] - lo[k] > hi[best] - lo[best] { k } else { best }); + let key = |v: usize| xyz.get(v).map_or(0.0, |p| p[axis]); + part.sort_by(|&a, &b| key(a).total_cmp(&key(b)).then(a.cmp(&b))); + let half = part.len() / 2; + for (i, &v) in part.iter().enumerate() { + if let Some(s) = side.get_mut(v) { + *s = if i < half { 1 } else { 2 }; + } + } + let touches = |v: usize, other: u8| adjacency.get(v).is_some_and(|row| row.iter().any(|w| side.get(*w) == Some(&other))); + let (left, right) = part.split_at(half); + let left_boundary: Vec = left.iter().copied().filter(|&v| touches(v, 2)).collect(); + let right_boundary: Vec = right.iter().copied().filter(|&v| touches(v, 1)).collect(); + let separator = if left_boundary.len() <= right_boundary.len() { left_boundary } else { right_boundary }; + for &v in &part { + if let Some(s) = side.get_mut(v) { + *s = 0; + } + } + // Mark the separator so the parts leave it out. + for &v in &separator { + if let Some(s) = side.get_mut(v) { + *s = 3; + } + } + let a: Vec = left.iter().copied().filter(|v| side.get(*v) != Some(&3)).collect(); + let b: Vec = right.iter().copied().filter(|v| side.get(*v) != Some(&3)).collect(); + for &v in &separator { + if let Some(s) = side.get_mut(v) { + *s = 0; + } + } + // Order: a, then b, then the separator (the stack runs in reverse). + stack.push(Task::Emit(separator)); + stack.push(Task::Split(b)); + stack.push(Task::Split(a)); + } + order +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn a_permutation_with_separators_last() { + // A 12 × 12 grid graph. + let m = 12; + let id = |i: usize, j: usize| i * m + j; + let xyz: Vec<[f64; 3]> = (0..m * m).map(|v| [(v / m) as f64, (v % m) as f64, 0.0]).collect(); + let mut adj = vec![Vec::new(); m * m]; + for i in 0..m { + for j in 0..m { + for (di, dj) in [(0i64, 1i64), (1, 0), (0, -1), (-1, 0), (0, 0)] { + let (a, b) = (i as i64 + di, j as i64 + dj); + if (0..m as i64).contains(&a) && (0..m as i64).contains(&b) { + adj[id(i, j)].push(id(a as usize, b as usize)); + } + } + } + } + let order = nested_dissection(&xyz, &adj); + let mut sorted = order.clone(); + sorted.sort_unstable(); + assert_eq!(sorted, (0..m * m).collect::>()); + // The top separator (last 12 entries) is one grid line. + let last: Vec = order[m * m - m..].to_vec(); + assert!(last.iter().all(|v| v / m == last[0] / m) || last.iter().all(|v| v % m == last[0] % m), "{last:?}"); + assert_eq!(order, nested_dissection(&xyz, &adj), "deterministic"); + } +} diff --git a/crates/solver/src/sparse.rs b/crates/solver/src/sparse.rs new file mode 100644 index 0000000..bd6a5ca --- /dev/null +++ b/crates/solver/src/sparse.rs @@ -0,0 +1,233 @@ +//! Sparse symmetric matrices in compressed rows, both triangles stored. +//! +//! Assembly works on the node graph: a node's row holds every node it shares an element with, +//! and each node pair is a dense 3 × 3 block, so the pattern is built once from the element +//! connectivity and the stiffness and mass share it. Sums run in a fixed order (elements by +//! label, then the element's own node order), so assembled values are the same bit for bit on +//! every run and platform. + +use crate::{Error, Result}; + +/// A square sparse matrix in compressed rows with sorted column indices. +#[derive(Debug, Clone, PartialEq)] +pub struct Csr { + pub n: usize, + pub ptr: Vec, + pub idx: Vec, + pub val: Vec, +} + +impl Csr { + /// The row's column indices and values. + pub fn row(&self, r: usize) -> (&[usize], &[f64]) { + let (a, b) = (self.ptr.get(r).copied().unwrap_or(0), self.ptr.get(r + 1).copied().unwrap_or(0)); + (self.idx.get(a..b).unwrap_or_default(), self.val.get(a..b).unwrap_or_default()) + } + + /// y = A x. + pub fn mul(&self, x: &[f64], y: &mut [f64]) { + for (r, yr) in y.iter_mut().enumerate().take(self.n) { + let (cols, vals) = self.row(r); + *yr = cols.iter().zip(vals).map(|(c, v)| v * x.get(*c).copied().unwrap_or(0.0)).sum(); + } + } + + pub fn diagonal(&self) -> Vec { + (0..self.n) + .map(|r| { + let (cols, vals) = self.row(r); + cols.binary_search(&r).ok().and_then(|i| vals.get(i)).copied().unwrap_or(0.0) + }) + .collect() + } + + /// The matrix on a subset of rows and columns, renumbered: `map[i]` is the new index of old + /// index i, or `None` to drop it. The new rows keep sorted columns. + pub fn restrict(&self, map: &[Option], n: usize) -> Csr { + let mut rows: Vec> = vec![Vec::new(); n]; + for (r, m) in map.iter().enumerate() { + let Some(nr) = *m else { continue }; + let (cols, vals) = self.row(r); + let Some(out) = rows.get_mut(nr) else { continue }; + for (c, v) in cols.iter().zip(vals) { + if let Some(Some(nc)) = map.get(*c) { + out.push((*nc, *v)); + } + } + out.sort_by_key(|e| e.0); + } + let mut ptr = Vec::with_capacity(n + 1); + ptr.push(0); + let mut idx = Vec::new(); + let mut val = Vec::new(); + for row in rows { + for (c, v) in row { + idx.push(c); + val.push(v); + } + ptr.push(idx.len()); + } + Csr { n, ptr, idx, val } + } + + /// `self + a·other` on the same pattern. + pub fn add_scaled(&self, a: f64, other: &Csr) -> Result { + if self.ptr != other.ptr || self.idx != other.idx { + return Err(Error::Internal("matrices with different patterns".into())); + } + let val = self.val.iter().zip(&other.val).map(|(x, y)| x + a * y).collect(); + Ok(Csr { n: self.n, ptr: self.ptr.clone(), idx: self.idx.clone(), val }) + } +} + +/// The node-block pattern of a mesh: for each node, the sorted nodes it shares an element with +/// (itself included). +#[derive(Debug, Clone)] +pub struct NodeGraph { + pub adjacency: Vec>, +} + +impl NodeGraph { + pub fn new(nodes: usize, elements: impl Iterator>) -> NodeGraph { + let mut adjacency: Vec> = vec![Vec::new(); nodes]; + for e in elements { + let e = e.as_ref(); + for &a in e { + if let Some(row) = adjacency.get_mut(a) { + row.extend_from_slice(e); + } + } + } + for row in &mut adjacency { + row.sort_unstable(); + row.dedup(); + } + NodeGraph { adjacency } + } + + /// Nonzeros of the dof-level matrix (3 × 3 per node pair). + pub fn entries(&self) -> usize { + self.adjacency.iter().map(|r| 9 * r.len()).sum() + } + + /// A zero matrix over every degree of freedom (3 per node) on this pattern. + pub fn zero_matrix(&self) -> Csr { + let n = 3 * self.adjacency.len(); + let mut ptr = Vec::with_capacity(n + 1); + ptr.push(0); + let mut idx = Vec::with_capacity(self.entries()); + for row in &self.adjacency { + for _ in 0..3 { + for &b in row { + idx.extend([3 * b, 3 * b + 1, 3 * b + 2]); + } + ptr.push(idx.len()); + } + } + let val = vec![0.0; idx.len()]; + Csr { n, ptr, idx, val } + } +} + +/// Where node b's 3 × 3 block starts in each of node a's three rows (the rows of a node share +/// one pattern, built by [`NodeGraph::zero_matrix`]): one search per node pair. +fn block_offsets(k: &Csr, a: usize, b: usize) -> Result<[usize; 3]> { + let (start, end) = (k.ptr.get(3 * a).copied().unwrap_or(0), k.ptr.get(3 * a + 1).copied().unwrap_or(0)); + let cols = k.idx.get(start..end).ok_or_else(|| Error::Internal("pattern".into()))?; + let at = cols.binary_search(&(3 * b)).map_err(|_| Error::Internal("entry outside the pattern".into()))?; + Ok([0, 1, 2].map(|r| k.ptr.get(3 * a + r).copied().unwrap_or(0) + at)) +} + +/// Add an element matrix (3n × 3n, row-major, degrees of freedom node by node) into a matrix +/// built by [`NodeGraph::zero_matrix`]; `nodes` are the element's global node indices. +pub fn add_block(k: &mut Csr, nodes: &[usize], data: &[f64]) -> Result<()> { + let width = 3 * nodes.len(); + if data.len() != width * width { + return Err(Error::Internal("element matrix size".into())); + } + for (i, &a) in nodes.iter().enumerate() { + for (j, &b) in nodes.iter().enumerate() { + let rows = block_offsets(k, a, b)?; + for (r, start) in rows.iter().enumerate() { + for c in 0..3 { + if let (Some(x), Some(v)) = (k.val.get_mut(start + c), data.get((3 * i + r) * width + 3 * j + c)) { + *x += v; + } + } + } + } + } + Ok(()) +} + +/// Add an n × n element matrix to each displacement component's block (consistent mass). +pub fn add_scalar(k: &mut Csr, nodes: &[usize], data: &[f64]) -> Result<()> { + let n = nodes.len(); + if data.len() != n * n { + return Err(Error::Internal("element matrix size".into())); + } + for (i, &a) in nodes.iter().enumerate() { + for (j, &b) in nodes.iter().enumerate() { + let rows = block_offsets(k, a, b)?; + let v = data.get(i * n + j).copied().unwrap_or(0.0); + for (c, start) in rows.iter().enumerate() { + if let Some(x) = k.val.get_mut(start + c) { + *x += v; + } + } + } + } + Ok(()) +} + +/// Add a per-node diagonal to each displacement component (lumped mass). +pub fn add_diagonal(k: &mut Csr, nodes: &[usize], diagonal: &[f64]) -> Result<()> { + for (&a, v) in nodes.iter().zip(diagonal) { + let rows = block_offsets(k, a, a)?; + for (c, start) in rows.iter().enumerate() { + if let Some(x) = k.val.get_mut(start + c) { + *x += v; + } + } + } + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn assemble_restrict_multiply() { + // Two "elements" sharing node 1. + let g = NodeGraph::new(3, [vec![0usize, 1], vec![1, 2]].into_iter()); + assert_eq!(g.adjacency, vec![vec![0, 1], vec![0, 1, 2], vec![1, 2]]); + let mut k = g.zero_matrix(); + let mut e = vec![0.0; 36]; + for i in 0..6 { + e[i * 6 + i] = 2.0; + e[i * 6 + (i + 3) % 6] = -1.0; + } + add_block(&mut k, &[0, 1], &e).unwrap(); + add_block(&mut k, &[1, 2], &e).unwrap(); + assert_eq!(k.diagonal(), vec![2.0, 2.0, 2.0, 4.0, 4.0, 4.0, 2.0, 2.0, 2.0]); + let x = vec![1.0; 9]; + let mut y = vec![0.0; 9]; + k.mul(&x, &mut y); + assert_eq!(y, vec![1.0; 3].into_iter().chain([2.0; 3]).chain([1.0; 3]).collect::>()); + // Keep the x components only, in reverse order. + let map: Vec> = (0..9).map(|i| (i % 3 == 0).then_some(2 - i / 3)).collect(); + let r = k.restrict(&map, 3); + assert_eq!(r.diagonal(), vec![2.0, 4.0, 2.0]); + assert_eq!(r.row(1).0, &[0, 1, 2]); + assert!(add_block(&mut k, &[0, 2], &e).is_err(), "nodes 0 and 2 share no element"); + // Scalar (per component) and diagonal additions land on each component. + let mut m = g.zero_matrix(); + add_scalar(&mut m, &[0, 1], &[2.0, 1.0, 1.0, 2.0]).unwrap(); + add_diagonal(&mut m, &[2], &[5.0]).unwrap(); + assert_eq!(m.diagonal(), vec![2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 5.0, 5.0, 5.0]); + assert_eq!(m.row(0).1, &[2.0, 0.0, 0.0, 1.0, 0.0, 0.0]); + let s = k.add_scaled(-1.0, &k).unwrap(); + assert!(s.val.iter().all(|v| *v == 0.0)); + } +} diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs new file mode 100644 index 0000000..567606c --- /dev/null +++ b/crates/solver/src/system.rs @@ -0,0 +1,218 @@ +//! From the resolved model to the discrete system: the solid elements with their materials and +//! formulations, the degrees of freedom, prescribed values and the load vector of a step. + +use modelcraft_fem::element::{Formulation, Solid}; +use modelcraft_fem::{Elasticity, elasticity, face}; +use modelcraft_model::{ElementKind, Load, Model, Resolved, SectionKind, Step, Target}; + +use crate::{Error, Result}; + +/// A solid element ready for the kernels. +#[derive(Debug, Clone)] +pub struct Elem { + pub label: u32, + pub kind: ElementKind, + /// Indices into [`Mesh::labels`]. + pub nodes: Vec, + pub xyz: Vec<[f64; 3]>, + pub formulation: Formulation, + pub d: Elasticity, + pub density: Option, +} + +impl Elem { + pub fn solid(&self) -> Result> { + Solid::new(self.kind, &self.xyz, self.formulation).map_err(|e| Error::element(self.label, e)) + } +} + +/// The analysed mesh: the nodes that carry degrees of freedom (those of the solid elements), +/// sorted by label, and the elements sorted by label. +#[derive(Debug, Clone)] +pub struct Mesh { + pub labels: Vec, + pub xyz: Vec<[f64; 3]>, + pub elements: Vec, +} + +impl Mesh { + pub fn node(&self, label: u32) -> Option { + self.labels.binary_search(&label).ok() + } + + pub fn element(&self, label: u32) -> Option<&Elem> { + self.elements.binary_search_by_key(&label, |e| e.label).ok().and_then(|i| self.elements.get(i)) + } + + pub fn dofs(&self) -> usize { + 3 * self.labels.len() + } + + /// Build it from the model and its resolved view, refusing what the solver cannot represent. + pub fn new(model: &Model, resolved: &Resolved, max_dofs: usize) -> Result { + if let Some(name) = model.couplings().keys().next() { + return Err(Error::NotSupported(format!( + "couplings (coupling `{name}`); the solver handles solid elements without reference points so far" + ))); + } + let mut elements = Vec::new(); + for (name, section) in model.sections() { + if !matches!(section.kind, SectionKind::Solid {}) { + return Err(Error::NotSupported(format!( + "{} sections (section `{name}`); the solver has solid elements so far", + section.kind.name() + ))); + } + let material = model + .material(§ion.material) + .ok_or_else(|| Error::Invalid(format!("section `{name}`: material `{}` is missing", section.material)))?; + let elastic = + material.elastic.as_ref().ok_or_else(|| Error::Invalid(format!("material `{}` has no elastic constants", section.material)))?; + let d = elasticity(elastic).map_err(|e| Error::Invalid(format!("material `{}`: {e}", section.material)))?; + let set = resolved.sets.get(§ion.set).ok_or_else(|| Error::Invalid(format!("set `{}` is not resolved", section.set)))?; + for &label in &set.elements { + let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; + let formulation = Formulation::from_section(el.kind, §ion.formulation).map_err(|e| Error::element(label, e))?; + elements.push((label, el.kind, el.nodes.clone(), formulation, d, material.density)); + } + } + elements.sort_by_key(|e| e.0); + let mut labels: Vec = elements.iter().flat_map(|e| e.2.iter().copied()).collect(); + labels.sort_unstable(); + labels.dedup(); + if 3 * labels.len() > max_dofs { + return Err(Error::Limit(format!("the model has {} degrees of freedom; the limit is {max_dofs}", 3 * labels.len()))); + } + let xyz: Vec<[f64; 3]> = labels + .iter() + .map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("node {l} is missing")))) + .collect::>()?; + let mut mesh = Mesh { labels, xyz, elements: Vec::with_capacity(elements.len()) }; + for (label, kind, nodes, formulation, d, density) 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 }); + } + Ok(mesh) + } +} + +/// The prescribed value of every degree of freedom (`None` = free) under a step's boundary +/// conditions. Rotational degrees of freedom do not exist on solid nodes and are ignored. +pub fn prescribed(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, values: bool) -> Result>> { + let mut out: Vec> = vec![None; mesh.dofs()]; + for name in &step.boundary_conditions { + let bc = model.boundary_condition(name).ok_or_else(|| Error::Invalid(format!("boundary condition `{name}` is missing")))?; + if let Target::ReferencePoint(rp) = bc.target() { + return Err(Error::NotSupported(format!("boundary conditions on reference points (`{name}` on `{rp}`)"))); + } + for &label in resolved.target_nodes(bc.target()) { + let Some(node) = mesh.node(label) else { continue }; + for dof in 1..=3u8 { + let Some(v) = bc.value(dof) else { continue }; + let v = if values { v } else { 0.0 }; + let slot = out.get_mut(3 * node + usize::from(dof - 1)).ok_or_else(|| Error::Internal("dof index".into()))?; + match *slot { + Some(old) if old != v => { + return Err(Error::Invalid(format!( + "node {label} direction {dof} is prescribed both {old} and {v} (boundary condition `{name}` and another)" + ))); + } + _ => *slot = Some(v), + } + } + } + } + Ok(out) +} + +fn add(f: &mut [f64], node: usize, v: [f64; 3]) { + for (k, c) in v.iter().enumerate() { + if let Some(x) = f.get_mut(3 * node + k) { + *x += c; + } + } +} + +/// The load vector of a static step: nodal forces and tractions as the model distributes them, +/// consistent pressure loads, gravity and body forces. +pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step) -> Result> { + let mut f = vec![0.0; mesh.dofs()]; + for name in &step.loads { + let load = model.load(name).ok_or_else(|| Error::Invalid(format!("load `{name}` is missing")))?; + match load { + Load::Force { target: Target::ReferencePoint(rp), .. } => { + return Err(Error::NotSupported(format!("loads on reference points (`{name}` on `{rp}`)"))); + } + Load::Force { .. } | Load::Traction { .. } => { + for (label, v) in resolved.nodal_loads.get(name).map(Vec::as_slice).unwrap_or_default() { + let node = mesh.node(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: node {label} belongs to no solid element")))?; + add(&mut f, node, *v); + } + } + Load::Moment { .. } => { + return Err(Error::NotSupported(format!("moments (`{name}`): solid nodes have no rotations; a coupled reference point is needed"))); + } + Load::Pressure { set, magnitude } => { + let faces = resolved.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default(); + for fc in faces { + let labels = model + .face_nodes(*fc) + .ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} face {} is missing", fc.element, fc.face)))?; + 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: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); + let forces = face::pressure(&xyz, *magnitude).map_err(|e| Error::element(fc.element, e))?; + for (i, v) in idx.iter().zip(forces) { + add(&mut f, *i, v); + } + } + } + Load::Gravity { acceleration, set } => { + let chosen: Option<&[u32]> = set.as_ref().map(|s| resolved.sets.get(s).map(|r| r.elements.as_slice()).unwrap_or_default()); + for e in &mesh.elements { + if chosen.is_some_and(|c| c.binary_search(&e.label).is_err()) { + continue; + } + let rho = e.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", e.label)))?; + body(&mut f, e, acceleration.map(|a| rho * a))?; + } + } + Load::BodyForce { set, vector } => { + for label in resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default() { + let e = mesh.element(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: element {label} is not a solid element")))?; + body(&mut f, e, *vector)?; + } + } + } + } + Ok(f) +} + +fn body(f: &mut [f64], e: &Elem, b: [f64; 3]) -> Result<()> { + let forces = e.solid()?.body_force(b).map_err(|err| Error::element(e.label, err))?; + for (node, v) in e.nodes.iter().zip(forces) { + add(f, *node, v); + } + Ok(()) +} + +/// Equation numbers of the free degrees of freedom, following a node order. +pub fn equations(order: &[usize], prescribed: &[Option]) -> (Vec>, usize) { + let mut map = vec![None; prescribed.len()]; + let mut next = 0; + for &node in order { + for c in 0..3 { + let dof = 3 * node + c; + if prescribed.get(dof).is_some_and(Option::is_none) { + if let Some(m) = map.get_mut(dof) { + *m = Some(next); + } + next += 1; + } + } + } + (map, next) +} diff --git a/crates/solver/tests/common/mod.rs b/crates/solver/tests/common/mod.rs new file mode 100644 index 0000000..e9be812 --- /dev/null +++ b/crates/solver/tests/common/mod.rs @@ -0,0 +1,179 @@ +//! Structured meshes generated in code for the analytic verification suite: a block of +//! hexahedra (hex8, hex20) or Kuhn-split tetrahedra (tet4, tet10) on [0, 1]³, mapped to the part +//! by a function, and helpers to define the model on it. +#![allow(dead_code)] + +use std::collections::BTreeMap; + +use modelcraft_model::{ + BoundaryCondition, Elastic, ElementFace, ElementKind, Formulation, Integration, Load, Material, Model, Procedure, Section, SectionKind, Set, + Step, Symmetry, Target, +}; + +pub struct Grid { + pub model: Model, + /// Parametric coordinates of every node. + pub param: BTreeMap, + pub kind: ElementKind, +} + +/// A structured mesh with `n` cells per direction. Nodes live on the doubled lattice +/// (2n + 1 points per direction) so quadratic elements get their mid-side nodes from the map +/// too (curved edges follow it). +pub fn grid(n: [usize; 3], kind: ElementKind, map: impl Fn([f64; 3]) -> [f64; 3]) -> Grid { + let size = n.map(|k| 2 * k + 1); + let label = |i: [usize; 3]| (1 + i[0] + size[0] * (i[1] + size[1] * i[2])) as u32; + let param_of = |i: [usize; 3]| [0, 1, 2].map(|k| i[k] as f64 / (2 * n[k]) as f64); + let mut model = Model::default(); + let mut param = BTreeMap::new(); + let mut node = |model: &mut Model, i: [usize; 3]| -> u32 { + let l = label(i); + if let std::collections::btree_map::Entry::Vacant(slot) = param.entry(l) { + let p = param_of(i); + model.set_node(l, map(p)).unwrap(); + slot.insert(p); + } + l + }; + let mut next = 1u32; + for cz in 0..n[2] { + for cy in 0..n[1] { + for cx in 0..n[0] { + let base = [2 * cx, 2 * cy, 2 * cz]; + let at = |d: [usize; 3]| [base[0] + d[0], base[1] + d[1], base[2] + d[2]]; + match kind { + ElementKind::Hex8 | ElementKind::Hex20 => { + let corners = [[0, 0, 0], [2, 0, 0], [2, 2, 0], [0, 2, 0], [0, 0, 2], [2, 0, 2], [2, 2, 2], [0, 2, 2]]; + let mut nodes: Vec = corners.iter().map(|c| node(&mut model, at(*c))).collect(); + if kind == ElementKind::Hex20 { + for e in kind.edges() { + let (a, b) = (corners[e[0]], corners[e[1]]); + nodes.push(node(&mut model, at([0, 1, 2].map(|k| (a[k] + b[k]) / 2)))); + } + } + model.set_element(next, kind, nodes).unwrap(); + next += 1; + } + _ => { + // Six tetrahedra around the cell's main diagonal (Kuhn), conforming across cells. + for perm in [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]] { + let mut v = [[0usize; 3]; 4]; + for s in 0..3 { + v[s + 1] = v[s]; + v[s + 1][perm[s]] = 2; + } + // Orientation from the parametric volume. + let p = v.map(|c| c.map(|x| x as f64)); + let e = |i: usize| [0, 1, 2].map(|k| p[i][k] - p[0][k]); + let (a, b, c) = (e(1), e(2), e(3)); + let vol = a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]); + if vol < 0.0 { + v.swap(1, 2); + } + let mut nodes: Vec = v.iter().map(|c| node(&mut model, at(*c))).collect(); + if kind == ElementKind::Tet10 { + for e in kind.edges() { + let (a, b) = (v[e[0]], v[e[1]]); + nodes.push(node(&mut model, at([0, 1, 2].map(|k| (a[k] + b[k]) / 2)))); + } + } + model.set_element(next, kind, nodes).unwrap(); + next += 1; + } + } + } + } + } + } + Grid { model, param, kind } +} + +pub fn near(a: f64, b: f64) -> bool { + (a - b).abs() < 1e-9 +} + +impl Grid { + pub fn nodes_where(&self, f: impl Fn([f64; 3]) -> bool) -> Vec { + self.param.iter().filter(|(_, p)| f(**p)).map(|(l, _)| *l).collect() + } + + /// Element faces whose nodes all satisfy `f` (on parametric coordinates). + pub fn faces_where(&self, f: impl Fn([f64; 3]) -> bool) -> Vec { + let mut out = Vec::new(); + for (label, e) in self.model.elements() { + for (i, face) in e.kind.faces().iter().enumerate() { + if face.iter().all(|k| f(self.param[&e.nodes[*k]])) { + out.push(ElementFace { element: *label, face: i as u8 + 1 }); + } + } + } + out + } + + /// Material, section over every element and an empty step of the given procedure. + pub fn define(&mut self, elastic: Elastic, density: Option, formulation: Formulation, procedure: Procedure) { + let all: Vec = self.model.elements().keys().copied().collect(); + self.model.add_set("all", Set { elements: all, ..Set::default() }).unwrap(); + self.model.add_material("Mat", Material { elastic: Some(elastic), density }).unwrap(); + self.model.add_section("solid", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap(); + self.model.add_step(Step::new("Step", procedure), None).unwrap(); + } + + pub fn node_set(&mut self, name: &str, f: impl Fn([f64; 3]) -> bool) { + let nodes = self.nodes_where(f); + self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap(); + } + + pub fn face_set(&mut self, name: &str, f: impl Fn([f64; 3]) -> bool) { + let faces = self.faces_where(f); + self.model.add_set(name, Set { faces, ..Set::default() }).unwrap(); + } + + pub fn bc(&mut self, name: &str, bc: BoundaryCondition) { + self.model.add_boundary_condition(name, bc).unwrap(); + self.model.set_boundary_condition_active("Step", name, true).unwrap(); + } + + pub fn load(&mut self, name: &str, load: Load) { + self.model.add_load(name, load).unwrap(); + self.model.set_load_active("Step", name, true).unwrap(); + } + + pub fn symmetry(&mut self, name: &str, plane: Symmetry, f: impl Fn([f64; 3]) -> bool) { + self.node_set(name, f); + self.bc(name, BoundaryCondition::Symmetry { target: Target::Set(name.into()), plane }); + } + + pub fn fixed(&mut self, name: &str, f: impl Fn([f64; 3]) -> bool) { + self.node_set(name, f); + self.bc(name, BoundaryCondition::Fixed { target: Target::Set(name.into()) }); + } +} + +pub fn iso(young: f64, poisson: f64) -> Elastic { + Elastic::Isotropic { young, poisson } +} + +pub fn full() -> Formulation { + Formulation::default() +} + +pub fn reduced() -> Formulation { + Formulation { integration: Integration::Reduced, ..Formulation::default() } +} + +pub fn incompatible() -> Formulation { + Formulation { incompatible_modes: true, ..Formulation::default() } +} + +/// The element kinds and formulations under test. +pub fn variants() -> Vec<(ElementKind, Formulation, &'static str)> { + vec![ + (ElementKind::Tet4, full(), "tet4"), + (ElementKind::Tet10, full(), "tet10"), + (ElementKind::Hex8, full(), "hex8"), + (ElementKind::Hex8, incompatible(), "hex8 incompatible"), + (ElementKind::Hex20, full(), "hex20"), + (ElementKind::Hex20, reduced(), "hex20 reduced"), + ] +} diff --git a/crates/solver/tests/verification.rs b/crates/solver/tests/verification.rs new file mode 100644 index 0000000..029dab0 --- /dev/null +++ b/crates/solver/tests/verification.rs @@ -0,0 +1,367 @@ +//! The analytic verification suite: every element and formulation against closed-form +//! solutions. Run with `--nocapture` to see the errors. + +mod common; + +use std::f64::consts::PI; + +use common::*; +use modelcraft_model::{BoundaryCondition, ElementKind, Load, Procedure, Symmetry, Target}; +use modelcraft_solver::{FrequencyOptions, Limits, MassKind, StaticResult, solve_frequency, solve_static}; + +fn static_solve(g: &Grid) -> StaticResult { + solve_static(&g.model, Some("Step"), &Limits::default()).unwrap() +} + +fn at(r: &StaticResult, label: u32) -> usize { + r.nodes.binary_search(&label).unwrap() +} + +/// Displacement patch test (MacNeal & Harder 1985 style, on a distorted 2 × 2 × 2 patch): a +/// linear displacement prescribed on the boundary is reproduced exactly inside, with exactly +/// constant strain and stress in every element. +#[test] +fn patch_test_every_element() { + let a = [[1e-3, 2e-4, -3e-4], [5e-4, -2e-3, 1e-4], [-1e-4, 3e-4, 7e-4]]; + let (young, nu) = (1000.0, 0.25); + let lambda = young * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); + let mu = young / (2.0 * (1.0 + nu)); + let strain = [a[0][0], a[1][1], a[2][2], a[0][1] + a[1][0], a[0][2] + a[2][0], a[1][2] + a[2][1]]; + let tr = strain[0] + strain[1] + strain[2]; + let stress = [ + lambda * tr + 2.0 * mu * strain[0], + lambda * tr + 2.0 * mu * strain[1], + lambda * tr + 2.0 * mu * strain[2], + mu * strain[3], + mu * strain[4], + mu * strain[5], + ]; + let field = |x: [f64; 3]| [0, 1, 2].map(|i| (0..3).map(|j| a[i][j] * x[j]).sum::()); + // Interior distortion that leaves the cube's boundary in place. + let map = |p: [f64; 3]| { + let b = (PI * p[0]).sin() * (PI * p[1]).sin() * (PI * p[2]).sin(); + [p[0] + 0.12 * b, p[1] - 0.08 * b, p[2] + 0.1 * b] + }; + for (kind, formulation, name) in variants() { + let mut g = grid([2, 2, 2], kind, map); + g.define(iso(young, nu), None, formulation, Procedure::Static {}); + let on_boundary = |p: [f64; 3]| p.iter().any(|c| near(*c, 0.0) || near(*c, 1.0)); + for (i, label) in g.nodes_where(on_boundary).into_iter().enumerate() { + let x = g.model.node(label).unwrap().xyz; + let u = field(x); + let set = format!("b{i}"); + g.model.add_set(&set, modelcraft_model::Set { nodes: vec![label], ..Default::default() }).unwrap(); + g.bc( + &set, + BoundaryCondition::Displacement { target: Target::Set(set.clone()), values: [Some(u[0]), Some(u[1]), Some(u[2]), None, None, None] }, + ); + } + let r = static_solve(&g); + let mut worst: f64 = 0.0; + for (label, u) in r.nodes.iter().zip(&r.displacement) { + let want = field(g.model.node(*label).unwrap().xyz); + worst = worst.max((0..3).map(|k| (u[k] - want[k]).abs()).fold(0.0, f64::max)); + } + let mut worst_stress: f64 = 0.0; + for e in &r.elements { + for p in &e.points { + worst_stress = worst_stress.max((0..6).map(|k| (p.stress[k] - stress[k]).abs()).fold(0.0, f64::max)); + } + } + for s in &r.stress { + worst_stress = worst_stress.max((0..6).map(|k| (s[k] - stress[k]).abs()).fold(0.0, f64::max)); + } + eprintln!("patch {name}: displacement error {worst:.1e}, stress error {worst_stress:.1e}"); + assert!(worst < 1e-12, "{name}: {worst}"); + assert!(worst_stress < 1e-10, "{name}: {worst_stress}"); + } +} + +/// Uniaxial tension by a consistent surface load on a distorted block with symmetry planes: +/// σxx = p exactly, ux = p L / E, the lateral contraction −ν p y / E, and the reactions balance +/// the load (the force form of the patch test, and a check of pressure loads). +#[test] +fn uniaxial_tension_by_pressure() { + let (young, nu, p) = (2000.0, 0.3, 5.0); + let (lx, ly, lz) = (3.0, 2.0, 1.5); + for (kind, formulation, name) in variants() { + let mut g = grid([2, 2, 2], kind, |q| { + let b = (PI * q[0]).sin() * (PI * q[1]).sin() * (PI * q[2]).sin(); + [lx * q[0] + 0.2 * b, ly * q[1] + 0.1 * b, lz * q[2] - 0.1 * b] + }); + g.define(iso(young, nu), None, formulation, Procedure::Static {}); + g.symmetry("x0", Symmetry::X, |q| near(q[0], 0.0)); + g.symmetry("y0", Symmetry::Y, |q| near(q[1], 0.0)); + g.symmetry("z0", Symmetry::Z, |q| near(q[2], 0.0)); + g.face_set("end", |q| near(q[0], 1.0)); + // A negative pressure pulls. + g.load("pull", Load::Pressure { set: "end".into(), magnitude: -p }); + let r = static_solve(&g); + let mut worst: f64 = 0.0; + for (label, u) in r.nodes.iter().zip(&r.displacement) { + let x = g.model.node(*label).unwrap().xyz; + let want = [p * x[0] / young, -nu * p * x[1] / young, -nu * p * x[2] / young]; + worst = worst.max((0..3).map(|k| (u[k] - want[k]).abs()).fold(0.0, f64::max)); + } + let s_err = r.stress.iter().map(|s| (s[0] - p).abs().max(s[1].abs()).max(s[5].abs())).fold(0.0, f64::max); + eprintln!("tension {name}: displacement error {worst:.1e}, stress error {s_err:.1e}"); + assert!(worst < 1e-11, "{name}: {worst}"); + assert!(s_err < 1e-9, "{name}: {s_err}"); + assert!((r.applied_total[0] - p * ly * lz).abs() < 1e-9 && (r.reaction_total[0] + p * ly * lz).abs() < 1e-9, "{name}"); + assert!(r.von_mises.iter().all(|v| (v - p).abs() < 1e-9)); + } +} + +/// Gravity is ρ g V, held by the reactions. +#[test] +fn gravity_and_body_force_totals() { + let mut g = grid([3, 2, 2], ElementKind::Hex20, |q| [4.0 * q[0], 2.0 * q[1], q[2]]); + g.define(iso(1000.0, 0.3), Some(2.5), full(), Procedure::Static {}); + g.fixed("base", |q| near(q[2], 0.0)); + g.load("g", Load::Gravity { acceleration: [0.0, 0.0, -9.81], set: None }); + g.load("b", Load::BodyForce { set: "all".into(), vector: [1.0, 0.0, 0.0] }); + let r = static_solve(&g); + let v = 8.0; + assert!((r.applied_total[2] + 2.5 * 9.81 * v).abs() < 1e-9, "{:?}", r.applied_total); + assert!((r.applied_total[0] - v).abs() < 1e-9); + for k in 0..3 { + assert!((r.applied_total[k] + r.reaction_total[k]).abs() < 1e-8, "{k}"); + } +} + +/// Cantilever L × b × h, clamped at x = 0, a shear traction on the free end. Timoshenko +/// deflection δ = P L³ / (3 E I) + P L / (κ G A), κ = 5/6, with ν = 0 so the clamp adds no +/// lateral restraint. +fn cantilever_tip(kind: ElementKind, formulation: modelcraft_model::Formulation, n: [usize; 3]) -> f64 { + let (l, b, h, young, p) = (10.0, 1.0, 1.0, 1000.0, 1.0); + let mut g = grid(n, kind, |q| [l * q[0], b * q[1], h * q[2]]); + g.define(iso(young, 0.0), None, formulation, Procedure::Static {}); + g.fixed("root", |q| near(q[0], 0.0)); + g.face_set("end", |q| near(q[0], 1.0)); + g.load("shear", Load::Traction { set: "end".into(), vector: [0.0, 0.0, -p / (b * h)] }); + let r = static_solve(&g); + let end = g.nodes_where(|q| near(q[0], 1.0)); + -end.iter().map(|l| r.displacement[at(&r, *l)][2]).sum::() / end.len() as f64 +} + +fn timoshenko() -> f64 { + let (l, b, h, young, p) = (10.0f64, 1.0f64, 1.0f64, 1000.0, 1.0); + let i = b * h.powi(3) / 12.0; + p * l.powi(3) / (3.0 * young * i) + p * l / (5.0 / 6.0 * young / 2.0 * b * h) +} + +#[test] +fn cantilever_tip_deflection_and_convergence() { + let exact = timoshenko(); + let euler_bernoulli = 4.0; + eprintln!("cantilever: Timoshenko {exact:.5}, Euler–Bernoulli {euler_bernoulli:.5}"); + let err = |d: f64| (d - exact).abs() / exact; + let study = |kind, f, meshes: &[[usize; 3]], name: &str| -> Vec { + meshes + .iter() + .map(|n| { + let d = cantilever_tip(kind, f, *n); + eprintln!(" {name} {n:?}: {d:.5} (error {:.2}%)", 100.0 * err(d)); + err(d) + }) + .collect() + }; + // Quadratic and enhanced elements are accurate on coarse meshes. + let hex20 = study(ElementKind::Hex20, full(), &[[10, 1, 1], [20, 2, 2]], "hex20"); + assert!(hex20[1] < 0.01); + let hex20r = study(ElementKind::Hex20, reduced(), &[[10, 2, 2]], "hex20 reduced"); + assert!(hex20r[0] < 0.01); + let im = study(ElementKind::Hex8, incompatible(), &[[10, 1, 1], [20, 2, 2]], "hex8 incompatible"); + assert!(im[0] < 0.03 && im[1] < 0.015); + let tet10 = study(ElementKind::Tet10, full(), &[[10, 2, 2]], "tet10"); + assert!(tet10[0] < 0.02); + // Linear elements lock in bending and converge under refinement. + let meshes = [[5, 1, 1], [10, 2, 2], [20, 4, 4], [40, 8, 8]]; + let hex8 = study(ElementKind::Hex8, full(), &meshes, "hex8"); + let tet4 = study(ElementKind::Tet4, full(), &meshes, "tet4"); + for (name, e) in [("hex8", &hex8), ("tet4", &tet4)] { + let rates: Vec = e.windows(2).map(|w| (w[0] / w[1]).log2()).collect(); + eprintln!(" {name} convergence rates {rates:.2?}"); + assert!(e.windows(2).all(|w| w[1] < w[0]), "{name} converges monotonically"); + assert!(*rates.last().unwrap() > 1.5, "{name} approaches second order: {rates:?}"); + } + assert!(hex8[3] < 0.04 && tet4[3] < 0.12, "{hex8:?} {tet4:?}"); +} + +/// Lamé's thick cylinder under internal pressure, plane strain (u_z = 0 on both ends), quarter +/// model: σ_r = A − B/r², σ_θ = A + B/r², A = p a²/(b² − a²), B = p a² b²/(b² − a²), +/// u_r = (1 + ν)/E ((1 − 2ν) A r + B / r). +#[test] +fn lame_thick_cylinder() { + let (a, b, t, p, young, nu) = (1.0f64, 2.0f64, 0.25, 10.0, 1000.0, 0.3); + let big_a = p * a * a / (b * b - a * a); + let big_b = p * a * a * b * b / (b * b - a * a); + let ur = |r: f64| (1.0 + nu) / young * ((1.0 - 2.0 * nu) * big_a * r + big_b / r); + let sigma_t = |r: f64| big_a + big_b / (r * r); + for (kind, formulation, n, tol_u, tol_s, name) in [ + (ElementKind::Hex20, full(), [6, 8, 1], 0.002, 0.01, "hex20"), + (ElementKind::Hex8, incompatible(), [12, 16, 1], 0.01, 0.04, "hex8 incompatible"), + (ElementKind::Tet10, full(), [6, 8, 1], 0.005, 0.03, "tet10"), + ] { + let mut g = grid(n, kind, |q| { + let (r, th) = (a + (b - a) * q[0], PI / 2.0 * q[1]); + [r * th.cos(), r * th.sin(), t * q[2]] + }); + g.define(iso(young, nu), None, formulation, Procedure::Static {}); + g.symmetry("y0", Symmetry::Y, |q| near(q[1], 0.0)); + g.symmetry("x0", Symmetry::X, |q| near(q[1], 1.0)); + g.node_set("ends", |q| near(q[2], 0.0) || near(q[2], 1.0)); + g.bc("ends", BoundaryCondition::Symmetry { target: Target::Set("ends".into()), plane: Symmetry::Z }); + g.face_set("bore", |q| near(q[0], 0.0)); + g.load("p", Load::Pressure { set: "bore".into(), magnitude: p }); + let r = static_solve(&g); + let mut eu: f64 = 0.0; + for (label, u) in r.nodes.iter().zip(&r.displacement) { + let x = g.model.node(*label).unwrap().xyz; + let rad = x[0].hypot(x[1]); + let got = (u[0] * x[0] + u[1] * x[1]) / rad; + eu = eu.max((got - ur(rad)).abs() / ur(a)); + } + // Hoop stress at the bore (nodal average), on the 45° line to stay away from the + // symmetry edges. + let bore: Vec = g.nodes_where(|q| near(q[0], 0.0) && near(q[1], 0.5)); + let st: f64 = bore + .iter() + .map(|l| { + let s = r.stress[at(&r, *l)]; + let x = g.model.node(*l).unwrap().xyz; + let (c, sn) = (x[0] / x[0].hypot(x[1]), x[1] / x[0].hypot(x[1])); + // σ_θ = s_xx sin² − 2 s_xy sin cos + s_yy cos². + s[0] * sn * sn - 2.0 * s[3] * sn * c + s[1] * c * c + }) + .sum::() + / bore.len() as f64; + let es = (st - sigma_t(a)).abs() / sigma_t(a); + eprintln!("lamé {name} {n:?}: u_r error {:.3}%, bore hoop stress {st:.4} vs {:.4} ({:.2}%)", 100.0 * eu, sigma_t(a), 100.0 * es); + assert!(eu < tol_u && es < tol_s, "{name}"); + } +} + +/// Kirsch: a plate with a hole of radius a under remote tension σ₀, plane stress, quarter model +/// with the plate's half-width 20a. The hoop stress at the top of the hole tends to 3σ₀ (the +/// finite width adds well under 1%). +#[test] +fn kirsch_plate_with_hole() { + let (a, w, s0) = (1.0f64, 20.0f64, 1.0); + let ratio: f64 = 40.0; + let map = |q: [f64; 3]| { + let th = PI / 2.0 * q[1]; + let outer = if th <= PI / 4.0 { [w, w * th.tan()] } else { [w / th.tan(), w] }; + // Geometric grading toward the hole. + let s = (ratio.powf(q[0]) - 1.0) / (ratio - 1.0); + let inner = [a * th.cos(), a * th.sin()]; + [inner[0] + s * (outer[0] - inner[0]), inner[1] + s * (outer[1] - inner[1]), 0.1 * q[2]] + }; + let mut kt = Vec::new(); + for (kind, formulation, n, name) in [ + (ElementKind::Hex20, full(), [6, 8, 1], "hex20 coarse"), + (ElementKind::Hex20, full(), [12, 16, 1], "hex20 fine"), + (ElementKind::Hex8, incompatible(), [24, 32, 1], "hex8 incompatible"), + ] { + let mut g = grid(n, kind, map); + g.define(iso(1000.0, 0.3), None, formulation, Procedure::Static {}); + g.symmetry("x0", Symmetry::X, |q| near(q[1], 1.0)); + g.symmetry("y0", Symmetry::Y, |q| near(q[1], 0.0)); + g.symmetry("z0", Symmetry::Z, |q| near(q[2], 0.0)); + g.face_set("right", |q| near(q[0], 1.0) && q[1] <= 0.5 + 1e-9); + g.load("pull", Load::Traction { set: "right".into(), vector: [s0, 0.0, 0.0] }); + let r = static_solve(&g); + let top: Vec = g.nodes_where(|q| near(q[0], 0.0) && near(q[1], 1.0)); + let k = top.iter().map(|l| r.stress[at(&r, *l)][0]).sum::() / top.len() as f64 / s0; + eprintln!("kirsch {name} {n:?}: Kt = {k:.4} (error vs 3: {:.2}%)", 100.0 * (k - 3.0).abs() / 3.0); + kt.push(k); + } + assert!((kt[1] - 3.0).abs() < (kt[0] - 3.0).abs(), "refinement approaches 3: {kt:?}"); + assert!((kt[1] - 3.0).abs() / 3.0 < 0.02, "{kt:?}"); + assert!((kt[2] - 3.0).abs() / 3.0 < 0.04, "{kt:?}"); +} + +/// Euler–Bernoulli frequencies f = (βL)² / (2π L²) √(EI/(ρA)). +fn beam_frequency(beta_l: f64, l: f64, young: f64, rho: f64, side: f64) -> f64 { + let i = side.powi(4) / 12.0; + let area = side * side; + beta_l * beta_l / (2.0 * PI * l * l) * (young * i / (rho * area)).sqrt() +} + +fn beam(kind: ElementKind, n: [usize; 3], clamped: bool, modes: u32) -> Grid { + let l = 20.0; + let mut g = grid(n, kind, |q| [l * q[0], q[1], q[2]]); + g.define(iso(1000.0, 0.0), Some(1.0), full(), Procedure::Frequency { modes, min_hz: None, max_hz: None }); + if clamped { + g.fixed("root", |q| near(q[0], 0.0)); + } + g +} + +#[test] +fn cantilever_natural_frequencies() { + let f1 = beam_frequency(1.875_104_07, 20.0, 1000.0, 1.0, 1.0); + let f2 = beam_frequency(4.694_091_13, 20.0, 1000.0, 1.0, 1.0); + for (kind, n, mass, name) in [ + (ElementKind::Hex20, [20, 1, 1], MassKind::Consistent, "hex20"), + (ElementKind::Hex20, [20, 1, 1], MassKind::Lumped, "hex20 lumped"), + (ElementKind::Tet10, [20, 2, 2], MassKind::Consistent, "tet10"), + ] { + let g = beam(kind, n, true, 4); + let r = solve_frequency(&g.model, Some("Step"), &FrequencyOptions { mass, shift: None }, &Limits::default()).unwrap(); + let f: Vec = r.modes.iter().map(|m| m.frequency).collect(); + eprintln!("cantilever modes {name}: {f:.6?} vs EB {f1:.6}, {f2:.6} ({} Lanczos vectors)", r.stats.lanczos_vectors); + assert_eq!(f.len(), 4); + // The square section bends alike in y and z: two pairs. + for (got, want, tol) in [(f[0], f1, 0.01), (f[1], f1, 0.01), (f[2], f2, 0.03), (f[3], f2, 0.03)] { + assert!((got - want).abs() / want < tol, "{name}: {got} vs {want}"); + } + assert!((r.total_mass - 20.0).abs() < 1e-9); + } +} + +#[test] +fn free_free_beam_has_six_rigid_body_modes() { + let g = beam(ElementKind::Hex20, [20, 1, 1], false, 8); + let r = solve_frequency(&g.model, Some("Step"), &FrequencyOptions::default(), &Limits::default()).unwrap(); + let f: Vec = r.modes.iter().map(|m| m.frequency).collect(); + let f1 = beam_frequency(4.730_040_74, 20.0, 1000.0, 1.0, 1.0); + eprintln!("free-free modes: {f:.6?} vs EB {f1:.6}"); + assert_eq!(f.len(), 8); + assert!(f[..6].iter().all(|x| *x < 1e-3 * f1), "{f:?}"); + assert!(f[6..].iter().all(|x| (x - f1).abs() / f1 < 0.015), "{f:?}"); + // The rigid translations carry the whole mass as effective mass. + for d in 0..3 { + let m: f64 = r.modes[..6].iter().map(|m| m.effective_mass[d]).sum(); + assert!((m - 20.0).abs() < 1e-6 * 20.0, "direction {d}: {m}"); + } + assert!(r.warnings.iter().any(|w| w.contains("rigid-body"))); +} + +#[test] +fn results_are_deterministic() { + let g = beam(ElementKind::Tet10, [6, 1, 1], true, 3); + let a = solve_frequency(&g.model, Some("Step"), &FrequencyOptions::default(), &Limits::default()).unwrap(); + let b = solve_frequency(&g.model, Some("Step"), &FrequencyOptions::default(), &Limits::default()).unwrap(); + assert_eq!(a, b); +} + +#[test] +fn errors_name_the_problem() { + // No boundary condition in a static step: singular, naming a node. + let mut g = grid([2, 1, 1], ElementKind::Hex8, |q| q); + g.define(iso(1000.0, 0.3), None, full(), Procedure::Static {}); + g.face_set("end", |q| near(q[0], 1.0)); + g.load("p", Load::Pressure { set: "end".into(), magnitude: 1.0 }); + let e = solve_static(&g.model, Some("Step"), &Limits::default()).unwrap_err().to_string(); + assert!(e.contains("singular at node"), "{e}"); + // Caps. + g.fixed("root", |q| near(q[0], 0.0)); + let tight = Limits { max_dofs: 10, ..Limits::default() }; + assert!(solve_static(&g.model, Some("Step"), &tight).unwrap_err().to_string().contains("limit")); + let tight = Limits { max_factor_entries: 5, ..Limits::default() }; + assert!(solve_static(&g.model, Some("Step"), &tight).unwrap_err().to_string().contains("limit")); + assert!(solve_static(&g.model, Some("Nope"), &Limits::default()).is_err()); + assert!(solve_frequency(&g.model, Some("Step"), &FrequencyOptions::default(), &Limits::default()).is_err(), "a static step"); + assert!(solve_static(&g.model, Some("Step"), &Limits::default()).is_ok()); +} diff --git a/docs/adr/0001-libraries.md b/docs/adr/0001-libraries.md index 1dd8be9..445af4b 100644 --- a/docs/adr/0001-libraries.md +++ b/docs/adr/0001-libraries.md @@ -26,8 +26,19 @@ predicates that `robust` lacks). `robust` (georust, MIT OR Apache-2.0) and `num- `num-rational` / `num-traits` (MIT OR Apache-2.0) are dev-dependencies only: differential test oracles, never linked into ModelCraft. -**Linear algebra:** to be decided in the solver task, among permissively licensed pure-Rust -crates (for example `faer`, MIT). Record the choice here. +**Linear algebra: our own.** *Decided 2026-10-13 (implicit solver task):* `modelcraft-solver` +has its own sparse symmetric storage, nested-dissection ordering (George 1973, by coordinate +bisection), up-looking sparse Cholesky (Davis 2006, elimination trees per Liu 1990) and +shift-invert block Lanczos with full reorthogonalization (Ericsson–Ruhe 1980, Grimes–Lewis–Simon +1994) on a Householder + implicit-QR dense eigensolver (Golub–Van Loan). `faer` (MIT) was the +alternative and is faster on large problems, but it selects SIMD kernels by CPU feature at run +time (AVX2, AVX-512, NEON), which changes floating-point reduction order between machines and +breaks the bit-for-bit determinism rule; it is also a large dependency tree for the browser +build. Our code is sequential with a fixed operation order, `unsafe`-free and wasm-safe, and its +factor size is known (and capped) before any value is computed. Cost if this is wrong: large +models (beyond a few hundred thousand equations) factor slower than they could; the remedy is a +supernodal factorization of our own, or `faer` behind the same `Factor` interface with a stated +tolerance instead of bit-identity, as the GPU exception already does. **Allowed:** MIT, BSD-2/3, Apache-2.0, ISC, Zlib, Unicode, MPL-2.0 (as an unmodified dependency), CC0 and public domain. diff --git a/docs/architecture.md b/docs/architecture.md index 693bf77..1a9e506 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -21,18 +21,20 @@ crates in strictly lower layers. | L0 | `modelcraft-model` | Solver-neutral FE model: nodes and elements by label, their geometry tags, sets, materials, sections, couplings, BCs, loads and steps, the model check and the resolved view the writers read. Knows no file format and no geometry type. | started | | L1 | `modelcraft-geometry` | The only door to SolidCraft: opens `.solidcraft` designs (rebuilt from their timeline, bodies with persistent face names), B-rep bodies, measures, tessellation into a face-tagged `SurfaceMesh`; face-tagged 3MF files read into the same type without the kernel (docs/formats/tagged-mesh.md). No SolidCraft type in its API. | started | | L1 | `modelcraft-mesh` | Surface preparation (weld, watertightness, feature edges), robust predicates, Delaunay tetrahedralization and refinement, Voronoi-based optimization, hex meshing, quality metrics. | started | -| L1 | `modelcraft-solver` | Element formulations, assembly, sparse linear static and modal solution. | planned | +| L1 | `modelcraft-fem` | Element kernels shared by every structural solver: the element table (tet4, tet10, hex8 full and incompatible modes, hex20 full and reduced), shape functions, quadrature, element stiffness, consistent and lumped mass, body and pressure loads, stresses and their extrapolation, isotropic and orthotropic elasticity. Pure functions, no solver loop. | started | | L1 | `modelcraft-decks` | Abaqus input, Nastran bulk data and LS-DYNA keyword readers and writers, to and from `model`. Text out, no file system; `deck.write` saves it. | started | | L1 | `modelcraft-render` | Display data without a GPU dependency: the visible faces, element edges and nodes of a model, the orbit camera, a CPU rasterizer for headless snapshots; contour maps later. | started | -| L2 | `modelcraft-doc` | The document: geometry, mesh, model, results, history. | planned | -| L3 | `modelcraft-io` | `.modelcraft` files (docs/formats/modelcraft.md): text model definition, binary mesh, geometry references with content hashes, the journal; atomic saves; autosave storage for crash recovery. Bytes in, bytes out; the file system only in `fs` and `recovery`. Result files later. | started | -| L4 | `modelcraft-engine` | Commands (`CommandSpec`), `Session::execute` with rollback and the panic guard, scripts. | started | -| L5 | `modelcraft-ui-egui`, `modelcraft-mcp` | Desktop UI (egui + wgpu viewport, model tree, commands, status bar) and control channel requests (docs/control-protocol.md); MCP server. | ui started | +| L2 | `modelcraft-solver` | Implicit solver: assembly from the resolved model, boundary conditions, sparse direct solution, linear static and modal analysis, stresses and reactions. The explicit (`modelcraft-explicit`) and lattice Boltzmann (`modelcraft-lbm`) solvers join this layer. | started | +| L3 | `modelcraft-fsi` | Fluid–structure coupling between the L2 solvers (later). | planned | +| L4 | `modelcraft-doc` | The document: geometry, mesh, model, results, history. | planned | +| L5 | `modelcraft-io` | `.modelcraft` files (docs/formats/modelcraft.md): text model definition, binary mesh, geometry references with content hashes, the journal; atomic saves; autosave storage for crash recovery. Bytes in, bytes out; the file system only in `fs` and `recovery`. Result files later. | started | +| L6 | `modelcraft-engine` | Commands (`CommandSpec`), `Session::execute` with rollback and the panic guard, scripts. | started | +| L7 | `modelcraft-ui-egui`, `modelcraft-mcp` | Desktop UI (egui + wgpu viewport, model tree, commands, status bar) and control channel requests (docs/control-protocol.md); MCP server. | ui started | | testkit | `modelcraft-testkit` | Test tooling, a dev-dependency only: the meshing test corpus (`tests/corpus/*/recipe.json`, analytic shapes) and its runner (`cargo xtask corpus`). | started | | apps | `modelcraft`, `modelcraft-cli`, `modelcraft-web` | Desktop, command line, browser. | desktop, cli started | Rules beyond the layer order: -- Nothing below L5 uses egui, eframe, winit, wgpu or rfd. +- Nothing below L7 uses egui, eframe, winit, wgpu or rfd. - Only `modelcraft-geometry` uses `solidcraft-*`, pinned by rev in the workspace `Cargo.toml`. SolidCraft's `[patch.crates-io]` (its fixed truck crates) does not reach dependents, so the workspace repeats it from the same rev; bump both together. The same goes for SolidCraft's @@ -120,6 +122,39 @@ element names, cards and keywords): - **Files.** A model deserializes by replaying its contents through the same checked methods the commands use, so a file cannot hold what commands would refuse. +## Implicit solver + +`modelcraft-solver` solves a model's linear static and frequency steps (`solve.static`, +`solve.frequency`; `solve.probe` reads the stored results). Methods and papers are cited in the +module docs; ADR 0001 records why the linear algebra is our own. + +1. **System.** `Model::resolve()`; each solid section's elements take their material's + elasticity matrix and the section's formulation (`modelcraft_fem::element::ELEMENT_TYPES`, one + row per element kind). Couplings, reference points, moments, shells and beams return "not + supported yet". The step's boundary conditions prescribe translations (rotations do not exist + on solid nodes); its loads become a load vector: forces and tractions from + `Resolved::nodal_loads`, pressure as consistent face loads, gravity and body forces through the + element shape functions. +2. **Assembly** on the node graph (3 × 3 blocks, elements in label order, so sums are + reproducible bit for bit). +3. **Ordering and elimination.** Nested dissection of the nodes by coordinate bisection; free + degrees of freedom numbered in that order; prescribed values moved to the right-hand side (no + penalty). +4. **Static solution.** Up-looking sparse Cholesky (the factor size is computed and capped before + any value), reactions K u − f, integration-point stresses extrapolated to the nodes and + averaged, von Mises. A zero pivot names the node and direction that can move freely. +5. **Frequency.** Shift-invert block Lanczos (block 6, full reorthogonalization) on + K − σM, σ = 0 for a supported structure and slightly negative for a free one, so rigid-body + modes come out at zero; mass-normalized shapes, participation factors and effective masses. +6. **Results** stay in the session per step, stamped with the session revision they were solved + at; `Session::result` hands out only those of the current revision, so an edit makes them + stale and undoing back to the solved state makes them current again. + +Verification (`crates/solver/tests/verification.rs`, analytic only): displacement and traction +patch tests for every element and formulation (exact to round-off), cantilever tip deflection +against Timoshenko with convergence rates of the linear elements, Lamé's cylinder, Kirsch's +plate with a hole (Kt → 3), cantilever and free–free beam frequencies with six rigid-body modes. + ## Commands Every user-visible action is a `CommandSpec` in `crates/engine/src/cmd/.rs` with a JSON diff --git a/docs/parity.md b/docs/parity.md index aac208f..03eff04 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -60,8 +60,8 @@ products are run, no documentation text is copied). Status: **done**, **started* | Decks | Nastran bulk data writer (`deck.write`, SOL 101 or 103 per deck; reports what Nastran cannot express) | M3 | | done | | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | done (implicit static and eigenvalue; beam sections planned) | | Decks | Readers keeping unknown keywords | M9 | | planned | -| Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | planned | -| Results | Deformed shape, contours, mode shapes, probes, reaction sums | M4 | | planned | +| Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, hex8, hex20) | +| Results | Deformed shape, contours, mode shapes, probes, reaction sums | M4 | | started (results in the session, `solve.probe`, reaction sums; display planned) | | Results | XY and path plots, cuts, display groups, free-body forces, animation | M9 | | planned | | Results | Readers for public result formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`) | M9 | | planned | | Results | Mesh convergence studies | M9 | | planned | diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index e07a719..bc58f97 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -953,3 +953,30 @@ the better pairing — this is why FSI should support both structure solvers. tagged as "benchmark", kept separate from the analytic suite). 8. **Quad shell meshing is missing from the roadmap.** X5 needs quad-dominant surface meshes; without them shells run on tri3 (BLT has a triangular variant, less accurate). Flag for the roadmap. + +## Part G: Implicit solver (M4), as implemented + +The papers and books `modelcraft-fem` and `modelcraft-solver` are written from (cited again in +the module docs): + +- Elements: O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, *The Finite Element Method: Its Basis + and Fundamentals*, 7th ed., 2013; T. J. R. Hughes, *The Finite Element Method*, Dover 2000. +- Incompatible modes: Wilson, Taylor, Doherty, Ghaboussi 1973; R. L. Taylor, P. J. Beresford, + E. L. Wilson, IJNME 10(6):1211–1219, 1976 (the det J₀/det J correction that passes the patch + test). +- Tetrahedron quadrature: A. H. Stroud, *Approximate Calculation of Multiple Integrals*, 1971 + (conical product rules). +- Mass lumping: E. Hinton, T. Rock, O. C. Zienkiewicz, Earthquake Eng. Struct. Dyn. + 4(3):245–249, 1976. +- Stress extrapolation: E. Hinton, J. S. Campbell, IJNME 8(3):461–480, 1974. +- Ordering: A. George, SIAM J. Numer. Anal. 10(2):345–363, 1973; H. D. Simon, Computing Systems + in Engineering 2:135–148, 1991. +- Sparse Cholesky: T. A. Davis, *Direct Methods for Sparse Linear Systems*, SIAM 2006; J. W. H. + Liu, SIAM J. Matrix Anal. Appl. 11(1):134–172, 1990. +- Eigenvalues: T. Ericsson, A. Ruhe, Math. Comp. 35(152):1251–1268, 1980; R. G. Grimes, J. G. + Lewis, H. D. Simon, SIAM J. Matrix Anal. Appl. 15(1):228–272, 1994; B. N. Parlett, *The + Symmetric Eigenvalue Problem*, SIAM 1998; G. H. Golub, C. F. Van Loan, *Matrix Computations*, + 4th ed., 2013. +- Verification: R. H. MacNeal, R. L. Harder, *A proposed standard set of problems to test finite + element accuracy*, Finite Elem. Anal. Des. 1(1):3–20, 1985 (patch tests); Lamé, Kirsch, + Timoshenko and Euler–Bernoulli closed forms. diff --git a/xtask/src/layers.rs b/xtask/src/layers.rs index 80b40ba..94927e3 100644 --- a/xtask/src/layers.rs +++ b/xtask/src/layers.rs @@ -14,8 +14,8 @@ pub enum Class { Layer(u8), /// Layer 0, and additionally may depend on no workspace crate at all. Standalone, - /// Test tooling: may depend on anything up to L5 (it sits at L6 for rule - /// purposes); other crates may use it only as a dev-dependency. + /// Test tooling: may depend on anything up to L7 (it sits at [`TESTKIT_LAYER`] + /// for rule purposes); other crates may use it only as a dev-dependency. Testkit, /// Binaries and build tooling: exempt from the rules. Exempt, @@ -26,7 +26,7 @@ impl Class { match self { Class::Layer(l) => Some(l), Class::Standalone => Some(0), - Class::Testkit => Some(6), + Class::Testkit => Some(TESTKIT_LAYER), Class::Exempt => None, } } @@ -35,17 +35,23 @@ impl Class { /// The layering table. Names are package names without the `modelcraft-` /// prefix. pub const TABLE: &[(&str, Class)] = &[ + // L0: the solver-neutral model. ("model", Class::Layer(0)), + // L1: geometry and math kernels. ("geometry", Class::Layer(1)), ("mesh", Class::Layer(1)), - ("solver", Class::Layer(1)), + ("fem", Class::Layer(1)), ("decks", Class::Layer(1)), ("render", Class::Layer(1)), - ("doc", Class::Layer(2)), - ("io", Class::Layer(3)), - ("engine", Class::Layer(4)), - ("ui-egui", Class::Layer(5)), - ("mcp", Class::Layer(5)), + // L2: solvers (the explicit and lattice Boltzmann solvers join here). + ("solver", Class::Layer(2)), + // L3: coupling between solvers (`fsi`, later). + // L4 and up: the document, files, commands, UI. + ("doc", Class::Layer(4)), + ("io", Class::Layer(5)), + ("engine", Class::Layer(6)), + ("ui-egui", Class::Layer(7)), + ("mcp", Class::Layer(7)), ("testkit", Class::Testkit), // apps and tooling ("modelcraft", Class::Exempt), @@ -73,7 +79,13 @@ pub const UI_CRATES: &[&str] = &["egui", "eframe", "egui-wgpu", "wgpu", "winit", pub const SOLIDCRAFT_GATEWAY: &str = "geometry"; /// First layer allowed to use UI crates. -pub const UI_MIN_LAYER: u8 = 5; +pub const UI_MIN_LAYER: u8 = 7; + +/// The layer the testkit sits at for rule purposes: above every layered crate. +pub const TESTKIT_LAYER: u8 = UI_MIN_LAYER + 1; + +/// The layer reported for apps and tooling when a layered crate depends on one. +const APPS_LAYER: u8 = TESTKIT_LAYER + 1; pub fn short_name(pkg: &str) -> &str { pkg.strip_prefix("modelcraft-").unwrap_or(pkg) @@ -118,7 +130,7 @@ pub enum Violation { Upward { krate: String, dep: String, from: u8, to: u8, kind: DepKind }, StandaloneHasWorkspaceDep { krate: String, dep: String }, TestkitAsNormalDep { krate: String }, - UiBelowL5 { krate: String, dep: String, layer: u8 }, + UiBelowUiLayer { krate: String, dep: String, layer: u8 }, SolidCraftOutsideGeometry { krate: String, dep: String }, } @@ -137,8 +149,8 @@ impl std::fmt::Display for Violation { Violation::TestkitAsNormalDep { krate } => { write!(f, "{krate}: modelcraft-testkit may only be a dev-dependency") } - Violation::UiBelowL5 { krate, dep, layer } => { - write!(f, "{krate} (L{layer}) depends on UI crate `{dep}`; UI toolkits are only allowed in L5+") + Violation::UiBelowUiLayer { krate, dep, layer } => { + write!(f, "{krate} (L{layer}) depends on UI crate `{dep}`; UI toolkits are only allowed in L{UI_MIN_LAYER}+") } Violation::SolidCraftOutsideGeometry { krate, dep } => { write!(f, "{krate} depends on `{dep}`; only modelcraft-{SOLIDCRAFT_GATEWAY} may use SolidCraft") @@ -183,7 +195,7 @@ pub fn check(crates: &[Crate]) -> Vec { krate: c.name.clone(), dep: d.name.clone(), from: layer, - to: if to == u8::MAX { 7 } else { to }, + to: to.min(APPS_LAYER), kind: d.kind, }); } @@ -192,7 +204,7 @@ pub fn check(crates: &[Crate]) -> Vec { } else if d.name.starts_with("solidcraft-") && short_name(&c.name) != SOLIDCRAFT_GATEWAY { out.push(Violation::SolidCraftOutsideGeometry { krate: c.name.clone(), dep: d.name.clone() }); } else if layer < UI_MIN_LAYER && is_ui_crate(&d.name) { - out.push(Violation::UiBelowL5 { krate: c.name.clone(), dep: d.name.clone(), layer }); + out.push(Violation::UiBelowUiLayer { krate: c.name.clone(), dep: d.name.clone(), layer }); } } } @@ -259,12 +271,14 @@ mod tests { #[test] fn upward_dependency_flagged() { let v = check(&[c("modelcraft-doc", &[("modelcraft-engine", Normal, true)])]); - assert!(matches!(v[..], [Violation::Upward { from: 2, to: 4, .. }])); + assert!(matches!(v[..], [Violation::Upward { from: 4, to: 6, .. }])); + let v = check(&[c("modelcraft-fem", &[("modelcraft-solver", Normal, true)])]); + assert!(matches!(v[..], [Violation::Upward { from: 1, to: 2, .. }])); } #[test] fn sideways_dependency_flagged() { - let v = check(&[c("modelcraft-mesh", &[("modelcraft-solver", Normal, true)])]); + let v = check(&[c("modelcraft-mesh", &[("modelcraft-fem", Normal, true)])]); assert!(matches!(v[..], [Violation::Upward { from: 1, to: 1, .. }])); } @@ -280,10 +294,26 @@ mod tests { } #[test] - fn ui_crates_below_l6_flagged() { + fn solvers_sit_above_the_kernels() { + let g = [ + c("modelcraft-fem", &[("modelcraft-model", Normal, true)]), + c("modelcraft-solver", &[("modelcraft-fem", Normal, true), ("modelcraft-mesh", Normal, true)]), + c("modelcraft-engine", &[("modelcraft-solver", Normal, true)]), + ]; + assert!(check(&g).is_empty(), "{:?}", check(&g)); + let v = check(&[c("modelcraft-solver", &[("modelcraft-engine", Normal, true)])]); + assert!(matches!(v[..], [Violation::Upward { from: 2, to: 6, .. }])); + let v = check(&[c("modelcraft-engine", &[("modelcraft-cli", Normal, true)])]); + assert!(matches!(v[..], [Violation::Upward { from: 6, to: APPS_LAYER, .. }])); + } + + #[test] + fn ui_crates_below_the_ui_layer_flagged() { for dep in ["egui", "eframe", "winit", "egui_kittest", "rfd", "bevy_ecs", "bevy"] { let v = check(&[c("modelcraft-engine", &[(dep, Normal, false)])]); - assert!(matches!(v[..], [Violation::UiBelowL5 { layer: 4, .. }]), "{dep}"); + assert!(matches!(v[..], [Violation::UiBelowUiLayer { layer: 6, .. }]), "{dep}"); + let v = check(&[c("modelcraft-solver", &[(dep, Normal, false)])]); + assert!(matches!(v[..], [Violation::UiBelowUiLayer { layer: 2, .. }]), "{dep}"); } assert!(check(&[c("modelcraft-engine", &[("egui_extras_not", Normal, false)])]).is_empty()); assert!(check(&[c("modelcraft-ui-egui", &[("winit", Normal, false)])]).is_empty());