diff --git a/Cargo.lock b/Cargo.lock index 2de351c..3d77841 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1880,6 +1880,7 @@ version = "0.1.0" dependencies = [ "log", "modelcraft-decks", + "modelcraft-explicit", "modelcraft-geometry", "modelcraft-io", "modelcraft-mesh", @@ -1891,6 +1892,16 @@ dependencies = [ "zip", ] +[[package]] +name = "modelcraft-explicit" +version = "0.1.0" +dependencies = [ + "modelcraft-fem", + "modelcraft-model", + "serde", + "thiserror 2.0.21", +] + [[package]] name = "modelcraft-fem" version = "0.1.0" diff --git a/Cargo.toml b/Cargo.toml index 93f2e69..cd9bf5c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -19,6 +19,7 @@ modelcraft-decks = { path = "crates/decks" } modelcraft-io = { path = "crates/io" } modelcraft-fem = { path = "crates/fem" } modelcraft-solver = { path = "crates/solver" } +modelcraft-explicit = { path = "crates/explicit" } modelcraft-engine = { path = "crates/engine" } modelcraft-testkit = { path = "crates/testkit" } modelcraft-render = { path = "crates/render" } diff --git a/ROADMAP.md b/ROADMAP.md index b232473..ada507b 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -28,7 +28,7 @@ hex meshes today and it gives the best elements where stresses are read. | 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 | | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; heat transfer | planned | | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned | -| M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting; large-deformation materials (objective rates, J2, Johnson–Cook, damage, erosion); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | planned | +| M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `modelcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (objective rates, J2, Johnson–Cook, damage, erosion); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | started | | M11 Equations of state and shells | Polynomial, ideal-gas, Mie–Grüneisen, JWL with programmed burn, tabulated; energy coupling; shock verification (Rankine–Hugoniot, Sod, Noh, Sedov); Belytschko–Lin–Tsay shells with self-contact on quad meshes | planned | | M12 CFD: lattice Boltzmann | D3Q19/D3Q27 on grids voxelised from the face-tagged surface; TRT, regularised and cumulant collision; interpolated walls, inlets and outlets; LES; octree refinement shared with the hex mesher; wgpu backend (deterministic per device, checked against the CPU) | planned | | M13 Fluid–structure interaction | Immersed-boundary coupling of the lattice Boltzmann fluid with the explicit or implicit structure; interface force transfer; sub-cycling; Aitken then interface quasi-Newton coupling; free-surface flow | planned | diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index f38fc8a..835cf19 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -216,6 +216,9 @@ fn check_expressible(model: &Model) -> Result<()> { )); } } + for step in model.steps().iter().filter(|s| matches!(s.procedure, Procedure::ExplicitDynamic(_))) { + bad.insert(format!("step `{}`: explicit dynamic steps (not written to Abaqus input yet)", step.name)); + } for (name, c) in model.couplings() { if c.kind == CouplingKind::Distributing && !(1..=3).all(|d| c.dofs.contains(d)) { bad.insert(format!("coupling `{name}`: an Abaqus distributing coupling always couples translations 1-3")); @@ -787,6 +790,9 @@ impl Deck { self.comment("loads have no effect on natural frequencies and are not written"); } } + Procedure::ExplicitDynamic(_) => { + return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: explicit dynamic steps", step.name) }); + } } self.output(step); self.line("*END STEP"); diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 66e8532..1f2c8ed 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -235,6 +235,9 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> { } let mut gravity_sets: BTreeSet> = BTreeSet::new(); for step in model.steps() { + if matches!(step.procedure, Procedure::ExplicitDynamic(_)) { + refused.insert(format!("step `{}`: explicit dynamic steps (this writer writes implicit analyses so far)", step.name)); + } let is_static = matches!(step.procedure, Procedure::Static {}); for b in &step.boundary_conditions { let Some(bc) = model.boundary_conditions().get(b) else { continue }; @@ -584,6 +587,9 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste d.card([Field::Blank, Field::Blank, Field::Blank, Field::Blank, Field::from(1)]); } } + Procedure::ExplicitDynamic(_) => { + return Err(Error::Unsupported { format: FORMAT, what: format!("step `{}`: explicit dynamic steps", step.name) }); + } } d.line("*CONTROL_TERMINATION"); d.header(&["endtim"]); diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index ecd19c9..7fe8226 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -234,6 +234,9 @@ struct Deck<'m> { /// The solution sequence of a list of steps, or why there is none. fn solution(steps: &[&Step]) -> std::result::Result { + if let Some(s) = steps.iter().find(|s| matches!(s.procedure, Procedure::ExplicitDynamic(_))) { + return Err(format!("step `{}`: explicit dynamic steps (Nastran's linear solutions have none)", s.name)); + } let statics = steps.iter().find(|s| matches!(s.procedure, Procedure::Static {})); let modes = steps.iter().find(|s| matches!(s.procedure, Procedure::Frequency { .. })); match (statics, modes) { @@ -581,8 +584,8 @@ impl<'m> Deck<'m> { comment(&mut out, "eigenvalue extraction"); } for sub in &self.subcases { - if let (Some(sid), Procedure::Frequency { modes, min_hz, max_hz }) = (sub.method, sub.step.procedure) { - Card::new("EIGRL").int(sid).opt(min_hz).opt(max_hz).int(modes).write(&mut out)?; + if let (Some(sid), Procedure::Frequency { modes, min_hz, max_hz }) = (sub.method, &sub.step.procedure) { + Card::new("EIGRL").int(sid).opt(*min_hz).opt(*max_hz).int(*modes).write(&mut out)?; } } out.push_str("ENDDATA\n"); diff --git a/crates/decks/tests/lsdyna.rs b/crates/decks/tests/lsdyna.rs index 8c75cee..ae82486 100644 --- a/crates/decks/tests/lsdyna.rs +++ b/crates/decks/tests/lsdyna.rs @@ -674,3 +674,15 @@ fn hostile_numbers_stay_within_their_fields() { let deck = read(&w.text); assert_eq!(deck.nodes[&1], [1e9, -1e9, 1e-300]); } + +/// Explicit dynamic steps are not written by any deck format yet: each says so by name. +#[test] +fn explicit_steps_are_refused_by_every_format() { + let mut m = cantilever(); + m.add_step(step("Impact", Procedure::ExplicitDynamic(ExplicitDynamic::new(1e-3)), &["clamp"], &[]), None).unwrap(); + for f in modelcraft_decks::FORMATS { + let e = f.write(&m).unwrap_err(); + assert!(matches!(e, Error::Unsupported { .. }), "{}: {e}", f.name); + assert!(e.to_string().contains("step `Impact`: explicit dynamic"), "{}: {e}", f.name); + } +} diff --git a/crates/engine/Cargo.toml b/crates/engine/Cargo.toml index 2c4bfb5..f3d0478 100644 --- a/crates/engine/Cargo.toml +++ b/crates/engine/Cargo.toml @@ -15,6 +15,7 @@ modelcraft-mesh = { workspace = true } modelcraft-io = { workspace = true } modelcraft-decks = { workspace = true } modelcraft-solver = { workspace = true } +modelcraft-explicit = { workspace = true } serde = { workspace = true } serde_json = { workspace = true } thiserror = { workspace = true } diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index 5ef2803..ebafb67 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -292,7 +292,13 @@ pub static CATALOGUE: &[Capability] = &[ soon("deck.read", "Read Deck", Solve, "Decks", "M9", "Read a deck, keeping keywords ModelCraft does not know"), ready("solve.static", "Static", Solve, "Implicit", "Linear static solve of a step with ModelCraft's own solver"), ready("solve.frequency", "Frequency", Solve, "Implicit", "Natural frequencies and mode shapes of a frequency step"), - soon("solve.explicit", "Explicit", Solve, "Explicit", "M10", "Explicit dynamics: central differences, contact, erosion"), + ready( + "solve.explicit", + "Explicit", + Solve, + "Explicit", + "Transient dynamics of an explicit step by central differences (contact and erosion coming in M10)", + ), soon("solve.blast", "Blast & Impact", Solve, "Explicit", "M14", "EOS-driven fluids with fluid–structure contact (ALE or SPH)"), soon("solve.flow", "Lattice Boltzmann", Solve, "Fluids", "M12", "Lattice Boltzmann flow on the CPU or the GPU"), soon("solve.fsi", "Fluid–Structure", Solve, "Fluids", "M13", "Coupled fluid and structure solve"), @@ -301,7 +307,14 @@ pub static CATALOGUE: &[Capability] = &[ soon("results.contour", "Contours", Results, "Display", "M4", "Contour plots of stress, strain and displacement"), soon("results.modes", "Mode Shapes", Results, "Display", "M4", "Mode shapes of a frequency step"), soon("results.animate", "Animate", Results, "Display", "M9", "Animate modes and time histories"), - ready("solve.probe", "Probe", Results, "Query", "Displacement, reaction and stress at nodes, or a mode shape"), + ready_more( + "solve.probe", + "Probe", + Results, + "Query", + "Displacement, reaction and stress at nodes, a mode shape, or an explicit field frame", + &["solve.history"], + ), soon("results.reactions", "Reactions", Results, "Query", "M4", "Reaction force and moment sums"), soon("results.free_body", "Free Body", Results, "Query", "M9", "Free-body forces across a cut"), soon("results.xy", "XY Plot", Results, "Plots", "M9", "Plot a quantity against time or another"), diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index e4c4310..ed2b3c6 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -159,7 +159,7 @@ definition_commands!( pub const STEP_CREATE: CommandSpec = CommandSpec { id: "model.step.create", label: "New Step", - params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, + params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal"}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, kind: CommandKind::Model, run: step_create, }; diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index b32489e..41c0e08 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -71,6 +71,8 @@ pub static ALL: &[CommandSpec] = &[ deck::WRITE, solve::STATIC, solve::FREQUENCY, + solve::EXPLICIT, + solve::HISTORY, solve::PROBE, ]; diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 9d9d7b8..2f3c33a 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -1,11 +1,12 @@ -//! `solve.*`: run the model's steps with the implicit solver, keep the results in the session, -//! and read them back. +//! `solve.*`: run the model's steps with the implicit or the explicit solver, keep the results +//! in the session, and read them back. -use modelcraft_solver::{FrequencyOptions, Limits, Location, MassKind, StepResult, solve_frequency, solve_static}; +use modelcraft_explicit::{ExplicitResult, solve_explicit}; +use modelcraft_solver::{FrequencyOptions, Limits, Location, MassKind, StepResult as ImplicitResult, solve_frequency, solve_static}; use serde_json::{Value, json}; -use super::{labels, opt}; -use crate::{CommandKind, CommandSpec, EngineError, Result, Session}; +use super::{labels, opt, req}; +use crate::{CommandKind, CommandSpec, EngineError, Result, Session, StepResult}; pub const STATIC: CommandSpec = CommandSpec { id: "solve.static", @@ -23,10 +24,26 @@ pub const FREQUENCY: CommandSpec = CommandSpec { run: run_frequency, }; +pub const EXPLICIT: CommandSpec = CommandSpec { + id: "solve.explicit", + label: "Solve Explicit Dynamic Step", + params: r#"{"step": "Impact", "max_dofs": 6000000, "max_increments": 1000000000, "max_output_values": 50000000} — transient response of an explicit dynamic step (the only one when `step` is left out) by central differences; keeps its histories and field frames in the session; the caps may only be lowered (the step's own increment cap applies too)"#, + kind: CommandKind::Model, + run: run_explicit, +}; + +pub const HISTORY: CommandSpec = CommandSpec { + id: "solve.history", + label: "Result Histories", + params: r#"{"step": "Impact", "every": 1} — the time histories of a solved explicit step: time, increment, energies, momentum, constraint force and the history nodes' displacement and velocity; `every` keeps one sample in so many (by default enough to return at most 10000)"#, + kind: CommandKind::Model, + run: history, +}; + 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); a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments)"#, + params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1, "frame": 0} — the stored result at nodes: displacement, reaction, stress and von Mises (static), the shape of mode `mode` (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default); a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments)"#, kind: CommandKind::Model, run: probe, }; @@ -34,14 +51,18 @@ pub const PROBE: CommandSpec = CommandSpec { /// Most nodes one probe lists. const MAX_PROBE: usize = 10_000; +/// Most history samples one `solve.history` returns. +const MAX_HISTORY_SAMPLES: usize = 10_000; + /// An optional cap that may only lower the default. -fn cap(p: &Value, key: &str, default: usize) -> Result { +fn cap + TryFrom + std::fmt::Display>(p: &Value, key: &str, default: T) -> Result { + let most: u64 = default.try_into().unwrap_or(u64::MAX); 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)) + .filter(|n| (1..=most).contains(n)) + .and_then(|n| T::try_from(n).ok()) .ok_or_else(|| EngineError::Params(format!("`{key}` must be a whole number from 1 to {default}"))), } } @@ -83,7 +104,7 @@ fn run_static(s: &mut Session, p: &Value) -> Result { })).collect::>(), "warnings": r.warnings, }); - s.store_result(StepResult::Static(r)); + s.store_result(StepResult::Implicit(ImplicitResult::Static(r))); Ok(summary) } @@ -107,21 +128,99 @@ fn run_frequency(s: &mut Session, p: &Value) -> Result { })).collect::>(), "warnings": r.warnings, }); - s.store_result(StepResult::Frequency(r)); + s.store_result(StepResult::Implicit(ImplicitResult::Frequency(r))); Ok(summary) } +fn run_explicit(s: &mut Session, p: &Value) -> Result { + use modelcraft_explicit::Limits as ExplicitLimits; + let limits = ExplicitLimits { + max_dofs: cap(p, "max_dofs", ExplicitLimits::MAX_DOFS)?, + max_increments: cap(p, "max_increments", ExplicitLimits::MAX_INCREMENTS)?, + max_output_values: cap(p, "max_output_values", ExplicitLimits::MAX_OUTPUT_VALUES)?, + }; + let step = opt::(p, "step")?; + let r = solve_explicit(&s.model, step.as_deref(), &limits)?; + let summary = explicit_summary(&r); + s.store_result(StepResult::Explicit(r)); + Ok(summary) +} + +fn explicit_summary(r: &ExplicitResult) -> Value { + json!({ + "step": r.step, + "dofs": 3 * r.nodes.len(), + "elements": r.elements.len(), + "increments": r.increments, + "time": r.time, + "completed": r.completed, + "initial_increment": r.initial_increment, + "controlling_elements": r.controlling, + "total_mass": r.total_mass, + "energy": r.energy, + "energy_balance": r.energy.relative_error(r.initial_kinetic), + "frames": r.frames.len(), + "history_samples": r.history.time.len(), + "warnings": r.warnings, + }) +} + +fn history(s: &mut Session, p: &Value) -> Result { + let step: String = req(p, "step")?; + let result = s.result(&step).ok_or_else(|| EngineError::Params(format!("step `{step}` has no result; solve it first")))?; + let StepResult::Explicit(r) = result.as_ref() else { + return Err(EngineError::Params(format!("step `{step}` is not an explicit dynamic step; it has no histories"))); + }; + let h = &r.history; + let samples = h.time.len(); + let every = match opt::(p, "every")? { + Some(0) => return Err(EngineError::Params("`every` must be at least 1".into())), + Some(k) => usize::try_from(k).unwrap_or(usize::MAX), + None => samples.div_ceil(MAX_HISTORY_SAMPLES).max(1), + }; + if samples.div_ceil(every) > MAX_HISTORY_SAMPLES { + return Err(EngineError::Params(format!("{samples} samples in steps of {every} are more than {MAX_HISTORY_SAMPLES}; raise `every`"))); + } + let pick = |v: &dyn Fn(usize) -> Value| -> Vec { (0..samples).step_by(every).map(v).collect() }; + let nodes: Vec = h + .nodes + .iter() + .enumerate() + .map(|(k, label)| { + json!({ + "node": label, + "displacement": pick(&|i| json!(h.displacement.get(i).and_then(|d| d.get(k)))), + "velocity": pick(&|i| json!(h.velocity.get(i).and_then(|d| d.get(k)))), + }) + }) + .collect(); + Ok(json!({ + "step": step, + "every": every, + "time": pick(&|i| json!(h.time.get(i))), + "increment": pick(&|i| json!(h.increment.get(i))), + "energies": pick(&|i| json!(h.energies.get(i))), + "momentum": pick(&|i| json!(h.momentum.get(i))), + "constraint_force": pick(&|i| json!(h.constraint_force.get(i))), + "nodes": nodes, + })) +} + 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 result = match result.as_ref() { + StepResult::Explicit(r) => return probe_frame(r, &step, &nodes, p), + StepResult::Implicit(r) => r, + }; let locate = |label: u32| { result .locate(label) .ok_or_else(|| EngineError::Params(format!("node {label} has no result (it belongs to no solved element or coupled reference point)"))) }; - let out: Vec = match result.as_ref() { - StepResult::Static(r) => nodes + let out: Vec = match result { + ImplicitResult::Static(r) => nodes .iter() .map(|&l| { Ok(match locate(l)? { @@ -139,7 +238,7 @@ fn probe(s: &mut Session, p: &Value) -> Result { }) }) .collect::>()?, - StepResult::Frequency(r) => { + ImplicitResult::Frequency(r) => { let mode = opt::(p, "mode")?.unwrap_or(1); let m = mode .checked_sub(1) @@ -161,6 +260,25 @@ fn probe(s: &mut Session, p: &Value) -> Result { Ok(json!({"step": step, "nodes": out})) } +/// Displacement and velocity at nodes in one field frame of an explicit result. +fn probe_frame(r: &ExplicitResult, step: &str, nodes: &[u32], p: &Value) -> Result { + let last = r.frames.len().saturating_sub(1); + let k = opt::(p, "frame")?.unwrap_or(last); + let f = r.frames.get(k).ok_or_else(|| EngineError::Params(format!("`frame` must be from 0 to {last}")))?; + if f.displacement.is_empty() { + return Err(EngineError::Params(format!("step `{step}` does not output displacements"))); + } + let out: Vec = nodes + .iter() + .map(|&l| { + let i = + r.nodes.binary_search(&l).map_err(|_| EngineError::Params(format!("node {l} has no result (it belongs to no solved element)")))?; + Ok(json!({"node": l, "displacement": f.displacement.get(i), "velocity": f.velocity.get(i)})) + }) + .collect::>()?; + Ok(json!({"step": step, "frame": k, "time": f.time, "nodes": out})) +} + #[cfg(test)] pub(crate) mod tests { use super::*; @@ -306,4 +424,211 @@ pub(crate) mod tests { assert!(s.execute("solve.probe", &p).is_err(), "{p}"); } } + + /// A cantilever L × 1 × 1 along x of `n[0]` × `n[1]` × `n[1]` hex8 elements, E = 1000, ν = 0, + /// ρ = 1, clamped at x = 0; one node set per cross-section (`x0`, `x1`, …) and the tip's + /// corner node `tip`; a frequency step `Modes` and the section `beam`. + fn explicit_cantilever(s: &mut Session, length: f64, n: [usize; 2]) { + use modelcraft_model::{ElementKind, Set}; + let [nx, ny] = n; + let label = |i: usize, j: usize, k: usize| (1 + i + (nx + 1) * (j + (ny + 1) * k)) as u32; + for k in 0..=ny { + for j in 0..=ny { + for i in 0..=nx { + s.model.set_node(label(i, j, k), [length * i as f64 / nx as f64, j as f64 / ny as f64, k as f64 / ny as f64]).unwrap(); + } + } + } + let corners = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]; + let mut e = 1; + for k in 0..ny { + for j in 0..ny { + for i in 0..nx { + s.model.set_element(e, ElementKind::Hex8, corners.iter().map(|c| label(i + c[0], j + c[1], k + c[2])).collect()).unwrap(); + e += 1; + } + } + } + for i in 0..=nx { + let nodes = (0..=ny).flat_map(|k| (0..=ny).map(move |j| label(i, j, k))).collect(); + s.model.add_set(&format!("x{i}"), Set { nodes, ..Set::default() }).unwrap(); + } + s.model.add_set("tip", Set { nodes: vec![label(nx, 0, 0)], ..Set::default() }).unwrap(); + s.model.add_set("all", Set { elements: (1..e).collect(), ..Set::default() }).unwrap(); + let script = json!([ + {"command": "model.material.create", "params": {"name": "Mat", "elastic": {"type": "isotropic", "young": 1000, "poisson": 0}, "density": 1}}, + {"command": "model.section.create", "params": {"name": "beam", "set": "all", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"incompatible_modes": true}}}, + {"command": "model.boundary_condition.create", "params": {"name": "clamp", "type": "fixed", "target": {"set": "x0"}}}, + {"command": "model.step.create", "params": {"name": "Modes", "procedure": {"type": "frequency", "modes": 1}, "boundary_conditions": ["clamp"]}}, + ]); + s.run_script(&script).unwrap(); + } + + /// Upward zero crossings of a sampled signal, linearly interpolated. + fn crossings(time: &[f64], x: &[f64]) -> Vec { + let mut out = Vec::new(); + for k in 1..x.len() { + let (a, b) = (x[k - 1], x[k]); + if a < 0.0 && b >= 0.0 { + out.push(time[k - 1] + (time[k] - time[k - 1]) * (-a) / (b - a)); + } + } + out + } + + /// Free vibration of a cantilever after an initial velocity in the shape of its static + /// deflection: the period of the explicit run (one-point hex8) against the first natural + /// frequency of the implicit modal solver (incompatible-mode hex8, lumped mass) on the same + /// mesh. The one-point element's mean strain misses the linear variation of the bending + /// strain through each element: with n elements through the depth it keeps 1 − 1/n² of the + /// bending energy, so without hourglass stiffness (viscous control) the period is longer by + /// (1 − 1/n²)^{−½} (3.3% for n = 4); the stiffness form restores part of it. + #[test] + fn explicit_cantilever_period_matches_the_modal_frequency() { + let (length, n) = (10.0, [40usize, 4usize]); + let mut s = Session::default(); + explicit_cantilever(&mut s, length, n); + let modes = s.execute("solve.frequency", &json!({"mass": "lumped"})).unwrap(); + let period = 1.0 / modes["modes"][0]["frequency"].as_f64().unwrap(); + // Euler–Bernoulli: f₁ = 1.875² / (2π L²) √(E I / (ρ A)), I/A = 1/12. + let euler = 2.0 * std::f64::consts::PI * length * length / 1.875104f64.powi(2) / (1000.0f64 / 12.0).sqrt(); + assert!((period / euler - 1.0).abs() < 0.01, "{period} vs {euler}"); + let section = json!({"name": "beam", "set": "all", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}); + s.execute("model.section.edit", §ion).unwrap(); + // v(x) ∝ x²(3L − x): the static deflection of a tip load, close to the first mode. + let velocities: Vec = (0..=n[0]) + .map(|i| { + let x = length * i as f64 / n[0] as f64; + json!({"set": format!("x{i}"), "velocity": [0, 0, 0.01 * x * x * (3.0 * length - x) / (2.0 * length.powi(3))]}) + }) + .collect(); + let one_point = period / (1.0 - 1.0 / (n[1] * n[1]) as f64).sqrt(); + let mut measured = Vec::new(); + for hourglass in ["viscous", "stiffness"] { + let step = json!({ + "type": "explicit_dynamic", "end_time": 4.2 * period, "history_interval": period / 400.0, "output_interval": period, + "history_nodes": "tip", "initial_velocities": velocities, "hourglass": {"type": hourglass} + }); + let name = format!("Swing-{hourglass}"); + s.execute("model.step.create", &json!({"name": name, "procedure": step, "boundary_conditions": ["clamp"]})).unwrap(); + let r = s.execute("solve.explicit", &json!({"step": name})).unwrap(); + let h = s.execute("solve.history", &json!({"step": name})).unwrap(); + let time: Vec = h["time"].as_array().unwrap().iter().map(|t| t.as_f64().unwrap()).collect(); + let uz: Vec = h["nodes"][0]["displacement"].as_array().unwrap().iter().map(|d| d[2].as_f64().unwrap()).collect(); + // The tip starts at zero moving up: it crosses zero upwards once per period after that. + let ups = crossings(&time, &uz); + let t = (ups[ups.len() - 1] - ups[0]) / (ups.len() - 1) as f64; + eprintln!( + "cantilever {n:?}, {hourglass} hourglass: explicit period {t:.5} vs modal {period:.5} ({:+.3}%), one-point estimate {one_point:.5} ({:+.3}%), Euler–Bernoulli {euler:.5}; energy balance {:.1e}, {} increments", + 100.0 * (t / period - 1.0), + 100.0 * (t / one_point - 1.0), + r["energy_balance"].as_f64().unwrap(), + r["increments"] + ); + assert!(ups.len() >= 4, "{ups:?}"); + assert!(r["energy_balance"].as_f64().unwrap() < 0.01, "{r}"); + assert!(r["warnings"].as_array().unwrap().iter().all(|w| w.as_str().unwrap().contains("unit system")), "{r}"); + measured.push(t); + } + // Viscous control adds no stiffness: the one-point estimate holds closely. + assert!((measured[0] / one_point - 1.0).abs() < 0.005, "{measured:?} vs {one_point}"); + // The stiffness form lies between it and the modal period, and within 3% of the latter. + assert!(measured[1] > period && measured[1] < measured[0], "{measured:?}"); + assert!((measured[1] / period - 1.0).abs() < 0.03, "{measured:?} vs {period}"); + } + + /// A small impact model through commands: a hex8 bar moving onto a rigid wall. + fn explicit_impact(s: &mut Session) { + explicit_cantilever(s, 4.0, [8, 1]); + let step = json!({ + "type": "explicit_dynamic", "end_time": 0.2, "history_interval": 0.002, "output_interval": 0.05, "history_nodes": "tip", + "initial_velocities": [{"set": "all_nodes", "velocity": [-0.1, 0, 0]}], + "rigid_walls": [{"set": "x0", "point": [0, 0, 0], "normal": [1, 0, 0]}] + }); + let nodes: Vec = s.model.nodes().keys().copied().collect(); + let script = json!([ + {"command": "model.set.create", "params": {"name": "all_nodes", "nodes": nodes}}, + {"command": "model.section.edit", "params": {"name": "beam", "set": "all", "material": "Mat", "kind": {"type": "solid"}, "formulation": {"integration": "reduced"}}}, + {"command": "model.step.create", "params": {"name": "Impact", "procedure": step}}, + ]); + s.run_script(&script).unwrap(); + } + + #[test] + fn explicit_results_are_kept_probed_and_guarded() { + let mut s = Session::default(); + explicit_impact(&mut s); + let journal = s.journal.len(); + let r = s.execute("solve.explicit", &json!({"step": "Impact"})).unwrap(); + assert_eq!((r["completed"].as_bool(), r["frames"].as_u64()), (Some(true), Some(5)), "{r}"); + // Eight elements, so eight controlling ones; on this coarse mesh the impact excites modes + // near the stability limit, whose energy the integer-step account measures only roughly. + assert_eq!(r["controlling_elements"].as_array().unwrap().len(), 8, "{r}"); + assert!(r["energy_balance"].as_f64().unwrap() < 0.1, "{r}"); + assert_eq!(s.journal.len(), journal, "solving does not change the model"); + let probe = s.execute("solve.probe", &json!({"step": "Impact", "nodes": [1, 9], "frame": 0})).unwrap(); + assert_eq!(probe["nodes"][1]["velocity"][0].as_f64(), Some(-0.1), "{probe}"); + let last = s.execute("solve.probe", &json!({"step": "Impact", "nodes": [9]})).unwrap(); + assert_eq!(last["frame"], 4, "{last}"); + let h = s.execute("solve.history", &json!({"step": "Impact", "every": 2})).unwrap(); + assert_eq!(h["time"].as_array().unwrap().len(), h["energies"].as_array().unwrap().len()); + assert!(h["nodes"][0]["velocity"].as_array().unwrap().len() > 5, "{h}"); + assert_eq!(s.solved_steps().collect::>(), ["Impact"]); + assert!(s.execute("solve.history", &json!({"step": "Modes"})).is_err(), "not solved"); + + let hostile = [ + json!(null), + json!({"step": 7}), + json!({"step": "Nope"}), + json!({"step": "Modes"}), + json!({"step": "Impact", "max_dofs": 0}), + json!({"step": "Impact", "max_dofs": 5}), + json!({"step": "Impact", "max_increments": 0}), + json!({"step": "Impact", "max_increments": 3}), + json!({"step": "Impact", "max_increments": -1}), + json!({"step": "Impact", "max_increments": 1e300}), + json!({"step": "Impact", "max_output_values": 1}), + json!({"step": "Impact", "every": 0}), + json!({"step": "Impact", "every": -3}), + json!({"step": "Impact", "every": 1e300}), + json!({"step": "Impact", "every": u64::MAX}), + json!({"step": "Impact", "nodes": [1], "frame": 99}), + json!({"step": "Impact", "nodes": [1], "frame": -1}), + json!({"step": "Impact", "nodes": [0, 4294967296u64]}), + ]; + for id in ["solve.explicit", "solve.history", "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}"); + } + } + } + // Hostile steps are refused when they are defined; solving one that a later edit makes + // impossible (a node behind its wall) is an error, not a panic. + for bad in [ + json!({"type": "explicit_dynamic", "end_time": -1}), + json!({"type": "explicit_dynamic", "end_time": 1, "dt_scale": 2}), + json!({"type": "explicit_dynamic", "end_time": 1, "max_increments": 0}), + json!({"type": "explicit_dynamic", "end_time": 1, "output_interval": 1e-300}), + json!({"type": "explicit_dynamic", "end_time": 1e308, "history_interval": 1e-308}), + json!({"type": "explicit_dynamic", "end_time": 1, "rigid_walls": [{"set": "x0", "point": [0, 0, 0], "normal": [0, 0, 0]}]}), + json!({"type": "explicit_dynamic", "end_time": 1, "initial_velocities": [{"set": "nope", "velocity": [1, 0, 0]}]}), + json!({"type": "explicit_dynamic", "end_time": 1, "hourglass": {"type": "magic"}}), + ] { + assert!(s.execute("model.step.create", &json!({"name": "Bad", "procedure": bad})).is_err(), "{bad}"); + } + let behind = json!({"type": "explicit_dynamic", "end_time": 0.1, "rigid_walls": [{"set": "x0", "point": [1, 0, 0], "normal": [1, 0, 0]}]}); + s.execute("model.step.create", &json!({"name": "Behind", "procedure": behind})).unwrap(); + assert!(s.execute("solve.explicit", &json!({"step": "Behind"})).unwrap_err().to_string().contains("behind")); + let long = json!({"type": "explicit_dynamic", "end_time": 1e6}); + s.execute("model.step.create", &json!({"name": "Long", "procedure": long})).unwrap(); + assert!(s.execute("solve.explicit", &json!({"step": "Long"})).unwrap_err().to_string().contains("increments")); + assert!(s.execute("model.set.delete", &json!({"name": "all_nodes"})).unwrap_err().to_string().contains("step `Impact`")); + // Defining steps edited the model, so the impact result is stale. + assert_eq!(s.solved_steps().count(), 0); + } } diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index 6f2f318..2fc9325 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -21,10 +21,10 @@ pub mod sample; use std::collections::BTreeMap; use std::sync::Arc; +use modelcraft_explicit::ExplicitResult; 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; @@ -47,6 +47,8 @@ pub enum EngineError { Deck(#[from] modelcraft_decks::Error), #[error(transparent)] Solver(#[from] modelcraft_solver::Error), + #[error(transparent)] + Explicit(#[from] modelcraft_explicit::Error), /// A command panicked; the model is unchanged. #[error("internal error in {command}: {message}")] Internal { command: String, message: String }, @@ -54,6 +56,30 @@ pub enum EngineError { pub type Result = std::result::Result; +/// The result of a solved step, from the implicit or the explicit solver. +#[derive(Debug, Clone, PartialEq)] +pub enum StepResult { + Implicit(modelcraft_solver::StepResult), + Explicit(ExplicitResult), +} + +impl StepResult { + pub fn step(&self) -> &str { + match self { + StepResult::Implicit(r) => r.step(), + StepResult::Explicit(r) => &r.step, + } + } + + /// The node labels the result's per-node values follow, sorted. + pub fn nodes(&self) -> &[u32] { + match self { + StepResult::Implicit(r) => r.nodes(), + StepResult::Explicit(r) => &r.nodes, + } + } +} + /// What a command works on, which decides how [`Session::execute`] records it. #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub enum CommandKind { diff --git a/crates/explicit/Cargo.toml b/crates/explicit/Cargo.toml new file mode 100644 index 0000000..8221ae0 --- /dev/null +++ b/crates/explicit/Cargo.toml @@ -0,0 +1,15 @@ +[package] +name = "modelcraft-explicit" +description = "ModelCraft's explicit dynamics solver: central-difference time integration with lumped mass, one-point hexahedra with hourglass control, nodal-pressure tetrahedra, bulk viscosity and energy accounting" +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/explicit/src/elements.rs b/crates/explicit/src/elements.rs new file mode 100644 index 0000000..a762570 --- /dev/null +++ b/crates/explicit/src/elements.rs @@ -0,0 +1,213 @@ +//! Element internal forces and critical time increments. +//! +//! **One-point hexahedron** (D. P. Flanagan, T. Belytschko, IJNME 17(5):679–706, 1981): the +//! stress of the mean strain gives the forces f_iI = σ_ij b_jI ([`modelcraft_fem::uniform`]). +//! The hourglass modes are resisted through the shape vectors γ_α: with the hourglass rates +//! q̇_iα = Σ_I v_iI γ_αI (or displacements q_iα = Σ_I u_iI γ_αI for the stiffness form), the +//! generalized forces Q_iα = C q̇_iα (viscous) or k q_iα (stiffness) give the nodal forces +//! f_iI = Σ_α Q_iα γ_αI. Because the γ_α are orthogonal to linear fields, these forces do no work +//! in a constant strain state. The coefficients are scaled from the dimensionless hourglass +//! coefficient κ with the element's size h = V / √(2 λ_max(Σ b bᵀ)) (the edge, for a cube), +//! density ρ and wave speed c: C = κ ρ c h² / 16 damps a cube's hourglass modes at a fraction κ +//! of critical at the frequency 2c/h, and k = κ ρ c² h / 16 puts their frequency at √κ · 2c/h. +//! (These normalizations are ours; the paper leaves the scale of the coefficient to the user.) +//! +//! **Average nodal pressure tetrahedron** (J. Bonet, A. J. Burton, *A simple average nodal +//! pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit +//! applications*, Commun. Numer. Methods Eng. 14(5):437–449, 1998): the volumetric strain of each +//! node is the volume-weighted mean of its tetrahedra's, and each tetrahedron takes the mean of +//! its nodes' as its own: σ = 2μ dev ε + K θ̄ I. The forces are those of the strain energy +//! Σ_e V_e μ ‖dev ε_e‖² + ½ K Σ_a V_a θ_a², so the element conserves energy, and it does not lock +//! when the material is nearly incompressible. +//! +//! **Bulk viscosity** (J. von Neumann, R. D. Richtmyer, J. Appl. Phys. 21(3):232–237, 1950; M. L. +//! Wilkins, J. Comput. Phys. 36(3):281–303, 1980): in compression (volumetric strain rate +//! ε̇ᵥ < 0) a pressure q = ρ L (C_q L ε̇ᵥ² − C_l c ε̇ᵥ) spreads shocks over a few elements and +//! damps the ringing behind wave fronts. It lowers the stable increment to +//! Δt = L / (Q + √(Q² + c²)), Q = C_l c + C_q L |ε̇ᵥ| (Wilkins 1980). +//! +//! **Critical increment.** Each element's stable length L = 2c/ω_max comes from the bound on its +//! highest frequency ([`modelcraft_fem::uniform::UniformStrain::frequency_bound`], plus the +//! hourglass stiffness when it is used), so without bulk viscosity its critical increment is +//! 2/ω_max (T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, 2014, §6.3; the element eigenvalue +//! inequality makes the smallest of them stable for the mesh). For the nodal-pressure +//! tetrahedron the averaging only lowers the volumetric energy (the mean of squares bounds the +//! square of the mean), so the plain tetrahedron's bound holds. + +use modelcraft_fem::material::stress; +use modelcraft_fem::uniform::UniformStrain; +use modelcraft_model::{BulkViscosity, HourglassKind}; + +use crate::system::{Material, System}; + +/// Scratch vectors and the element state of the latest evaluation. +#[derive(Debug, Clone)] +pub struct Workspace { + /// Internal forces per degree of freedom, by origin. + pub stress_force: Vec, + pub hourglass_force: Vec, + pub viscous_force: Vec, + /// Per element (in [`System::elements`] order): the material stress and the strain. + pub stress: Vec<[f64; 6]>, + pub strain: Vec<[f64; 6]>, + /// Per node: the mean volumetric strain of its tetrahedra (a volume-weighted sum until it is + /// divided by `tet_volume`). + theta: Vec, + tet_volume: Vec, +} + +impl Workspace { + pub fn new(sys: &System) -> Workspace { + let n = sys.labels.len(); + let mut tet_volume = vec![0.0; n]; + for t in &sys.tets { + for node in t.nodes { + if let Some(v) = tet_volume.get_mut(node) { + *v += t.kinematics.volume; + } + } + } + let dofs = sys.dofs(); + let elements = sys.elements.len(); + Workspace { + stress_force: vec![0.0; dofs], + hourglass_force: vec![0.0; dofs], + viscous_force: vec![0.0; dofs], + stress: vec![[0.0; 6]; elements], + strain: vec![[0.0; 6]; elements], + theta: vec![0.0; n], + tet_volume, + } + } +} + +pub(crate) fn gather(x: &[f64], nodes: &[usize; N]) -> [[f64; 3]; N] { + nodes.map(|n| [0, 1, 2].map(|k| x.get(3 * n + k).copied().unwrap_or(0.0))) +} + +/// Add a node's vector into a per-degree-of-freedom vector. +pub(crate) fn add_at(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; + } + } +} + +fn scatter(f: &mut [f64], nodes: &[usize; N], v: &[[f64; 3]; N]) { + for (n, vi) in nodes.iter().zip(v) { + add_at(f, *n, vi); + } +} + +/// The bulk viscosity pressure (positive in compression) for a volumetric strain rate. +pub fn bulk_pressure(rate: f64, density: f64, length: f64, wave_speed: f64, b: &BulkViscosity) -> f64 { + if rate < 0.0 { density * length * (b.quadratic * length * rate * rate - b.linear * wave_speed * rate) } else { 0.0 } +} + +/// The critical increment of an element of stable length `length` at a volumetric strain rate +/// (Wilkins 1980). +pub fn critical_increment(rate: f64, length: f64, wave_speed: f64, b: &BulkViscosity) -> f64 { + let q = if rate < 0.0 { b.linear * wave_speed + b.quadratic * length * rate.abs() } else { 0.0 }; + length / (q + (q * q + wave_speed * wave_speed).sqrt()) +} + +/// What the hexahedron and the tetrahedron share once their stress is known: the stress and bulk +/// viscosity forces, the element state and the critical increment. +struct Finish<'a, const N: usize> { + kinematics: &'a UniformStrain, + nodes: &'a [usize; N], + index: usize, + length: f64, +} + +impl Finish<'_, N> { + fn apply(&self, ws: &mut Workspace, mat: &Material, v: &[f64], sigma: [f64; 6], strain: [f64; 6], bulk: &BulkViscosity) -> f64 { + let k = self.kinematics; + let rate = k.divergence(&gather(v, self.nodes)); + let q = bulk_pressure(rate, mat.density, self.length, mat.wave_speed, bulk); + scatter(&mut ws.stress_force, self.nodes, &k.forces(&sigma)); + if q != 0.0 { + scatter(&mut ws.viscous_force, self.nodes, &k.pressure_forces(q)); + } + if let Some(s) = ws.stress.get_mut(self.index) { + *s = sigma; + } + if let Some(s) = ws.strain.get_mut(self.index) { + *s = strain; + } + critical_increment(rate, self.length, mat.wave_speed, bulk) + } +} + +/// Evaluate the internal forces at displacements `u` and velocities `v` into `ws` and return the +/// smallest critical increment of the elements. +pub fn internal_forces(sys: &System, u: &[f64], v: &[f64], bulk: &BulkViscosity, ws: &mut Workspace) -> f64 { + for f in [&mut ws.stress_force, &mut ws.hourglass_force, &mut ws.viscous_force] { + f.fill(0.0); + } + let mut critical = f64::INFINITY; + for h in &sys.hexes { + let Some(mat) = sys.materials.get(h.material) else { continue }; + let ue = gather(u, &h.nodes); + let strain = h.kinematics.strain(&ue); + let sigma = stress(&mat.elasticity, &strain); + // Hourglass forces on the displacement (stiffness) or velocity (viscous) modes. + let modes = match sys.hourglass { + HourglassKind::Viscous => gather(v, &h.nodes), + HourglassKind::Stiffness => ue, + }; + let mut f = [[0.0; 3]; 8]; + for g in &h.gamma { + let mut generalized = [0.0; 3]; + for (x, gk) in modes.iter().zip(g) { + for i in 0..3 { + generalized[i] += x[i] * gk; + } + } + for (fk, gk) in f.iter_mut().zip(g) { + for i in 0..3 { + fk[i] += h.hourglass * generalized[i] * gk; + } + } + } + scatter(&mut ws.hourglass_force, &h.nodes, &f); + let finish = Finish { kinematics: &h.kinematics, nodes: &h.nodes, index: h.index, length: h.length }; + critical = critical.min(finish.apply(ws, mat, v, sigma, strain, bulk)); + } + if !sys.tets.is_empty() { + // Pass 1: each tetrahedron's strain, and each node's mean volumetric strain. + ws.theta.fill(0.0); + for t in &sys.tets { + let strain = t.kinematics.strain(&gather(u, &t.nodes)); + for node in t.nodes { + if let Some(x) = ws.theta.get_mut(node) { + *x += t.kinematics.volume * (strain[0] + strain[1] + strain[2]); + } + } + if let Some(s) = ws.strain.get_mut(t.index) { + *s = strain; + } + } + for (x, v) in ws.theta.iter_mut().zip(&ws.tet_volume) { + if *v > 0.0 { + *x /= v; + } + } + // Pass 2: each tetrahedron's pressure from the mean of its nodes' volumetric strains. + for t in &sys.tets { + let Some(mat) = sys.materials.get(t.material) else { continue }; + let strain = ws.strain.get(t.index).copied().unwrap_or_default(); + let averaged = t.nodes.iter().map(|n| ws.theta.get(*n).copied().unwrap_or(0.0)).sum::() / 4.0; + let own = strain[0] + strain[1] + strain[2]; + let mut sigma = stress(&mat.elasticity, &strain); + let shift = mat.bulk.unwrap_or(0.0) * (averaged - own); + for s in sigma.iter_mut().take(3) { + *s += shift; + } + let finish = Finish { kinematics: &t.kinematics, nodes: &t.nodes, index: t.index, length: t.length }; + critical = critical.min(finish.apply(ws, mat, v, sigma, strain, bulk)); + } + } + critical +} diff --git a/crates/explicit/src/integrate.rs b/crates/explicit/src/integrate.rs new file mode 100644 index 0000000..fabe939 --- /dev/null +++ b/crates/explicit/src/integrate.rs @@ -0,0 +1,343 @@ +//! The time loop: central differences with variable increments, kinematic constraints, the +//! energy account and output. +//! +//! One increment (Belytschko, Liu, Moran, Elkhodary 2014, Box 6.1, with Δt^n = ½(Δt^{n−½} + +//! Δt^{n+½})): +//! 1. u^n = u^{n−1} + Δt^{n−½} v^{n−½}; t^n = t^{n−1} + Δt^{n−½}; +//! 2. internal forces f^n from u^n (stresses, hourglass stiffness) and v^{n−½} (hourglass and +//! bulk viscosity), and the next increment Δt^{n+½} = scale · min_e Δt_e; +//! 3. v^{n+½} = v^{n−½} + Δt^n M⁻¹(f_ext − f^n), then the constraints: rigid walls set the normal +//! velocity of a node that would cross its wall so it lands on it (a frictionless kinematic +//! contact condition: T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier, W. +//! Kanoknukulchai, *A finite element method for a class of contact-impact problems*, CMAME +//! 8(3):249–276, 1976), and prescribed degrees of freedom take their velocity. The change a +//! constraint makes is the impulse of its reaction, r^n = M Δv / Δt^n; +//! 4. energies: the works of the internal and external forces by the trapezoidal rule, +//! W += ½ Δu · (f^{n−1} + f^n) (exactly the strain energy for a linear elastic material), +//! the constraints' work Δt^n r^n · ½(v^{n−½} + v^{n+½}) (exactly the kinetic energy their +//! impulses add or remove), and the kinetic energy of v^n = v^{n−½} + ½Δt^{n−½} a^n. +//! +//! The run starts from u = 0 and the initial velocities with the half increment +//! v^{½} = v^0 + ½Δt^{½} a^0. + +use modelcraft_model::{ExplicitDynamic, Model, OutputVariable, Procedure, Step}; + +use crate::elements::{Workspace, internal_forces}; +use crate::result::{ControllingElement, Energies, ExplicitResult, Frame, History}; +use crate::system::System; +use crate::{Error, Limits, Result}; + +/// How many elements the result names as controlling the time increment. +pub const CONTROLLING: usize = 10; + +/// Relative energy balance error above which a run warns. +pub const BALANCE_TOLERANCE: f64 = 0.01; + +/// Hourglass energy as a fraction of the largest internal energy of the run above which it warns +/// (the internal energy may fall back to zero, as in a rebound). +pub const HOURGLASS_TOLERANCE: f64 = 0.1; + +/// The step named `name`, or the model's only explicit dynamic step. +fn select_step<'m>(model: &'m Model, name: Option<&str>) -> Result<(&'m Step, &'m ExplicitDynamic)> { + let step = model.select_step(name, Procedure::is_dynamic, "explicit dynamic")?; + match &step.procedure { + Procedure::ExplicitDynamic(e) => Ok((step, e)), + _ => Err(Error::Internal("an explicit dynamic step".into())), + } +} + +/// Times at a fixed interval: whether `t` has reached the next one. +struct Schedule { + interval: f64, + next: f64, +} + +impl Schedule { + fn due(&mut self, t: f64) -> bool { + if t < self.next * (1.0 - 1e-12) { + return false; + } + self.next = ((t / self.interval).floor() + 1.0) * self.interval; + true + } +} + +/// What the run writes. +struct Output { + displacement: bool, + stress: bool, + strain: bool, +} + +/// Velocity constraints over an increment of length `dt` (rigid walls when dt > 0, then +/// prescribed degrees of freedom). +fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64) { + if dt > 0.0 { + for w in &sys.walls { + for &node in &w.nodes { + let x = sys.xyz.get(node).copied().unwrap_or_default(); + let at = |k: usize| 3 * node + k; + let gap: f64 = (0..3).map(|k| (x[k] + u.get(at(k)).copied().unwrap_or(0.0) - w.point[k]) * w.normal[k]).sum(); + let normal: f64 = (0..3).map(|k| v.get(at(k)).copied().unwrap_or(0.0) * w.normal[k]).sum(); + if gap + dt * normal < 0.0 { + let change = -gap / dt - normal; + for k in 0..3 { + if let Some(c) = v.get_mut(at(k)) { + *c += change * w.normal[k]; + } + } + } + } + } + } + for (c, p) in v.iter_mut().zip(&sys.prescribed) { + if let Some(p) = p { + *c = *p; + } + } +} + +/// v_out = v_in + dt M⁻¹ (f_ext − f_int). +fn accelerate(sys: &System, ws: &Workspace, mass: &[f64], v_in: &[f64], dt: f64, v_out: &mut [f64]) { + let internal = ws.stress_force.iter().zip(&ws.hourglass_force).zip(&ws.viscous_force); + for ((((out, v), m), e), ((s, h), q)) in v_out.iter_mut().zip(v_in).zip(mass).zip(&sys.external).zip(internal) { + *out = v + dt * ((e - (s + h + q)) / m); + } +} + +/// The constraint reactions M (v_constrained − v_free) / Δt, into `r`. +fn impulses(mass: &[f64], free: &[f64], constrained: &[f64], dt: f64, r: &mut [f64]) { + for (r, (m, (f, c))) in r.iter_mut().zip(mass.iter().zip(free.iter().zip(constrained))) { + *r = if c == f { 0.0 } else { m * (c - f) / dt }; + } +} + +fn per_node(x: &[f64]) -> Vec<[f64; 3]> { + x.as_chunks::<3>().0.to_vec() +} + +/// Solve an explicit dynamic step: the step named `step`, or the model's only one. +pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Result { + let resolved = model.resolve()?; + let (step, settings) = select_step(model, step)?; + let sys = System::new(model, &resolved, step, settings, limits.max_dofs)?; + let end = settings.end_time; + let scale = settings.dt_scale; + let bulk = settings.bulk_viscosity; + let max_increments = settings.max_increments.min(limits.max_increments); + let mut warnings = sys.warnings.clone(); + + // The controlling elements: the smallest critical increments at rest. + let mut critical: Vec = sys + .hexes + .iter() + .map(|h| (h.label, h.length, h.material)) + .chain(sys.tets.iter().map(|t| (t.label, t.length, t.material))) + .map(|(label, length, m)| ControllingElement { label, increment: length / sys.materials.get(m).map_or(f64::INFINITY, |m| m.wave_speed) }) + .collect(); + critical.sort_by(|a, b| a.increment.total_cmp(&b.increment).then(a.label.cmp(&b.label))); + critical.truncate(CONTROLLING); + let smallest = critical + .first() + .map(|c| c.increment) + .filter(|d| *d > 0.0 && d.is_finite()) + .ok_or_else(|| Error::Invalid("no element has a finite, positive critical time increment (check the material constants)".into()))?; + let estimate = end / (scale * smallest); + if estimate > max_increments as f64 { + let first = critical.first().map_or(0, |c| c.label); + return Err(Error::Limit(format!( + "the step needs about {estimate:.3e} increments of {:.3e} (element {first} controls it), more than the cap of {max_increments}", + scale * smallest + ))); + } + + // Output sizes. + let output = Output { + displacement: step.output.contains(&OutputVariable::Displacement), + stress: step.output.contains(&OutputVariable::Stress), + strain: step.output.contains(&OutputVariable::Strain), + }; + let (nodes, elements) = (sys.labels.len(), sys.elements.len()); + let frames = (end / settings.output_interval()).ceil() + 2.0; + let samples = (end / settings.history_interval()).ceil() + 2.0; + let per_frame = (if output.displacement { 6 * nodes } else { 0 }) + 6 * elements * (usize::from(output.stress) + usize::from(output.strain)); + let per_sample = 16 + 6 * sys.history_nodes.len(); + let values = frames * per_frame as f64 + samples * per_sample as f64; + if values > limits.max_output_values as f64 { + return Err(Error::Limit(format!( + "the output would hold {values:.3e} numbers; the limit is {} (lengthen the output or history interval, or request less output)", + limits.max_output_values + ))); + } + + let n = sys.dofs(); + let mass: Vec = (0..n).map(|d| sys.mass.get(d / 3).copied().unwrap_or(0.0)).collect(); + if mass.iter().any(|m| !(*m > 0.0 && m.is_finite())) { + return Err(Error::Internal("a node has no positive mass".into())); + } + let total_mass: f64 = sys.mass.iter().sum(); + let mut ws = Workspace::new(&sys); + let mut previous = Workspace::new(&sys); + let mut u = vec![0.0; n]; + let mut du = vec![0.0; n]; + let v0 = sys.initial_velocity.clone(); + let initial_kinetic = 0.5 * v0.iter().zip(&mass).map(|(v, m)| m * v * v).sum::(); + + // Start: forces at rest with the initial velocities, then the half increment. + let crit = internal_forces(&sys, &u, &v0, &bulk, &mut ws); + let mut dt_half = (scale * crit).min(end); + let initial_increment = dt_half; + let mut v_half = vec![0.0; n]; + accelerate(&sys, &ws, &mass, &v0, 0.5 * dt_half, &mut v_half); + let mut free = v_half.clone(); + constrain(&sys, &u, &mut v_half, dt_half); + let mut reaction = vec![0.0; n]; + impulses(&mass, &free, &v_half, 0.5 * dt_half, &mut reaction); + // The constraints' work is that of their impulses on the mean of the velocities before + // and after the increment, Δt r · ½(v^{n−½} + v^{n+½}): exact for the kinetic energy a + // kinematic constraint removes or adds (a node stopped by a wall loses ½ m v²). + let impulse_work = |r: &[f64], before: &[f64], after: &[f64], dt: f64| -> f64 { + dt * r.iter().zip(before.iter().zip(after)).map(|(r, (a, b))| r * 0.5 * (a + b)).sum::() + }; + let mut energy = Energies { kinetic: initial_kinetic, constraint: impulse_work(&reaction, &v0, &v_half, 0.5 * dt_half), ..Energies::default() }; + + let mut history = History { nodes: sys.history_nodes.iter().filter_map(|i| sys.labels.get(*i).copied()).collect(), ..History::default() }; + let mut out_frames = Vec::new(); + let sample = |history: &mut History, t: f64, dt: f64, e: Energies, u: &[f64], v: &[f64], r: &[f64]| { + history.time.push(t); + history.increment.push(dt); + history.energies.push(e); + let sum = |x: &[f64], weight: &dyn Fn(usize) -> f64| { + let mut s = [0.0; 3]; + for (node, c) in x.as_chunks::<3>().0.iter().enumerate() { + for k in 0..3 { + s[k] += weight(node) * c[k]; + } + } + s + }; + history.momentum.push(sum(v, &|node| sys.mass.get(node).copied().unwrap_or(0.0))); + history.constraint_force.push(sum(r, &|_| 1.0)); + let pick = |x: &[f64]| sys.history_nodes.iter().map(|i| [0, 1, 2].map(|k| x.get(3 * i + k).copied().unwrap_or(0.0))).collect(); + history.displacement.push(pick(u)); + history.velocity.push(pick(v)); + }; + let frame = |t: f64, u: &[f64], v: &[f64], ws: &Workspace| Frame { + time: t, + displacement: if output.displacement { per_node(u) } else { Vec::new() }, + velocity: if output.displacement { per_node(v) } else { Vec::new() }, + stress: if output.stress { ws.stress.clone() } else { Vec::new() }, + strain: if output.strain { ws.strain.clone() } else { Vec::new() }, + }; + let mut histories = Schedule { interval: settings.history_interval(), next: settings.history_interval() }; + let mut frames_due = Schedule { interval: settings.output_interval(), next: settings.output_interval() }; + sample(&mut history, 0.0, dt_half, energy, &u, &v0, &reaction); + out_frames.push(frame(0.0, &u, &v0, &ws)); + + let mut t = 0.0; + let mut increments = 0u64; + let mut peak_internal: f64 = 0.0; + let mut v_now = v0; + let mut completed = true; + let mut v_new = vec![0.0; n]; + loop { + if increments >= max_increments { + completed = false; + warnings.push(format!("stopped at time {t:.6e} after {increments} increments (the cap); the step ends at {end:.6e}")); + break; + } + for ((x, dx), v) in u.iter_mut().zip(du.iter_mut()).zip(&v_half) { + *dx = dt_half * v; + *x += *dx; + } + t += dt_half; + increments += 1; + std::mem::swap(&mut ws, &mut previous); + let crit = internal_forces(&sys, &u, &v_half, &bulk, &mut ws); + let remaining = end - t; + let last = remaining <= 1e-12 * end; + let dt_next = if last { 0.0 } else { (scale * crit).min(remaining) }; + let dt_n = 0.5 * (dt_half + dt_next); + accelerate(&sys, &ws, &mass, &v_half, dt_n, &mut v_new); + free.clone_from(&v_new); + constrain(&sys, &u, &mut v_new, dt_next); + impulses(&mass, &free, &v_new, dt_n, &mut reaction); + // One pass over the degrees of freedom: v^n = v^{n−½} + ½Δt^{n−½} a^n with the constrained + // a^n = (v^{n+½} − v^{n−½}) / Δt^n, and the works (trapezoidal over the displacement for + // the internal and external forces). + let weight = 0.5 * dt_half / dt_n; + let [mut internal, mut hourglass, mut viscous, mut external, mut constraint, mut kinetic] = [0.0; 6]; + let forces = ws.stress_force.iter().zip(&ws.hourglass_force).zip(&ws.viscous_force); + let before = previous.stress_force.iter().zip(&previous.hourglass_force).zip(&previous.viscous_force); + let state = v_now.iter_mut().zip(&v_half).zip(&v_new).zip(&mass).zip(&du).zip(&sys.external).zip(&reaction); + for ((((((((now, vh), vn), m), dx), e), r), ((s, h), q)), ((s0, h0), q0)) in state.zip(forces).zip(before) { + *now = vh + weight * (vn - vh); + internal += dx * (s0 + s); + hourglass += dx * (h0 + h); + viscous += dx * (q0 + q); + external += dx * e; + constraint += r * 0.5 * (vh + vn); + kinetic += m * *now * *now; + } + energy.internal += 0.5 * internal; + energy.hourglass += 0.5 * hourglass; + energy.viscous += 0.5 * viscous; + energy.external += external; + energy.constraint += dt_n * constraint; + energy.kinetic = 0.5 * kinetic; + peak_internal = peak_internal.max(energy.internal.abs()); + energy.balance = energy.kinetic + energy.internal + energy.hourglass + energy.viscous - energy.external - energy.constraint - initial_kinetic; + if !(energy.balance.is_finite() && crit.is_finite()) { + return Err(Error::Unstable(format!("the solution is not finite at time {t:.6e} (increment {increments})"))); + } + if histories.due(t) || last { + sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction); + } + if frames_due.due(t) || last { + out_frames.push(frame(t, &u, &v_now, &ws)); + } + std::mem::swap(&mut v_half, &mut v_new); + dt_half = dt_next; + if last { + break; + } + } + if !completed { + if history.time.last() != Some(&t) { + sample(&mut history, t, dt_half, energy, &u, &v_now, &reaction); + } + if out_frames.last().map(|f| f.time) != Some(t) { + out_frames.push(frame(t, &u, &v_now, &ws)); + } + } + + let error = energy.relative_error(initial_kinetic); + if error > BALANCE_TOLERANCE { + warnings.push(format!("the energy balance is off by {:.2}% of the largest energy (kinetic, internal or work)", 100.0 * error)); + } + // Round-off hourglass energy in a model at rest or in rigid motion is not worth a warning. + let floor = + 1e-9 * [initial_kinetic, energy.kinetic, peak_internal, energy.external.abs(), energy.constraint.abs()].into_iter().fold(0.0, f64::max); + if energy.hourglass.abs() > HOURGLASS_TOLERANCE * peak_internal && energy.hourglass.abs() > floor { + warnings.push(format!( + "hourglass energy is {:.1}% of the largest internal energy; refine the mesh or change the hourglass control", + 100.0 * energy.hourglass.abs() / peak_internal.max(f64::MIN_POSITIVE) + )); + } + Ok(ExplicitResult { + step: step.name.clone(), + nodes: sys.labels.clone(), + elements: sys.elements.clone(), + frames: out_frames, + history, + controlling: critical, + initial_increment, + increments, + time: t, + completed, + total_mass, + energy, + initial_kinetic, + warnings, + }) +} diff --git a/crates/explicit/src/lib.rs b/crates/explicit/src/lib.rs new file mode 100644 index 0000000..4dbe072 --- /dev/null +++ b/crates/explicit/src/lib.rs @@ -0,0 +1,92 @@ +//! ModelCraft's explicit dynamics solver: transient response of a model's explicit dynamic step +//! by central-difference time integration. +//! +//! The method is the explicit central-difference scheme of T. Belytschko, W. K. Liu, B. Moran, +//! K. Elkhodary, *Nonlinear Finite Elements for Continua and Structures*, 2nd ed., Wiley 2014 +//! (ch. 6, the flowchart of Box 6.1): with a diagonal (lumped) mass M, +//! a = M⁻¹(f_ext − f_int), v^{n+½} = v^{n−½} + Δt^n a^n, u^{n+1} = u^n + Δt^{n+½} v^{n+½}, +//! so a time increment is one pass over the elements and a division per degree of freedom. +//! +//! - **System** ([`system`]): the model is resolved; one-point hexahedra and nodal-pressure +//! tetrahedra with their materials, row-sum lumped masses, prescribed velocities, rigid walls, +//! initial velocities and the constant load vector of the step. +//! - **Elements** ([`elements`]): the uniform-strain hexahedron with Flanagan–Belytschko +//! hourglass control (viscous or stiffness), the Bonet–Burton average nodal pressure +//! tetrahedron, von Neumann–Richtmyer/Wilkins bulk viscosity, and the element critical time +//! increments. +//! - **Integration** ([`integrate`]): the time loop, kinematic constraints, energy accounting +//! with a balance check, history samples and field frames ([`result`]). +//! +//! This first version is geometrically linear (small displacements and strains; the elements' +//! gradients are those of the undeformed mesh) with linear elastic materials, which is what the +//! analytic verification problems need. Large deformation kinematics with objective stress rates +//! and inelastic materials come next (docs/solvers/literature.md, Part F, X2). +//! +//! Everything is sequential, with elements in label order and no parallel reductions, so the +//! same model gives the same numbers bit for bit; sizes, increments and output are capped by +//! [`Limits`]. +#![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] +#![forbid(unsafe_code)] +// Small fixed-size vector kernels read as the textbook formulas with index loops. +#![allow(clippy::needless_range_loop)] + +pub mod elements; +pub mod integrate; +pub mod result; +pub mod system; + +pub use integrate::solve_explicit; +pub use result::{ControllingElement, Energies, ExplicitResult, Frame, History}; + +#[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, …). + #[error("element {label}: {source}")] + Element { label: u32, source: modelcraft_fem::Error }, + #[error("limit reached: {0}")] + Limit(String), + /// The solution stopped being finite. + #[error("the solution became unstable: {0}")] + Unstable(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 and length of a run. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct Limits { + /// Degrees of freedom of the mesh. + pub max_dofs: usize, + /// Time increments; the step's own cap applies too, whichever is lower. + pub max_increments: u64, + /// Numbers kept in field frames and histories together (8 bytes each). + pub max_output_values: usize, +} + +impl Limits { + pub const MAX_DOFS: usize = 6_000_000; + pub const MAX_INCREMENTS: u64 = modelcraft_model::steps::MAX_INCREMENTS; + pub const MAX_OUTPUT_VALUES: usize = 50_000_000; +} + +impl Default for Limits { + fn default() -> Limits { + Limits { max_dofs: Limits::MAX_DOFS, max_increments: Limits::MAX_INCREMENTS, max_output_values: Limits::MAX_OUTPUT_VALUES } + } +} diff --git a/crates/explicit/src/result.rs b/crates/explicit/src/result.rs new file mode 100644 index 0000000..f4519f7 --- /dev/null +++ b/crates/explicit/src/result.rs @@ -0,0 +1,99 @@ +//! What an explicit dynamic run keeps: histories sampled at the history interval and field +//! frames at the output interval. + +use serde::Serialize; + +/// The energy account at one time (Belytschko, Liu, Moran, Elkhodary 2014, §6.2.3). The works of +/// the internal and external forces are integrated with the trapezoidal rule over each +/// increment's displacement; the constraints' with their impulses (see [`crate::integrate`]). +#[derive(Debug, Clone, Copy, Default, PartialEq, Serialize)] +pub struct Energies { + /// ½ vᵀ M v. + pub kinetic: f64, + /// Work of the material stresses: the strain energy for an elastic material. + pub internal: f64, + /// Work of the hourglass forces (stored in the stiffness form, dissipated in the viscous one). + pub hourglass: f64, + /// Work of the bulk viscosity (dissipated). + pub viscous: f64, + /// Work of the loads. + pub external: f64, + /// Work of the boundary conditions and rigid walls on the model. + pub constraint: f64, + /// kinetic + internal + hourglass + viscous − external − constraint − the initial kinetic + /// energy: zero when energy is conserved. + pub balance: f64, +} + +impl Energies { + /// The balance relative to the largest energy in the account. + pub fn relative_error(&self, initial_kinetic: f64) -> f64 { + let scale = [initial_kinetic, self.kinetic, self.internal + self.hourglass + self.viscous, self.external.abs(), self.constraint.abs()] + .into_iter() + .fold(0.0, f64::max); + if scale > 0.0 { self.balance.abs() / scale } else { 0.0 } + } +} + +/// The fields at one time. Vectors that the step's output does not request are empty. +#[derive(Debug, Clone, PartialEq)] +pub struct Frame { + pub time: f64, + /// Per node in [`ExplicitResult::nodes`] order. + pub displacement: Vec<[f64; 3]>, + pub velocity: Vec<[f64; 3]>, + /// Per element in [`ExplicitResult::elements`] order (one-point elements: one value each): + /// xx, yy, zz, xy, xz, yz, engineering shear strains. + pub stress: Vec<[f64; 6]>, + pub strain: Vec<[f64; 6]>, +} + +/// Quantities sampled through the run. +#[derive(Debug, Clone, Default, PartialEq)] +pub struct History { + pub time: Vec, + /// The time increment taken at each sample. + pub increment: Vec, + pub energies: Vec, + /// Total momentum M v. + pub momentum: Vec<[f64; 3]>, + /// Sum of the forces the boundary conditions and rigid walls exert on the model. + pub constraint_force: Vec<[f64; 3]>, + /// The nodes with histories, and per sample their displacement and velocity (in this order). + pub nodes: Vec, + pub displacement: Vec>, + pub velocity: Vec>, +} + +/// An element that limits the stable time increment. +#[derive(Debug, Clone, Copy, PartialEq, Serialize)] +pub struct ControllingElement { + pub label: u32, + /// Its critical increment at rest (no bulk viscosity), before the scale factor. + pub increment: f64, +} + +/// The result of an explicit dynamic step. +#[derive(Debug, Clone, PartialEq)] +pub struct ExplicitResult { + pub step: String, + /// Node labels with degrees of freedom, sorted. + pub nodes: Vec, + /// Element labels, sorted. + pub elements: Vec, + pub frames: Vec, + pub history: History, + /// The elements with the smallest critical increments (at most ten), smallest first. + pub controlling: Vec, + /// The first time increment taken (the scale times the smallest critical increment). + pub initial_increment: f64, + pub increments: u64, + /// The time reached: the step's end time unless the run stopped at the increment cap. + pub time: f64, + pub completed: bool, + pub total_mass: f64, + /// The energy account at the end. + pub energy: Energies, + pub initial_kinetic: f64, + pub warnings: Vec, +} diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs new file mode 100644 index 0000000..9d2ed7e --- /dev/null +++ b/crates/explicit/src/system.rs @@ -0,0 +1,420 @@ +//! From the resolved model to the explicit system: elements with their materials and stable +//! lengths, lumped masses, prescribed velocities, rigid walls, initial velocities and the load +//! vector of an explicit dynamic step. + +use modelcraft_fem::face; +use modelcraft_fem::material::Elasticity; +use modelcraft_fem::uniform::{Hex8, Tet4, modulus_bound}; +use modelcraft_model::{Elastic, ElementKind, ExplicitDynamic, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, Target}; + +use crate::elements::add_at; +use crate::{Error, Result}; + +/// Most nodes whose histories a step records. +pub const MAX_HISTORY_NODES: usize = 1000; + +/// A material as the elements use it. +#[derive(Debug, Clone)] +pub struct Material { + pub name: String, + pub elasticity: Elasticity, + pub density: f64, + /// Dilatational wave speed √(M/ρ), with M the largest normal stiffness (λ + 2μ when + /// isotropic). + pub wave_speed: f64, + /// The bound of [`modelcraft_fem::uniform::modulus_bound`]. + pub modulus: f64, + /// Bulk modulus, for the nodal-pressure tetrahedron (isotropic materials). + pub bulk: Option, +} + +/// A one-point hexahedron. +#[derive(Debug, Clone)] +pub struct Hex { + pub label: u32, + /// Position in [`System::elements`]. + pub index: usize, + pub nodes: [usize; 8], + pub material: usize, + pub kinematics: Hex8, + pub gamma: [[f64; 8]; 4], + /// The hourglass coefficient: a damping (viscous) or a stiffness (stiffness form), see + /// [`crate::elements`]. + pub hourglass: f64, + /// The element's stable length: its critical increment without bulk viscosity is + /// `length / wave_speed`. + pub length: f64, +} + +/// A nodal-pressure tetrahedron. +#[derive(Debug, Clone)] +pub struct Tet { + pub label: u32, + pub index: usize, + pub nodes: [usize; 4], + pub material: usize, + pub kinematics: Tet4, + pub length: f64, +} + +/// A fixed, frictionless rigid plane. +#[derive(Debug, Clone)] +pub struct Wall { + pub point: [f64; 3], + /// Unit normal, towards the side the nodes stay on. + pub normal: [f64; 3], + pub nodes: Vec, +} + +/// The discrete system of an explicit dynamic step. +#[derive(Debug, Clone)] +pub struct System { + /// Node labels with degrees of freedom (those of the elements), sorted. + pub labels: Vec, + pub xyz: Vec<[f64; 3]>, + /// Element labels, sorted. + pub elements: Vec, + pub materials: Vec, + pub hexes: Vec, + pub tets: Vec, + /// Row-sum lumped mass per node. + pub mass: Vec, + /// Prescribed velocity per degree of freedom (`None` = free). + pub prescribed: Vec>, + pub walls: Vec, + /// Velocity per degree of freedom at the start (prescribed velocities included). + pub initial_velocity: Vec, + /// Constant external force per degree of freedom. + pub external: Vec, + /// Nodes whose histories are recorded. + pub history_nodes: Vec, + pub hourglass: HourglassKind, + pub warnings: Vec, +} + +fn invalid(m: impl Into) -> Error { + Error::Invalid(m.into()) +} + +/// A material's constants for the explicit elements. +fn material(model: &Model, name: &str) -> Result { + let m = model.material(name).ok_or_else(|| invalid(format!("material `{name}` is missing")))?; + let elastic = m.elastic.as_ref().ok_or_else(|| invalid(format!("material `{name}` has no elastic constants")))?; + let density = m + .density + .filter(|d| *d > 0.0 && d.is_finite()) + .ok_or_else(|| invalid(format!("material `{name}` has no density, and explicit dynamics needs one")))?; + let d = modelcraft_fem::elasticity(elastic).map_err(|e| invalid(format!("material `{name}`: {e}")))?; + let stiffest = d[0][0].max(d[1][1]).max(d[2][2]); + let bulk = match *elastic { + Elastic::Isotropic { young, poisson } => Some(young / (3.0 * (1.0 - 2.0 * poisson))), + Elastic::EngineeringConstants { .. } => None, + }; + Ok(Material { + name: name.to_string(), + elasticity: d, + density, + wave_speed: (stiffest / density).sqrt(), + modulus: modulus_bound(elastic, &d), + bulk, + }) +} + +/// The largest absolute row sum of the hourglass vectors' Gram matrix Σ_I γ_αI γ_βI / v_I (an +/// upper bound of its largest eigenvalue, exact for parallelepipeds where the γ_α are +/// orthogonal). +fn hourglass_gram_bound(gamma: &[[f64; 8]; 4], nodal_volume: &[f64; 8]) -> f64 { + let mut worst: f64 = 0.0; + for a in gamma { + let row: f64 = gamma.iter().map(|b| a.iter().zip(b).zip(nodal_volume).map(|((x, y), v)| x * y / v).sum::().abs()).sum(); + worst = worst.max(row); + } + worst +} + +/// The section's elements of a kind the explicit solver has, or why not. +fn check_formulation(kind: ElementKind, section_name: &str, f: &modelcraft_model::Formulation) -> Result<()> { + let refuse = |why: &str| Err(Error::NotSupported(format!("section `{section_name}`: {why}"))); + if f.hybrid { + return refuse("hybrid elements in explicit dynamics"); + } + match kind { + ElementKind::Hex8 if f.incompatible_modes || f.integration == Integration::Full => refuse( + "fully integrated or incompatible-mode hex8 in explicit dynamics (explicit hex8 elements are one-point with hourglass control: set the section's integration to reduced)", + ), + ElementKind::Hex8 => Ok(()), + ElementKind::Tet4 if f.integration == Integration::Reduced => refuse("reduced integration of tet4 elements"), + ElementKind::Tet4 => Ok(()), + other => refuse(&format!("{} elements in explicit dynamics (it has one-point hex8 and nodal-pressure tet4 so far)", other.name())), + } +} + +impl System { + pub fn new(model: &Model, resolved: &Resolved, step: &Step, settings: &ExplicitDynamic, max_dofs: usize) -> Result { + if let Some(name) = model.couplings().keys().next() { + return Err(Error::NotSupported(format!("couplings (coupling `{name}`) in explicit dynamics"))); + } + // Elements by label, with their section's material. + let mut materials: Vec = Vec::new(); + let mut chosen: Vec<(u32, ElementKind, Vec, usize)> = Vec::new(); + for (name, section) in model.sections() { + if !matches!(section.kind, SectionKind::Solid {}) { + return Err(Error::NotSupported(format!( + "{} sections (section `{name}`); explicit dynamics has solid elements so far", + section.kind.name() + ))); + } + let mi = match materials.iter().position(|m| m.name == section.material) { + Some(i) => i, + None => { + materials.push(material(model, §ion.material)?); + materials.len() - 1 + } + }; + let set = resolved.sets.get(§ion.set).ok_or_else(|| invalid(format!("set `{}` is not resolved", section.set)))?; + for &label in &set.elements { + let el = model.element(label).ok_or_else(|| invalid(format!("element {label} is missing")))?; + check_formulation(el.kind, name, §ion.formulation)?; + if el.kind == ElementKind::Tet4 && materials.get(mi).is_some_and(|m| m.bulk.is_none()) { + return Err(Error::NotSupported(format!( + "nodal-pressure tet4 elements of the orthotropic material `{}` (section `{name}`)", + section.material + ))); + } + chosen.push((label, el.kind, el.nodes.clone(), mi)); + } + } + chosen.sort_by_key(|e| e.0); + let mut labels: Vec = chosen.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()))); + } + if labels.is_empty() { + return Err(invalid("the model has no elements to analyse")); + } + let xyz: Vec<[f64; 3]> = + labels.iter().map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| invalid(format!("node {l} is missing")))).collect::>()?; + let n = labels.len(); + let mut sys = System { + elements: chosen.iter().map(|e| e.0).collect(), + materials, + hexes: Vec::new(), + tets: Vec::new(), + mass: vec![0.0; n], + prescribed: vec![None; 3 * n], + walls: Vec::new(), + initial_velocity: vec![0.0; 3 * n], + external: vec![0.0; 3 * n], + history_nodes: Vec::new(), + hourglass: settings.hourglass.kind, + warnings: resolved.warnings.iter().map(ToString::to_string).collect(), + labels, + xyz, + }; + for (index, (label, kind, nodes, mi)) in chosen.into_iter().enumerate() { + let idx: Vec = nodes.iter().map(|l| sys.node(*l).ok_or_else(|| Error::Internal(format!("node {l}")))).collect::>()?; + let mat = sys.materials.get(mi).ok_or_else(|| Error::Internal("material index".into()))?; + let (rho, c) = (mat.density, mat.wave_speed); + let at = |e| Error::element(label, e); + match kind { + ElementKind::Hex8 => { + let nodes: [usize; 8] = idx.try_into().map_err(|_| invalid(format!("element {label} does not have 8 nodes")))?; + let x = nodes.map(|i| sys.xyz.get(i).copied().unwrap_or_default()); + let k = Hex8::hex8(&x).map_err(at)?; + let gamma = k.hourglass_vectors(&x); + // h: the edge of the cube with the element's volume and gradient norm. + let h = k.volume / (2.0 * k.gradient_norm()).sqrt(); + let kappa = settings.hourglass.coefficient; + let mut omega2 = k.frequency_bound(mat.modulus, rho); + let hourglass = match settings.hourglass.kind { + HourglassKind::Viscous => kappa * rho * c * h * h / 16.0, + HourglassKind::Stiffness => { + let stiffness = kappa * rho * c * c * h / 16.0; + omega2 += stiffness * hourglass_gram_bound(&gamma, &k.nodal_volume) / rho; + stiffness + } + }; + sys.lump(&nodes, &k.nodal_volume, rho); + sys.hexes.push(Hex { label, index, nodes, material: mi, kinematics: k, gamma, hourglass, length: 2.0 * c / omega2.sqrt() }); + } + ElementKind::Tet4 => { + let nodes: [usize; 4] = idx.try_into().map_err(|_| invalid(format!("element {label} does not have 4 nodes")))?; + let x = nodes.map(|i| sys.xyz.get(i).copied().unwrap_or_default()); + let k = Tet4::tet4(&x).map_err(at)?; + let omega2 = k.frequency_bound(mat.modulus, rho); + sys.lump(&nodes, &k.nodal_volume, rho); + sys.tets.push(Tet { label, index, nodes, material: mi, kinematics: k, length: 2.0 * c / omega2.sqrt() }); + } + other => return Err(Error::Internal(format!("{} element kind", other.name()))), + } + } + sys.prescribe(model, resolved, step, settings.end_time)?; + sys.load(model, resolved, step)?; + sys.conditions(resolved, settings)?; + Ok(sys) + } + + /// Add an element's row-sum lumped masses. + fn lump(&mut self, nodes: &[usize], nodal_volume: &[f64], density: f64) { + for (node, v) in nodes.iter().zip(nodal_volume) { + if let Some(m) = self.mass.get_mut(*node) { + *m += density * v; + } + } + } + + pub fn node(&self, label: u32) -> Option { + self.labels.binary_search(&label).ok() + } + + pub fn dofs(&self) -> usize { + 3 * self.labels.len() + } + + /// Prescribed velocities from the step's boundary conditions: zero for held degrees of + /// freedom, value / end time for a displacement reached over the step. + fn prescribe(&mut self, model: &Model, resolved: &Resolved, step: &Step, end_time: f64) -> Result<()> { + for name in &step.boundary_conditions { + let bc = model.boundary_condition(name).ok_or_else(|| 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) = self.node(label) else { continue }; + for dof in 1..=3u8 { + let Some(value) = bc.value(dof) else { continue }; + let velocity = value / end_time; + if !velocity.is_finite() { + return Err(invalid(format!("boundary condition `{name}`: the velocity {value} / {end_time} is not finite"))); + } + let slot = self.prescribed.get_mut(3 * node + usize::from(dof - 1)).ok_or_else(|| Error::Internal("dof index".into()))?; + match *slot { + Some(old) if old != velocity => { + return Err(invalid(format!( + "node {label} direction {dof} moves at both {old} and {velocity} (boundary condition `{name}` and another)" + ))); + } + _ => *slot = Some(velocity), + } + } + } + } + Ok(()) + } + + /// The constant load vector: forces and tractions as the model distributes them, consistent + /// pressure loads, gravity and body forces with the lumped (row-sum) weights. + fn load(&mut self, model: &Model, resolved: &Resolved, step: &Step) -> Result<()> { + let mut f = std::mem::take(&mut self.external); + for name in &step.loads { + let load = model.load(name).ok_or_else(|| invalid(format!("load `{name}` is missing")))?; + let off_mesh = |l: u32| invalid(format!("load `{name}`: node {l} belongs to no element")); + match load { + Load::Force { target: Target::ReferencePoint(rp), .. } => { + return Err(Error::NotSupported(format!("loads on reference points (`{name}` on `{rp}`)"))); + } + Load::Moment { .. } => { + return Err(Error::NotSupported(format!("moments (`{name}`): solid nodes have no rotations"))); + } + Load::Force { .. } | Load::Traction { .. } => { + for (label, v) in resolved.nodal_loads.get(name).map(Vec::as_slice).unwrap_or_default() { + let node = self.node(*label).ok_or_else(|| off_mesh(*label))?; + add_at(&mut f, node, v); + } + } + Load::Pressure { set, magnitude } => { + for fc in resolved.sets.get(set).map(|s| s.faces.as_slice()).unwrap_or_default() { + let labels = model + .face_nodes(*fc) + .ok_or_else(|| invalid(format!("load `{name}`: element {} face {} is missing", fc.element, fc.face)))?; + let idx: Vec = labels.iter().map(|l| self.node(*l).ok_or_else(|| off_mesh(*l))).collect::>()?; + let xyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| self.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_at(&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()); + self.each_element(chosen, |nodes, weights, rho| { + for (node, w) in nodes.iter().zip(weights) { + add_at(&mut f, *node, &acceleration.map(|a| rho * w * a)); + } + }); + } + Load::BodyForce { set, vector } => { + let chosen = resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default(); + self.each_element(Some(chosen), |nodes, weights, _| { + for (node, w) in nodes.iter().zip(weights) { + add_at(&mut f, *node, &vector.map(|b| w * b)); + } + }); + } + } + } + self.external = f; + Ok(()) + } + + /// Visit the elements (all, or those listed): their nodes, nodal volumes and density. + fn each_element(&self, only: Option<&[u32]>, mut visit: impl FnMut(&[usize], &[f64], f64)) { + let wanted = |label: u32| only.is_none_or(|l| l.binary_search(&label).is_ok()); + let density = |m: usize| self.materials.get(m).map_or(0.0, |m| m.density); + for h in self.hexes.iter().filter(|h| wanted(h.label)) { + visit(&h.nodes, &h.kinematics.nodal_volume, density(h.material)); + } + for t in self.tets.iter().filter(|t| wanted(t.label)) { + visit(&t.nodes, &t.kinematics.nodal_volume, density(t.material)); + } + } + + /// The analysed nodes of a resolved set. + fn nodes_of(&self, resolved: &Resolved, set: &str) -> Vec { + resolved.sets.get(set).map(|s| s.nodes.as_slice()).unwrap_or_default().iter().filter_map(|l| self.node(*l)).collect() + } + + /// Initial velocities, rigid walls and history nodes of the step. + fn conditions(&mut self, resolved: &Resolved, settings: &ExplicitDynamic) -> Result<()> { + for iv in &settings.initial_velocities { + for node in self.nodes_of(resolved, &iv.set) { + for (k, c) in iv.velocity.iter().enumerate() { + if let Some(v) = self.initial_velocity.get_mut(3 * node + k) { + *v = *c; + } + } + } + } + // A prescribed velocity holds from the start. + for (v, p) in self.initial_velocity.iter_mut().zip(&self.prescribed) { + if let Some(p) = p { + *v = *p; + } + } + let size = self.xyz.iter().flatten().fold(0.0f64, |m, c| m.max(c.abs())).max(1.0); + for w in &settings.rigid_walls { + let len = w.normal.iter().map(|c| c * c).sum::().sqrt(); + let normal = w.normal.map(|c| c / len); + let nodes = self.nodes_of(resolved, &w.set); + for &i in &nodes { + let x = self.xyz.get(i).copied().unwrap_or_default(); + let gap: f64 = (0..3).map(|k| (x[k] - w.point[k]) * normal[k]).sum(); + if gap < -1e-9 * size { + let label = self.labels.get(i).copied().unwrap_or(0); + return Err(invalid(format!("node {label} of set `{}` starts behind its rigid wall (gap {gap:e})", w.set))); + } + } + self.walls.push(Wall { point: w.point, normal, nodes }); + } + if let Some(set) = &settings.history_nodes { + let nodes = self.nodes_of(resolved, set); + if nodes.len() > MAX_HISTORY_NODES { + return Err(Error::Limit(format!("set `{set}` has {} nodes; histories are recorded for at most {MAX_HISTORY_NODES}", nodes.len()))); + } + self.history_nodes = nodes; + } + Ok(()) + } +} diff --git a/crates/explicit/tests/verification.rs b/crates/explicit/tests/verification.rs new file mode 100644 index 0000000..3746a79 --- /dev/null +++ b/crates/explicit/tests/verification.rs @@ -0,0 +1,350 @@ +//! Analytic verification of the explicit solver (docs/solvers/literature.md, A10): elastic waves +//! in a bar (front at c·t, stress σ = ρ c v), a bar striking a rigid wall (contact for 2L/c, +//! rebound at the impact speed), the energy balance and bit-for-bit repeatability. +//! +//! The bar has Poisson's ratio 0, so a uniaxial stress wave travels at the rod speed c = √(E/ρ) +//! without lateral inertia effects (d'Alembert's solution of the one-dimensional wave equation). + +use modelcraft_explicit::{ExplicitResult, Limits, solve_explicit}; +use modelcraft_model::*; + +const YOUNG: f64 = 1.0e4; +const DENSITY: f64 = 1.0; +const LENGTH: f64 = 10.0; +const WIDTH: f64 = 0.1; + +fn wave_speed() -> f64 { + (YOUNG / DENSITY).sqrt() +} + +/// A bar along x from 0 to LENGTH, WIDTH × WIDTH, with `cells` cells along it: one hex8 per +/// cell, or six tet4 (Kuhn's split of the cube), with ν = 0; its sets, and an explicit step. +struct Bar { + model: Model, + kind: ElementKind, + cells: usize, +} + +impl Bar { + fn new(kind: ElementKind, cells: usize, step: ExplicitDynamic) -> Bar { + let mut m = Model::default(); + let h = LENGTH / cells as f64; + let label = |i: usize, j: usize, k: usize| (1 + i + (cells + 1) * (j + 2 * k)) as u32; + for k in 0..2 { + for j in 0..2 { + for i in 0..=cells { + m.set_node(label(i, j, k), [i as f64 * h, j as f64 * WIDTH, k as f64 * WIDTH]).unwrap(); + } + } + } + let corners = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]; + let mut next = 1; + for i in 0..cells { + let at = |c: [usize; 3]| label(i + c[0], c[1], c[2]); + if kind == ElementKind::Hex8 { + m.set_element(next, kind, corners.iter().map(|c| at(*c)).collect()).unwrap(); + next += 1; + continue; + } + 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]] = 1; + } + let x = v.map(|c| [c[0] as f64 * h, c[1] as f64 * WIDTH, c[2] as f64 * WIDTH]); + let e = |a: usize| [0, 1, 2].map(|k| x[a][k] - x[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); + } + m.set_element(next, kind, v.iter().map(|c| at(*c)).collect()).unwrap(); + next += 1; + } + } + let all: Vec = m.elements().keys().copied().collect(); + m.add_set("all", Set { elements: all, ..Set::default() }).unwrap(); + let end_nodes = |i: usize| (0..4).map(|n| label(i, n % 2, n / 2)).collect::>(); + m.add_set("start", Set { nodes: end_nodes(0), ..Set::default() }).unwrap(); + m.add_set("end", Set { nodes: end_nodes(cells), ..Set::default() }).unwrap(); + let mut faces = Vec::new(); + for (l, e) in m.elements() { + for (f, local) in e.kind.faces().iter().enumerate() { + if local.iter().all(|n| m.node(e.nodes[*n]).unwrap().xyz[0] == LENGTH) { + faces.push(ElementFace { element: *l, face: f as u8 + 1 }); + } + } + } + m.add_set("end_face", Set { faces, ..Set::default() }).unwrap(); + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: Some(DENSITY) }).unwrap(); + let integration = if kind == ElementKind::Hex8 { Integration::Reduced } else { Integration::Full }; + let formulation = Formulation { integration, ..Formulation::default() }; + m.add_section("bar", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap(); + m.add_step(Step::new("Dynamic", Procedure::ExplicitDynamic(step)), None).unwrap(); + Bar { model: m, kind, cells } + } + + fn per_cell(&self) -> usize { + if self.kind == ElementKind::Hex8 { 1 } else { 6 } + } + + /// σ_xx per cell (the mean over its elements) and the cell's centre. + fn axial_stress(&self, r: &ExplicitResult, frame: usize) -> Vec<(f64, f64)> { + let f = &r.frames[frame]; + let h = LENGTH / self.cells as f64; + (0..self.cells) + .map(|i| { + let s: f64 = f.stress[i * self.per_cell()..(i + 1) * self.per_cell()].iter().map(|s| s[0]).sum(); + ((i as f64 + 0.5) * h, s / self.per_cell() as f64) + }) + .collect() + } + + /// The axial velocity of the nodes on the bar's axis line y = z = 0, by x. + fn axial_velocity(&self, r: &ExplicitResult, frame: usize) -> Vec<(f64, f64)> { + let f = &r.frames[frame]; + let h = LENGTH / self.cells as f64; + (0..=self.cells).map(|i| (i as f64 * h, f.velocity[r.nodes.binary_search(&(1 + i as u32)).unwrap()][0])).collect() + } + + fn activate(&mut self, bc: Option<(&str, BoundaryCondition)>, load: Option<(&str, Load)>) { + if let Some((name, bc)) = bc { + self.model.add_boundary_condition(name, bc).unwrap(); + self.model.set_boundary_condition_active("Dynamic", name, true).unwrap(); + } + if let Some((name, load)) = load { + self.model.add_load(name, load).unwrap(); + self.model.set_load_active("Dynamic", name, true).unwrap(); + } + } +} + +/// Where a profile first falls below half of `level` going along x (linear interpolation). +fn front(profile: &[(f64, f64)], level: f64) -> f64 { + for w in profile.windows(2) { + let (a, b) = (w[0].1 / level, w[1].1 / level); + if a >= 0.5 && b < 0.5 { + return w[0].0 + (a - 0.5) / (a - b) * (w[1].0 - w[0].0); + } + } + f64::NAN +} + +fn mean_between(profile: &[(f64, f64)], from: f64, to: f64) -> f64 { + let inside: Vec = profile.iter().filter(|(x, _)| *x >= from && *x <= to).map(|(_, v)| *v).collect(); + inside.iter().sum::() / inside.len() as f64 +} + +/// A pressure step on the free end of a bar sends a compressive wave in −x: front at L − c t, +/// stress −p, particle velocity −p/(ρc), so σ = ρ c v. +#[test] +fn pressure_wave_in_a_hex_bar() { + let c = wave_speed(); + let t = 0.6 * LENGTH / c; + let mut step = ExplicitDynamic::new(t); + step.output_interval = Some(t / 4.0); + let mut bar = Bar::new(ElementKind::Hex8, 200, step); + let p = 1.0; + bar.activate(None, Some(("push", Load::Pressure { set: "end_face".into(), magnitude: p }))); + let r = solve_explicit(&bar.model, None, &Limits::default()).unwrap(); + assert!(r.completed, "{:?}", r.warnings); + let last = r.frames.len() - 1; + assert!((r.frames[last].time - t).abs() < 1e-12 * t); + // Measured from the loaded end, the profiles are those of a wave travelling in +x. + let flip = |v: Vec<(f64, f64)>| -> Vec<(f64, f64)> { v.into_iter().rev().map(|(x, s)| (LENGTH - x, s)).collect() }; + let stress = flip(bar.axial_stress(&r, last)); + let velocity = flip(bar.axial_velocity(&r, last)); + let h = LENGTH / 200.0; + let x_front = front(&stress, -p); + let x_vfront = front(&velocity, -p / (DENSITY * c)); + let sigma = mean_between(&stress, 0.1 * c * t, 0.8 * c * t); + let v = mean_between(&velocity, 0.1 * c * t, 0.8 * c * t); + eprintln!( + "hex8 pressure wave: front {x_front:.4} and {x_vfront:.4} vs c t = {:.4} (h = {h}); σ {sigma:.5} vs −p; v {v:.6} vs {:.6}; ρ c v = {:.5}; energy error {:.2e}; Δt {:.3e}, {} increments", + c * t, + -p / (DENSITY * c), + DENSITY * c * v, + r.energy.relative_error(r.initial_kinetic), + r.initial_increment, + r.increments + ); + assert!((x_front - c * t).abs() < 2.0 * h, "stress front {x_front} vs {}", c * t); + assert!((x_vfront - c * t).abs() < 2.0 * h, "velocity front {x_vfront} vs {}", c * t); + assert!((sigma + p).abs() < 0.01 * p, "σ = {sigma}"); + assert!((DENSITY * c * v - sigma).abs() < 0.01 * p, "ρ c v = {} vs σ = {sigma}", DENSITY * c * v); + // Work of the pressure = kinetic + strain energy (+ the viscous dissipation at the front). + assert!(r.energy.external > 0.0 && r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy); + assert!(r.warnings.iter().all(|w| w.contains("unit system")), "{:?}", r.warnings); + // The exact energy input is p A (c t)(p / E) = p² A c t / E; half is kinetic, half strain. + let input = p * p * WIDTH * WIDTH * c * t / YOUNG; + eprintln!("hex8 pressure wave: work {:.6e} vs p² A c t / E = {input:.6e}", r.energy.external); + assert!((r.energy.external - input).abs() < 0.005 * input, "{} vs {input}", r.energy.external); +} + +/// A velocity v prescribed at one end of a tetrahedral bar (nodal-pressure tet4): front at c t, +/// stress −ρ c v behind it. +#[test] +fn velocity_wave_in_a_tet_bar() { + let c = wave_speed(); + let t = 0.6 * LENGTH / c; + let v = 0.01; + let mut step = ExplicitDynamic::new(t); + step.output_interval = Some(t / 2.0); + let mut bar = Bar::new(ElementKind::Tet4, 200, step); + let push = BoundaryCondition::Displacement { target: Target::Set("start".into()), values: [Some(v * t), None, None, None, None, None] }; + bar.activate(Some(("push", push)), None); + let r = solve_explicit(&bar.model, None, &Limits::default()).unwrap(); + let last = r.frames.len() - 1; + let stress = bar.axial_stress(&r, last); + let h = LENGTH / 200.0; + let sigma_exact = -DENSITY * c * v; + let x_front = front(&stress, sigma_exact); + let sigma = mean_between(&stress, 0.1 * c * t, 0.8 * c * t); + eprintln!( + "tet4 velocity wave: front {x_front:.4} vs c t = {:.4} (h = {h}); σ {sigma:.6} vs −ρcv = {sigma_exact:.6}; energy error {:.2e}; Δt {:.3e}", + c * t, + r.energy.relative_error(r.initial_kinetic), + r.initial_increment + ); + assert!((x_front - c * t).abs() < 3.0 * h, "front {x_front} vs {}", c * t); + assert!((sigma - sigma_exact).abs() < 0.02 * sigma_exact.abs(), "σ = {sigma}"); + assert!(r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy); + // The boundary condition does the work: ρ c v² A t (kinetic and strain energy). + let input = DENSITY * c * v * v * WIDTH * WIDTH * t; + eprintln!("tet4 velocity wave: work {:.6e} vs ρ c v² A t = {input:.6e}", r.energy.constraint); + assert!((r.energy.constraint - input).abs() < 0.01 * input, "{} vs {input}", r.energy.constraint); +} + +/// A bar moving at v₀ strikes a rigid wall: it stays in contact for 2L/c, then leaves at v₀ in +/// the opposite direction with its kinetic energy (elastic impact). The comparison is with the +/// elastic solution, so the linear bulk viscosity (which dissipates a few percent of the energy +/// at the wave fronts) is off; the kinematic wall stops the first layer of nodes dead, an O(h) +/// loss that the constraint work accounts for. +#[test] +fn bar_rebounds_from_a_rigid_wall() { + let c = wave_speed(); + let v0 = 0.01; + let end = 3.0 * LENGTH / c; + for (kind, linear) in [(ElementKind::Hex8, 0.0), (ElementKind::Tet4, 0.0), (ElementKind::Hex8, 0.06)] { + let mut step = ExplicitDynamic::new(end); + step.history_interval = Some(end / 3000.0); + step.output_interval = Some(end); + step.initial_velocities = vec![InitialVelocity { set: "all_nodes".into(), velocity: [-v0, 0.0, 0.0] }]; + step.rigid_walls = vec![RigidWall { set: "start".into(), point: [0.0; 3], normal: [1.0, 0.0, 0.0] }]; + step.bulk_viscosity.linear = linear; + let mut bar = Bar::new(kind, 200, ExplicitDynamic::new(end)); + let nodes: Vec = bar.model.nodes().keys().copied().collect(); + bar.model.add_set("all_nodes", Set { nodes, ..Set::default() }).unwrap(); + bar.model.replace_step(Step::new("Dynamic", Procedure::ExplicitDynamic(step))).unwrap(); + let r = solve_explicit(&bar.model, Some("Dynamic"), &Limits::default()).unwrap(); + let hist = &r.history; + let peak = DENSITY * c * v0 * WIDTH * WIDTH; + let force: Vec = hist.constraint_force.iter().map(|f| f[0]).collect(); + // Contact lasts until the wall force drops for good. + let release = hist.time.iter().zip(&force).filter(|(_, f)| **f > 0.01 * peak).map(|(t, _)| *t).fold(0.0, f64::max); + let rebound = hist.momentum.last().unwrap()[0] / r.total_mass; + let error = r.energy.relative_error(r.initial_kinetic); + let kinetic = r.energy.kinetic / r.initial_kinetic; + let mid = hist.time.iter().position(|t| *t >= LENGTH / c).unwrap(); + eprintln!( + "{} impact (linear q {linear}): contact {release:.5} vs 2L/c = {:.5}; wall force at L/c {:.4e} vs ρcv₀A = {peak:.4e}; rebound {rebound:.6} vs {v0}; final KE/KE₀ {kinetic:.4}; energy balance {error:.2e}; {:?}", + kind.name(), + 2.0 * LENGTH / c, + force[mid], + r.energy + ); + assert!(r.warnings.iter().all(|w| w.contains("unit system")), "{:?}", r.warnings); + assert!(error < 0.01, "{} {:?}", kind.name(), r.energy); + // The first layer of nodes (half an element's mass) stops dead on the wall. + let stop = 0.5 * (0.5 * DENSITY * WIDTH * WIDTH * LENGTH / 200.0) * v0 * v0; + assert!((r.energy.constraint + stop).abs() < 1e-6 * stop, "{} {:?}", kind.name(), r.energy); + if linear > 0.0 { + assert!(r.energy.viscous > 0.01 * r.initial_kinetic, "{:?}", r.energy); + continue; + } + assert!((release - 2.0 * LENGTH / c).abs() < 0.01 * 2.0 * LENGTH / c, "{} contact {release}", kind.name()); + assert!((rebound - v0).abs() < 0.01 * v0, "{} rebound {rebound}", kind.name()); + assert!(kinetic > 0.99, "{} kinetic {kinetic}", kind.name()); + // While in contact the plateau force is ρ c v₀ A. + assert!((force[mid] - peak).abs() < 0.03 * peak, "{} force {} vs {peak}", kind.name(), force[mid]); + } +} + +/// The same model gives the same numbers, bit for bit. +#[test] +fn runs_are_bit_identical() { + let end = 0.3 * LENGTH / wave_speed(); + let mut step = ExplicitDynamic::new(end); + step.hourglass = Hourglass { kind: HourglassKind::Stiffness, coefficient: 0.1 }; + let mut bar = Bar::new(ElementKind::Hex8, 40, step); + bar.activate( + None, + Some(("push", Load::Force { target: Target::Set("end".into()), vector: [0.0, 0.02, -0.01], distribution: Distribution::Total })), + ); + let a = solve_explicit(&bar.model, None, &Limits::default()).unwrap(); + let b = solve_explicit(&bar.model, None, &Limits::default()).unwrap(); + assert_eq!(a, b); + assert!(a.increments > 10); +} + +/// Hostile and unsupported models are errors, never panics; the caps hold. +#[test] +fn caps_and_unsupported_models_are_errors() { + let end = 0.5 * LENGTH / wave_speed(); + let bar = Bar::new(ElementKind::Hex8, 20, ExplicitDynamic::new(end)); + let tight = |f: fn(&mut Limits)| { + let mut l = Limits::default(); + f(&mut l); + l + }; + let e = solve_explicit(&bar.model, None, &tight(|l| l.max_increments = 3)).unwrap_err().to_string(); + assert!(e.contains("increments") && e.contains("element"), "{e}"); + assert!(solve_explicit(&bar.model, None, &tight(|l| l.max_dofs = 10)).is_err()); + assert!(solve_explicit(&bar.model, None, &tight(|l| l.max_output_values = 10)).is_err()); + assert!(solve_explicit(&bar.model, Some("Nope"), &Limits::default()).is_err()); + + // A step capped below what it needs stops early, says so and keeps what it computed. + let mut capped = ExplicitDynamic::new(end); + capped.max_increments = 1_000_000; + let mut short = Bar::new(ElementKind::Hex8, 20, capped); + let r = solve_explicit(&short.model, None, &Limits::default()).unwrap(); + assert!(r.completed && r.increments > 5); + let needed = r.increments; + let Procedure::ExplicitDynamic(mut s) = short.model.step("Dynamic").unwrap().procedure.clone() else { unreachable!() }; + s.max_increments = needed; + s.dt_scale = 0.89; + short.model.replace_step(Step::new("Dynamic", Procedure::ExplicitDynamic(s))).unwrap(); + let r = solve_explicit(&short.model, None, &Limits::default()); + // Either refused up front (the estimate exceeds the cap) or stopped at the cap with a warning. + match r { + Ok(r) => assert!(!r.completed && r.warnings.iter().any(|w| w.contains("cap")) && r.time < end, "{r:?}"), + Err(e) => assert!(e.to_string().contains("cap"), "{e}"), + } + + // Fully integrated hex8, orthotropic tets and missing densities are refused by name. + let mut full = Bar::new(ElementKind::Hex8, 4, ExplicitDynamic::new(end)); + full.model + .replace_section( + "bar", + Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation: Formulation::default() }, + ) + .unwrap(); + assert!(solve_explicit(&full.model, None, &Limits::default()).unwrap_err().to_string().contains("not supported yet")); + let mut ortho = Bar::new(ElementKind::Tet4, 4, ExplicitDynamic::new(end)); + let constants = Elastic::EngineeringConstants { young: [YOUNG; 3], poisson: [0.0; 3], shear: [YOUNG / 2.0; 3] }; + ortho.model.replace_material("Mat", Material { elastic: Some(constants), density: Some(DENSITY) }).unwrap(); + assert!(solve_explicit(&ortho.model, None, &Limits::default()).unwrap_err().to_string().contains("orthotropic")); + let mut light = Bar::new(ElementKind::Tet4, 4, ExplicitDynamic::new(end)); + light.model.replace_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: None }).unwrap(); + assert!(solve_explicit(&light.model, None, &Limits::default()).is_err()); + // A node behind its wall at the start. + let mut behind = ExplicitDynamic::new(end); + behind.rigid_walls = vec![RigidWall { set: "end".into(), point: [LENGTH + 1.0, 0.0, 0.0], normal: [-1.0, 0.0, 0.0] }]; + let mut wall = Bar::new(ElementKind::Hex8, 4, ExplicitDynamic::new(end)); + wall.model.replace_step(Step::new("Dynamic", Procedure::ExplicitDynamic(behind.clone()))).unwrap(); + assert!(solve_explicit(&wall.model, None, &Limits::default()).is_ok(), "the wall is beyond the end node, which is on its allowed side"); + behind.rigid_walls[0].point = [LENGTH - 1.0, 0.0, 0.0]; + wall.model.replace_step(Step::new("Dynamic", Procedure::ExplicitDynamic(behind))).unwrap(); + assert!(solve_explicit(&wall.model, None, &Limits::default()).unwrap_err().to_string().contains("behind")); +} diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index 672fa00..750b57d 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -9,6 +9,8 @@ //! - [`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; +//! - [`uniform`]: uniform-strain (one-point) hexahedron and tetrahedron kinematics, hourglass +//! shape vectors and element frequency bounds for explicit dynamics; //! - [`dense`]: the small dense linear algebra inside an element. //! //! Conventions: three translational degrees of freedom per node, ordered node by node; stress @@ -24,6 +26,7 @@ pub mod face; pub mod material; pub mod quadrature; pub mod shape; +pub mod uniform; pub use element::{ELEMENT_TYPES, ElementType, Family, Formulation, PointStress, Solid, element_type}; pub use material::{Elasticity, elasticity, von_mises}; diff --git a/crates/fem/src/uniform.rs b/crates/fem/src/uniform.rs new file mode 100644 index 0000000..a577fb1 --- /dev/null +++ b/crates/fem/src/uniform.rs @@ -0,0 +1,358 @@ +//! Uniform-strain (one-point) kinematics of the eight-node hexahedron and the four-node +//! tetrahedron, for explicit dynamics. +//! +//! **Hexahedron** (D. P. Flanagan, T. Belytschko, *A uniform strain hexahedron and quadrilateral +//! with orthogonal hourglass control*, IJNME 17(5):679–706, 1981): the strain is the element's +//! volume average, computed with the integrated gradients b_iI = ∫ ∂N_I/∂x_i dV, so +//! ε̄ = (1/V) Σ_I sym(u_I ⊗ b_I) and the nodal forces of a stress are f_iI = σ_ij b_jI. The +//! mean strain is exact for every linear displacement field on any hexahedron (the element +//! passes the patch test). The products det J · ∂N/∂x are polynomials of degree two at most in +//! each natural coordinate, so the 2 × 2 × 2 Gauss rule integrates b and V exactly. +//! +//! The one-point element has twelve hourglass modes that the mean strain does not see: in each +//! direction, the four base vectors h₁ = ηζ, h₂ = ξζ, h₃ = ξη, h₄ = ξηζ at the corners. The +//! hourglass shape vectors γ_α = h_α − (1/V) Σ_J (h_αJ x_iJ) b_i (Flanagan & Belytschko 1981) +//! are orthogonal to every linear velocity field, so forces resisting γ_α leave the constant +//! strain modes alone and the patch test passes on distorted elements too. +//! +//! **Tetrahedron**: the linear tetrahedron's strain is uniform, b_I = V ∇N_I. +//! +//! **Critical frequency.** The element eigenvalue inequality (T. Belytschko, W. K. Liu, B. Moran, +//! K. Elkhodary, *Nonlinear Finite Elements for Continua and Structures*, 2nd ed., Wiley 2014, +//! §6.3) bounds the mesh's highest frequency by the largest element frequency. For the uniform +//! strain stiffness K = (1/V) B̃ᵀ D B̃ we bound it in closed form, in the spirit of D. P. +//! Flanagan, T. Belytschko, *Eigenvalues and stable time steps for the uniform strain hexahedron +//! and quadrilateral*, J. Appl. Mech. 51(1):35–40, 1984: with G = Σ_I u_I ⊗ b_I, +//! uᵀ K u ≤ (k/V) ‖G‖² ≤ (k/V) λ_max(Σ_I b_I b_Iᵀ) ‖u‖², where k bounds the material's +//! energy per squared displacement gradient ([`modulus_bound`]); with the lumped masses m_I the +//! element frequency satisfies ω² ≤ k λ_max(Σ b bᵀ) / (V min m_I). The bound is exact for a +//! cube (its dilatation mode) and conservative for distorted elements. + +use modelcraft_model::Elastic; + +use crate::dense::inverse3; +use crate::material::Elasticity; +use crate::quadrature::Rule; +use crate::shape::{self, HEX_CORNERS}; +use crate::{Error, Result}; + +/// The integrated gradients, volume and nodal volumes of a uniform-strain element with `N` nodes. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct UniformStrain { + /// b_I = ∫ ∇N_I dV. + pub b: [[f64; 3]; N], + pub volume: f64, + /// ∫ N_I dV: the row sums of the consistent mass per unit density (row-sum lumping). + pub nodal_volume: [f64; N], +} + +/// The one-point hexahedron. +pub type Hex8 = UniformStrain<8>; +/// The linear tetrahedron. +pub type Tet4 = UniformStrain<4>; + +fn degenerate(what: &str) -> Error { + Error::Degenerate(format!("a {what} element is inverted or degenerate (its Jacobian is not positive)")) +} + +impl UniformStrain<8> { + /// The uniform-strain hexahedron of corners `xyz` (the model's hex8 node order). + pub fn hex8(xyz: &[[f64; 3]; 8]) -> Result { + modelcraft_model::check_coordinates(xyz.as_flattened()).map_err(|e| Error::Invalid(e.to_string()))?; + let mut b = [[0.0; 3]; 8]; + let mut nodal_volume = [0.0; 8]; + let mut volume = 0.0; + for p in Rule::Hex(2).points() { + let s = shape::hex8(p.xi); + let mut j = [[0.0; 3]; 3]; + for (d, x) in s.derivatives().iter().zip(xyz) { + for a in 0..3 { + for c in 0..3 { + j[a][c] += d[a] * x[c]; + } + } + } + let (inv, det) = inverse3(&j).ok_or_else(|| degenerate("hex8"))?; + 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(degenerate("hex8")); + } + let w = det * p.weight; + volume += w; + for (((bi, vi), d), n) in b.iter_mut().zip(nodal_volume.iter_mut()).zip(s.derivatives()).zip(s.values()) { + for c in 0..3 { + bi[c] += w * (inv[c][0] * d[0] + inv[c][1] * d[1] + inv[c][2] * d[2]); + } + *vi += w * n; + } + } + Ok(UniformStrain { b, volume, nodal_volume }) + } + + /// The hourglass shape vectors γ_α (α = 1..4) of the element with corners `xyz`. + pub fn hourglass_vectors(&self, xyz: &[[f64; 3]; 8]) -> [[f64; 8]; 4] { + let mut gamma = [[0.0; 8]; 4]; + for (alpha, g) in gamma.iter_mut().enumerate() { + let h: [f64; 8] = HEX_CORNERS.map(|c| match alpha { + 0 => c[1] * c[2], + 1 => c[0] * c[2], + 2 => c[0] * c[1], + _ => c[0] * c[1] * c[2], + }); + // (h · x_i) for each coordinate direction i. + let hx = [0, 1, 2].map(|i| h.iter().zip(xyz).map(|(hj, x)| hj * x[i]).sum::()); + for (k, gk) in g.iter_mut().enumerate() { + *gk = h[k] - (hx[0] * self.b[k][0] + hx[1] * self.b[k][1] + hx[2] * self.b[k][2]) / self.volume; + } + } + gamma + } +} + +impl UniformStrain<4> { + /// The linear tetrahedron of corners `xyz` (the model's tet4 node order). + pub fn tet4(xyz: &[[f64; 3]; 4]) -> Result { + modelcraft_model::check_coordinates(xyz.as_flattened()).map_err(|e| Error::Invalid(e.to_string()))?; + let [x0, x1, x2, x3] = *xyz; + let e = |x: [f64; 3]| [x[0] - x0[0], x[1] - x0[1], x[2] - x0[2]]; + // J[a][c] = ∂x_c/∂ξ_a: the edges from corner 1. + let j = [e(x1), e(x2), e(x3)]; + let (inv, det) = inverse3(&j).ok_or_else(|| degenerate("tet4"))?; + 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(degenerate("tet4")); + } + let volume = det / 6.0; + let s = shape::tet4([0.25; 3]); + let mut b = [[0.0; 3]; 4]; + for (bi, d) in b.iter_mut().zip(s.derivatives()) { + for c in 0..3 { + bi[c] = volume * (inv[c][0] * d[0] + inv[c][1] * d[1] + inv[c][2] * d[2]); + } + } + Ok(UniformStrain { b, volume, nodal_volume: [volume / 4.0; 4] }) + } +} + +impl UniformStrain { + /// The displacement gradient (1/V) Σ_I u_I ⊗ b_I: `[i][j]` = ∂u_i/∂x_j. + pub fn gradient(&self, u: &[[f64; 3]; N]) -> [[f64; 3]; 3] { + let mut g = [[0.0; 3]; 3]; + for (ui, bi) in u.iter().zip(&self.b) { + for i in 0..3 { + for j in 0..3 { + g[i][j] += ui[i] * bi[j]; + } + } + } + g.map(|r| r.map(|v| v / self.volume)) + } + + /// The mean strain of nodal displacements (or its rate, of nodal velocities): xx, yy, zz, + /// xy, xz, yz with engineering shears. + pub fn strain(&self, u: &[[f64; 3]; N]) -> [f64; 6] { + let g = self.gradient(u); + [g[0][0], g[1][1], g[2][2], g[0][1] + g[1][0], g[0][2] + g[2][0], g[1][2] + g[2][1]] + } + + /// The trace of the mean gradient (1/V) Σ_I u_I · b_I: the volumetric strain of nodal + /// displacements, or its rate of nodal velocities. + pub fn divergence(&self, u: &[[f64; 3]; N]) -> f64 { + let mut d = [0.0; 3]; + for (ui, bi) in u.iter().zip(&self.b) { + for i in 0..3 { + d[i] += ui[i] * bi[i]; + } + } + (d[0] + d[1] + d[2]) / self.volume + } + + /// The nodal forces of a uniform pressure p (positive in compression): f_I = −p b_I. + pub fn pressure_forces(&self, p: f64) -> [[f64; 3]; N] { + self.b.map(|bi| bi.map(|c| -p * c)) + } + + /// The nodal forces that balance a uniform stress (xx, yy, zz, xy, xz, yz): f_iI = σ_ij b_jI, + /// the internal force ∫ Bᵀ σ dV. + pub fn forces(&self, s: &[f64; 6]) -> [[f64; 3]; N] { + let sigma = [[s[0], s[3], s[4]], [s[3], s[1], s[5]], [s[4], s[5], s[2]]]; + self.b.map(|bi| [0, 1, 2].map(|i| sigma[i][0] * bi[0] + sigma[i][1] * bi[1] + sigma[i][2] * bi[2])) + } + + /// The largest eigenvalue of Σ_I b_I b_Iᵀ. + pub fn gradient_norm(&self) -> f64 { + let mut m = [[0.0; 3]; 3]; + for bi in &self.b { + for i in 0..3 { + for j in 0..3 { + m[i][j] += bi[i] * bi[j]; + } + } + } + largest_eigenvalue(&m) + } + + /// An upper bound of the element's highest natural frequency squared (ω²) with the + /// row-sum lumped mass, for a material whose [`modulus_bound`] is `modulus`. + pub fn frequency_bound(&self, modulus: f64, density: f64) -> f64 { + let lightest = self.nodal_volume.iter().copied().fold(f64::INFINITY, f64::min); + modulus * self.gradient_norm() / (self.volume * density * lightest) + } +} + +/// The largest eigenvalue of a symmetric 3 × 3 matrix, in closed form (O. K. Smith, *Eigenvalues +/// of a symmetric 3 × 3 matrix*, Commun. ACM 4(4):168, 1961). +pub fn largest_eigenvalue(a: &[[f64; 3]; 3]) -> f64 { + let off = a[0][1] * a[0][1] + a[0][2] * a[0][2] + a[1][2] * a[1][2]; + let q = (a[0][0] + a[1][1] + a[2][2]) / 3.0; + let p2 = (a[0][0] - q).powi(2) + (a[1][1] - q).powi(2) + (a[2][2] - q).powi(2) + 2.0 * off; + if !(p2 > 1e-30 * q * q) { + return a[0][0].max(a[1][1]).max(a[2][2]); + } + let p = (p2 / 6.0).sqrt(); + let mut b = *a; + for (i, row) in b.iter_mut().enumerate() { + for (j, v) in row.iter_mut().enumerate() { + *v = (*v - if i == j { q } else { 0.0 }) / p; + } + } + let r = (crate::dense::det3(&b) / 2.0).clamp(-1.0, 1.0); + q + 2.0 * p * (r.acos() / 3.0).cos() +} + +/// A bound k of a material's strain energy density per squared displacement gradient: +/// εᵀ D ε ≤ k ‖∇u‖² for every displacement gradient. Isotropic: λ(tr ∇u)² + 2μ‖sym ∇u‖² ≤ +/// max(3λ + 2μ, 2μ)‖∇u‖². Otherwise the largest eigenvalue of D (bounded by its largest +/// absolute row sum) times 2, because the engineering shear strains satisfy γ² ≤ 2(g_ij² + g_ji²). +pub fn modulus_bound(elastic: &Elastic, d: &Elasticity) -> f64 { + match *elastic { + Elastic::Isotropic { young, poisson } => { + let lambda = young * poisson / ((1.0 + poisson) * (1.0 - 2.0 * poisson)); + let mu = young / (2.0 * (1.0 + poisson)); + (3.0 * lambda + 2.0 * mu).max(2.0 * mu) + } + Elastic::EngineeringConstants { .. } => 2.0 * d.iter().map(|row| row.iter().map(|v| v.abs()).sum::()).fold(0.0, f64::max), + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::material::isotropic; + + /// A distorted hexahedron (every face warped) and a tetrahedron. + fn distorted() -> [[f64; 3]; 8] { + [[0.0, 0.0, 0.0], [2.1, 0.2, -0.1], [2.3, 1.9, 0.2], [-0.2, 1.7, 0.1], [0.1, -0.1, 1.8], [2.0, 0.1, 2.2], [2.4, 2.2, 1.9], [0.2, 1.8, 2.1]] + } + + const TET: [[f64; 3]; 4] = [[0.1, 0.0, 0.2], [1.3, 0.1, 0.0], [0.2, 1.1, 0.1], [0.3, 0.2, 0.9]]; + + #[test] + fn mean_strain_is_exact_for_linear_fields() { + let a = [[0.003, -0.001, 0.002], [0.004, -0.002, 0.0005], [-0.0015, 0.001, 0.0025]]; + let field = |x: &[f64; 3]| [0, 1, 2].map(|i| 0.7 + a[i][0] * x[0] + a[i][1] * x[1] + a[i][2] * x[2]); + let expected = [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 hex = distorted(); + let h = Hex8::hex8(&hex).unwrap(); + let t = Tet4::tet4(&TET).unwrap(); + for (got, b, x) in + [(h.strain(&hex.map(|x| field(&x))), h.b.as_slice(), hex.as_slice()), (t.strain(&TET.map(|x| field(&x))), t.b.as_slice(), TET.as_slice())] + { + assert!(got.iter().zip(expected).all(|(g, e)| (g - e).abs() < 1e-15), "{got:?}"); + // Σ_I b_iI = 0 and Σ_I b_iI x_jI = V δ_ij. + let v: f64 = (0..3).map(|i| b.iter().zip(x).map(|(bi, xi)| bi[i] * xi[i]).sum::()).sum::() / 3.0; + for i in 0..3 { + assert!(b.iter().map(|bi| bi[i]).sum::().abs() < 1e-14); + for j in 0..3 { + let s: f64 = b.iter().zip(x).map(|(bi, xi)| bi[i] * xi[j]).sum(); + assert!((s - if i == j { v } else { 0.0 }).abs() < 1e-14, "{i} {j}: {s}"); + } + } + } + assert!((h.nodal_volume.iter().sum::() - h.volume).abs() < 1e-14); + // The volume of the distorted hexahedron agrees with the fem element's. + let solid = crate::Solid::new(modelcraft_model::ElementKind::Hex8, &hex, crate::Formulation::Full).unwrap(); + assert!((solid.volume().unwrap() - h.volume).abs() < 1e-14); + let mut inverted = hex; + inverted.swap(0, 1); + inverted.swap(4, 5); + assert!(Hex8::hex8(&inverted).is_err()); + assert!(Tet4::tet4(&[TET[1], TET[0], TET[2], TET[3]]).is_err()); + assert!(Hex8::hex8(&[[f64::NAN; 3]; 8]).is_err()); + } + + #[test] + fn hourglass_vectors_are_orthogonal_to_linear_fields() { + let hex = distorted(); + let h = Hex8::hex8(&hex).unwrap(); + for g in h.hourglass_vectors(&hex) { + assert!(g.iter().sum::().abs() < 1e-14); + for i in 0..3 { + let s: f64 = g.iter().zip(&hex).map(|(gk, x)| gk * x[i]).sum(); + assert!(s.abs() < 1e-13, "{s}"); + } + } + // On a parallelepiped the shape vectors are the base vectors themselves. + let cube: [[f64; 3]; 8] = HEX_CORNERS.map(|c| [c[0] + 0.3 * c[1], c[1], c[2] + 0.2 * c[0]]); + let hc = Hex8::hex8(&cube).unwrap(); + let g = hc.hourglass_vectors(&cube); + assert!(g[3].iter().zip(HEX_CORNERS).all(|(v, c)| (v - c[0] * c[1] * c[2]).abs() < 1e-14), "{g:?}"); + } + + /// The largest eigenvalue of M⁻¹K by power iteration on M^{-1/2} K M^{-1/2}. + fn largest_frequency(e: &UniformStrain, d: &Elasticity, density: f64) -> f64 { + let n = 3 * N; + let mut k = vec![0.0; n * n]; + for a in 0..n { + let mut u = [[0.0; 3]; N]; + u[a / 3][a % 3] = 1.0; + let f = e.forces(&crate::material::stress(d, &e.strain(&u))); + for b in 0..n { + k[b * n + a] = f[b / 3][b % 3]; + } + } + let m: Vec = (0..n).map(|i| density * e.nodal_volume[i / 3]).collect(); + let mut x: Vec = (0..n).map(|i| 1.0 + (i as f64 * 0.37).sin()).collect(); + let mut lambda = 0.0; + for _ in 0..20_000 { + let y: Vec = (0..n).map(|i| (0..n).map(|j| k[i * n + j] * x[j] / m[j].sqrt()).sum::() / m[i].sqrt()).collect(); + let norm = y.iter().map(|v| v * v).sum::().sqrt(); + lambda = x.iter().zip(&y).map(|(a, b)| a * b).sum::() / x.iter().map(|v| v * v).sum::(); + x = y.iter().map(|v| v / norm).collect(); + } + lambda + } + + #[test] + fn frequency_bound_holds_and_is_exact_for_a_cube() { + let cube: [[f64; 3]; 8] = HEX_CORNERS.map(|c| c.map(|v| 0.5 * (v + 1.0))); + for nu in [0.0, 0.3, 0.45] { + let elastic = Elastic::Isotropic { young: 1000.0, poisson: nu }; + let d = isotropic(1000.0, nu).unwrap(); + let k = modulus_bound(&elastic, &d); + let h = Hex8::hex8(&cube).unwrap(); + let (bound, actual) = (h.frequency_bound(k, 2.0), largest_frequency(&h, &d, 2.0)); + assert!((actual - bound).abs() < 1e-6 * bound, "ν = {nu}: {actual} vs {bound}"); + let hex = distorted(); + let h = Hex8::hex8(&hex).unwrap(); + let t = Tet4::tet4(&TET).unwrap(); + for (bound, actual) in + [(h.frequency_bound(k, 2.0), largest_frequency(&h, &d, 2.0)), (t.frequency_bound(k, 2.0), largest_frequency(&t, &d, 2.0))] + { + assert!(actual <= bound * (1.0 + 1e-9) && actual > 0.3 * bound, "ν = {nu}: {actual} vs {bound}"); + } + } + let ortho = Elastic::EngineeringConstants { young: [140e3, 10e3, 12e3], poisson: [0.3, 0.25, 0.4], shear: [5e3, 4e3, 3.8e3] }; + let d = crate::material::elasticity(&ortho).unwrap(); + let h = Hex8::hex8(&distorted()).unwrap(); + assert!(largest_frequency(&h, &d, 1.0) <= h.frequency_bound(modulus_bound(&ortho, &d), 1.0)); + } + + #[test] + fn closed_form_eigenvalue() { + assert_eq!(largest_eigenvalue(&[[2.0, 0.0, 0.0], [0.0, 5.0, 0.0], [0.0, 0.0, 1.0]]), 5.0); + // [[2, 1, 0], [1, 2, 0], [0, 0, 1]] has eigenvalues 3, 1, 1. + assert!((largest_eigenvalue(&[[2.0, 1.0, 0.0], [1.0, 2.0, 0.0], [0.0, 0.0, 1.0]]) - 3.0).abs() < 1e-14); + assert!((largest_eigenvalue(&[[4.0, 1.0, 2.0], [1.0, 3.0, 0.5], [2.0, 0.5, 5.0]]) - 6.831_254_430_698_204).abs() < 1e-9); + } +} diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs index aedd06a..b1a7797 100644 --- a/crates/model/src/check.rs +++ b/crates/model/src/check.rs @@ -244,8 +244,14 @@ impl Model { let needs_mass = self .steps .iter() - .find(|s| matches!(s.procedure, Procedure::Frequency { .. })) - .map(|s| format!("step `{}` extracts frequencies", s.name)) + .find_map(|s| { + let why = match s.procedure { + Procedure::Frequency { .. } => "extracts frequencies", + Procedure::ExplicitDynamic(_) => "is dynamic", + Procedure::Static {} => return None, + }; + Some(format!("step `{}` {why}", s.name)) + }) .or_else(|| self.loads.iter().find(|(_, l)| matches!(l, Load::Gravity { .. })).map(|(n, _)| format!("load `{n}` is gravity"))); if let Some(why) = needs_mass { let used: BTreeSet<&str> = self.sections.values().map(|s| s.material.as_str()).collect(); @@ -340,13 +346,18 @@ impl Model { } for s in &self.steps { let subject = format!("step `{}`", s.name); - match s.procedure { + let prescribes_motion = s + .boundary_conditions + .iter() + .filter_map(|b| self.boundary_conditions.get(b)) + .any(|b| (1..=6).any(|dof| b.value(dof).is_some_and(|v| v != 0.0))); + for set in s.procedure.sets() { + if set_of(set).nodes.is_empty() { + d.error(&subject, format!("set `{set}` has no nodes")); + } + } + match &s.procedure { Procedure::Static {} => { - let prescribes_motion = s - .boundary_conditions - .iter() - .filter_map(|b| self.boundary_conditions.get(b)) - .any(|b| (1..=6).any(|dof| b.value(dof).is_some_and(|v| v != 0.0))); if s.loads.is_empty() && !prescribes_motion { d.warning(&subject, "nothing loads the model in this step"); } @@ -356,8 +367,15 @@ impl Model { d.warning(&subject, "loads have no effect on natural frequencies"); } } + Procedure::ExplicitDynamic(e) => { + let moves = e.initial_velocities.iter().any(|v| v.velocity.iter().any(|c| *c != 0.0)); + if s.loads.is_empty() && !prescribes_motion && !moves { + d.warning(&subject, "nothing loads or moves the model in this step"); + } + } } - if s.boundary_conditions.is_empty() { + // A free body is a natural case of transient dynamics (flight, impact). + if s.boundary_conditions.is_empty() && !s.procedure.is_dynamic() { d.warning(&subject, "no boundary condition is active, so the model is free to move as a rigid body"); } } diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index 3420775..7260d59 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -56,7 +56,7 @@ pub use materials::{Elastic, Material}; pub use mesh::{Element, ElementFace, ElementKind, Node, check_coordinates}; pub use sections::{BeamProfile, Formulation, Integration, Section, SectionKind}; pub use sets::Set; -pub use steps::{OutputVariable, Procedure, Step}; +pub use steps::{BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, OutputVariable, Procedure, RigidWall, Step}; pub use units::UnitSystem; /// The largest node or element label. Nastran's GRID and element ids stop at 99 999 999 (eight diff --git a/crates/model/src/sets.rs b/crates/model/src/sets.rs index af62d71..fff5207 100644 --- a/crates/model/src/sets.rs +++ b/crates/model/src/sets.rs @@ -116,6 +116,7 @@ impl Model { users .extend(self.boundary_conditions.iter().filter(|(_, b)| b.target().set() == Some(set)).map(|(n, _)| format!("boundary condition `{n}`"))); users.extend(self.loads.iter().filter(|(_, l)| l.set() == Some(set)).map(|(n, _)| format!("load `{n}`"))); + users.extend(self.steps.iter().filter(|s| s.procedure.sets().any(|x| x == set)).map(|s| format!("step `{}`", s.name))); users } diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index 52fb47d..da8f77c 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -10,7 +10,7 @@ use std::collections::BTreeSet; use serde::{Deserialize, Serialize}; -use crate::{Error, Model, Result, check_name, context, finite, require}; +use crate::{Error, Model, Result, check_name, context, finite, positive, require}; /// Most steps in a model. pub const MAX_STEPS: usize = 1000; @@ -18,8 +18,20 @@ pub const MAX_STEPS: usize = 1000; /// Most modes a frequency step extracts. pub const MAX_MODES: u32 = 100_000; +/// Most time increments an explicit dynamic step may ask for. +pub const MAX_INCREMENTS: u64 = 1_000_000_000; + +/// Most field frames an explicit dynamic step writes (end time over output interval). +pub const MAX_FRAMES: f64 = 10_000.0; + +/// Most history samples an explicit dynamic step records (end time over history interval). +pub const MAX_HISTORY_SAMPLES: f64 = 1_000_000.0; + +/// Most initial velocities or rigid walls an explicit dynamic step lists. +pub const MAX_STEP_ITEMS: usize = 10_000; + /// What a step solves. -#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] #[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] pub enum Procedure { /// Linear static equilibrium under the active loads and boundary conditions. @@ -33,6 +45,233 @@ pub enum Procedure { #[serde(default, skip_serializing_if = "Option::is_none")] max_hz: Option, }, + /// Transient dynamics by explicit time integration up to `end_time`. + ExplicitDynamic(ExplicitDynamic), +} + +/// The settings of an explicit dynamic step. +/// +/// Loads act in full from the start of the step; a displacement boundary condition with a value +/// is reached at constant velocity over the step (value / `end_time`), a zero value holds the +/// node. The step starts from the undeformed model at rest, apart from `initial_velocities`. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct ExplicitDynamic { + /// Duration of the step (time units of the unit system). + pub end_time: f64, + /// Fraction of the stable time increment taken (0 < scale ≤ 1). + #[serde(default = "default_dt_scale")] + pub dt_scale: f64, + /// Cap on the time increments; the run stops early (and says so) when it is reached. + #[serde(default = "default_max_increments")] + pub max_increments: u64, + /// Time between field frames (displacement, velocity, stress); by default a twentieth of the + /// step. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub output_interval: Option, + /// Time between history samples (energies, momentum, constraint forces and the histories of + /// `history_nodes`); by default a thousandth of the step. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub history_interval: Option, + /// A set whose nodes' displacement and velocity are recorded as histories. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub history_nodes: Option, + /// Velocities of the nodes of sets at the start of the step (a node in several sets takes + /// the last one listed). + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub initial_velocities: Vec, + /// Fixed, frictionless rigid planes the nodes of sets cannot pass through. + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub rigid_walls: Vec, + /// Hourglass control of one-point hexahedra. + #[serde(default)] + pub hourglass: Hourglass, + /// Artificial bulk viscosity. + #[serde(default)] + pub bulk_viscosity: BulkViscosity, +} + +fn default_dt_scale() -> f64 { + 0.9 +} + +fn default_max_increments() -> u64 { + 10_000_000 +} + +impl ExplicitDynamic { + /// A step of `end_time` with every other setting at its default. + pub fn new(end_time: f64) -> ExplicitDynamic { + ExplicitDynamic { + end_time, + dt_scale: default_dt_scale(), + max_increments: default_max_increments(), + output_interval: None, + history_interval: None, + history_nodes: None, + initial_velocities: Vec::new(), + rigid_walls: Vec::new(), + hourglass: Hourglass::default(), + bulk_viscosity: BulkViscosity::default(), + } + } + + /// The time between field frames. + pub fn output_interval(&self) -> f64 { + self.output_interval.unwrap_or(self.end_time / 20.0) + } + + /// The time between history samples. + pub fn history_interval(&self) -> f64 { + self.history_interval.unwrap_or(self.end_time / 1000.0) + } + + /// The sets the step refers to. + pub fn sets(&self) -> impl Iterator { + self.history_nodes + .iter() + .map(String::as_str) + .chain(self.initial_velocities.iter().map(|v| v.set.as_str())) + .chain(self.rigid_walls.iter().map(|w| w.set.as_str())) + } + + pub fn validate(&self) -> Result<()> { + let end = positive("end time", self.end_time)?; + if !(self.dt_scale > 0.0 && self.dt_scale <= 1.0) { + return Err(Error::Invalid(format!("the time increment scale is in (0, 1], got {}", self.dt_scale))); + } + if !(1..=MAX_INCREMENTS).contains(&self.max_increments) { + return Err(Error::Invalid(format!("an explicit step takes 1 to {MAX_INCREMENTS} increments, got {}", self.max_increments))); + } + for (what, interval, most) in + [("output interval", self.output_interval, MAX_FRAMES), ("history interval", self.history_interval, MAX_HISTORY_SAMPLES)] + { + if let Some(i) = interval { + positive(what, i)?; + if end / i > most { + return Err(Error::Invalid(format!("the {what} {i} gives more than {most} samples over the step"))); + } + } + } + if self.initial_velocities.len() > MAX_STEP_ITEMS || self.rigid_walls.len() > MAX_STEP_ITEMS { + return Err(Error::Limit(format!("a step lists at most {MAX_STEP_ITEMS} initial velocities and rigid walls"))); + } + for v in &self.initial_velocities { + for c in v.velocity { + finite("initial velocity", c)?; + } + } + for w in &self.rigid_walls { + for c in w.point.iter().chain(&w.normal) { + finite("rigid wall", *c)?; + } + if !(w.normal.iter().map(|c| c * c).sum::() > 0.0) { + return Err(Error::Invalid("a rigid wall's normal must not be zero".into())); + } + } + self.hourglass.validate()?; + self.bulk_viscosity.validate() + } +} + +/// A velocity given to the nodes of a set at the start of an explicit step. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct InitialVelocity { + pub set: String, + pub velocity: [f64; 3], +} + +/// A fixed, frictionless rigid plane through `point`. Its `normal` points to the side the nodes +/// of `set` stay on. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct RigidWall { + pub set: String, + pub point: [f64; 3], + pub normal: [f64; 3], +} + +/// How one-point hexahedra resist their hourglass (zero-energy) modes. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Serialize, Deserialize)] +#[serde(rename_all = "snake_case")] +pub enum HourglassKind { + /// Forces proportional to the hourglass velocities. + #[default] + Viscous, + /// Forces proportional to the hourglass displacements. + Stiffness, +} + +/// Hourglass control: the kind and a dimensionless coefficient. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Hourglass { + #[serde(rename = "type", default)] + pub kind: HourglassKind, + #[serde(default = "default_hourglass_coefficient")] + pub coefficient: f64, +} + +/// Largest hourglass coefficient. +pub const MAX_HOURGLASS_COEFFICIENT: f64 = 0.25; + +fn default_hourglass_coefficient() -> f64 { + 0.1 +} + +impl Default for Hourglass { + fn default() -> Hourglass { + Hourglass { kind: HourglassKind::Viscous, coefficient: default_hourglass_coefficient() } + } +} + +impl Hourglass { + pub fn validate(&self) -> Result<()> { + if !(0.0..=MAX_HOURGLASS_COEFFICIENT).contains(&self.coefficient) { + return Err(Error::Invalid(format!("the hourglass coefficient is in 0..={MAX_HOURGLASS_COEFFICIENT}, got {}", self.coefficient))); + } + Ok(()) + } +} + +/// Artificial bulk viscosity in compression: a pressure `ρ L (quadratic L ε̇ᵥ² − linear c ε̇ᵥ)` +/// for a volumetric strain rate ε̇ᵥ < 0, element size L and wave speed c. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct BulkViscosity { + #[serde(default = "default_linear_viscosity")] + pub linear: f64, + #[serde(default = "default_quadratic_viscosity")] + pub quadratic: f64, +} + +/// Largest bulk viscosity coefficient. +pub const MAX_BULK_VISCOSITY: f64 = 10.0; + +fn default_linear_viscosity() -> f64 { + 0.06 +} + +fn default_quadratic_viscosity() -> f64 { + 1.5 +} + +impl Default for BulkViscosity { + fn default() -> BulkViscosity { + BulkViscosity { linear: default_linear_viscosity(), quadratic: default_quadratic_viscosity() } + } +} + +impl BulkViscosity { + pub fn validate(&self) -> Result<()> { + for (what, v) in [("linear", self.linear), ("quadratic", self.quadratic)] { + if !(0.0..=MAX_BULK_VISCOSITY).contains(&v) { + return Err(Error::Invalid(format!("the {what} bulk viscosity coefficient is in 0..={MAX_BULK_VISCOSITY}, got {v}"))); + } + } + Ok(()) + } } impl Procedure { @@ -41,10 +280,28 @@ impl Procedure { match self { Procedure::Static {} => "static", Procedure::Frequency { .. } => "frequency", + Procedure::ExplicitDynamic(_) => "explicit_dynamic", } } + /// Whether the procedure is transient dynamics (it needs mass; a free body is legitimate). + pub fn is_dynamic(&self) -> bool { + matches!(self, Procedure::ExplicitDynamic(_)) + } + + /// The sets the procedure's own settings refer to. + pub fn sets(&self) -> impl Iterator { + let explicit = match self { + Procedure::ExplicitDynamic(e) => Some(e), + Procedure::Static {} | Procedure::Frequency { .. } => None, + }; + explicit.into_iter().flat_map(ExplicitDynamic::sets) + } + pub fn validate(&self) -> Result<()> { + if let Procedure::ExplicitDynamic(e) = self { + return e.validate(); + } if let Procedure::Frequency { modes, min_hz, max_hz } = *self { if !(1..=MAX_MODES).contains(&modes) { return Err(Error::Invalid(format!("a frequency step extracts 1 to {MAX_MODES} modes, got {modes}"))); @@ -68,6 +325,7 @@ impl Procedure { match self { Procedure::Static {} => OutputVariable::ALL.into_iter().collect(), Procedure::Frequency { .. } => [OutputVariable::Displacement].into_iter().collect(), + Procedure::ExplicitDynamic(_) => [OutputVariable::Displacement, OutputVariable::Stress].into_iter().collect(), } } } @@ -121,6 +379,29 @@ impl Model { self.steps.iter().find(|s| s.name == name) } + /// The step named `name`, or the only step whose procedure is `wanted` when no name is + /// given; it must be such a step (`what` names the kind in messages, e.g. "static"). + pub fn select_step(&self, name: Option<&str>, wanted: fn(&Procedure) -> bool, what: &str) -> Result<&Step> { + let step = match name { + Some(n) => self.step(n).ok_or_else(|| Error::Invalid(format!("step `{n}` does not exist")))?, + None => { + let matching: Vec<&Step> = self.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) + } + fn step_index(&self, name: &str) -> Result { self.steps.iter().position(|s| s.name == name).ok_or_else(|| Error::NotFound(format!("step `{name}`"))) } @@ -131,6 +412,9 @@ impl Model { fn check_step(&self, step: &Step) -> Result<()> { step.procedure.validate().map_err(|e| context("step", &step.name, e))?; + for set in step.procedure.sets() { + require(&self.sets, "set", set).map_err(|e| context("step", &step.name, e))?; + } for b in &step.boundary_conditions { require(&self.boundary_conditions, "boundary condition", b).map_err(|e| context("step", &step.name, e))?; } @@ -230,4 +514,50 @@ mod tests { m.remove_step("Modes").unwrap(); assert_eq!(m.steps().len(), 1); } + + #[test] + fn explicit_steps_validate_and_hold_their_sets() { + let mut m = Model::default(); + m.add_set("front", Set { nodes: vec![1], ..Set::default() }).unwrap(); + let json = serde_json::json!({ + "type": "explicit_dynamic", "end_time": 1e-3, "history_nodes": "front", + "initial_velocities": [{"set": "front", "velocity": [-10, 0, 0]}], + "rigid_walls": [{"set": "front", "point": [0, 0, 0], "normal": [1, 0, 0]}], + "hourglass": {"type": "stiffness"} + }); + let p: Procedure = serde_json::from_value(json).unwrap(); + let Procedure::ExplicitDynamic(e) = &p else { panic!("{p:?}") }; + assert_eq!((e.dt_scale, e.hourglass.kind, e.hourglass.coefficient), (0.9, HourglassKind::Stiffness, 0.1)); + assert_eq!(e.bulk_viscosity, BulkViscosity { linear: 0.06, quadratic: 1.5 }); + assert!((e.output_interval() - 5e-5).abs() < 1e-18 && (e.history_interval() - 1e-6).abs() < 1e-18); + let back: Procedure = serde_json::from_value(serde_json::to_value(&p).unwrap()).unwrap(); + assert_eq!(back, p); + m.add_step(Step::new("Impact", p.clone()), None).unwrap(); + assert!(matches!(m.remove_set("front"), Err(Error::InUse(_))), "the step refers to the set"); + assert!(m.add_step(Step::new("Impact", p), None).is_err(), "names stay unique"); + let missing = serde_json::json!({"type": "explicit_dynamic", "end_time": 1.0, "history_nodes": "nope"}); + assert!(m.add_step(Step::new("Bad", serde_json::from_value(missing).unwrap()), None).is_err()); + + let base = ExplicitDynamic::new(1.0); + let bad = [ + ExplicitDynamic { end_time: 0.0, ..base.clone() }, + ExplicitDynamic { end_time: f64::NAN, ..base.clone() }, + ExplicitDynamic { dt_scale: 1.5, ..base.clone() }, + ExplicitDynamic { dt_scale: 0.0, ..base.clone() }, + ExplicitDynamic { max_increments: 0, ..base.clone() }, + ExplicitDynamic { output_interval: Some(1e-9), ..base.clone() }, + ExplicitDynamic { history_interval: Some(-1.0), ..base.clone() }, + ExplicitDynamic { rigid_walls: vec![RigidWall { set: "front".into(), point: [0.0; 3], normal: [0.0; 3] }], ..base.clone() }, + ExplicitDynamic { + initial_velocities: vec![InitialVelocity { set: "front".into(), velocity: [f64::INFINITY, 0.0, 0.0] }], + ..base.clone() + }, + ExplicitDynamic { hourglass: Hourglass { kind: HourglassKind::Viscous, coefficient: 2.0 }, ..base.clone() }, + ExplicitDynamic { bulk_viscosity: BulkViscosity { linear: -1.0, quadratic: 1.0 }, ..base.clone() }, + ]; + for b in bad { + assert!(Procedure::ExplicitDynamic(b.clone()).validate().is_err(), "{b:?}"); + } + assert!(serde_json::from_value::(serde_json::json!({"type": "explicit_dynamic", "end_time": 1, "dt": 2})).is_err()); + } } diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index edf7951..01a619b 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -121,29 +121,6 @@ pub struct FrequencyOptions { 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, @@ -218,7 +195,7 @@ fn budget(limits: &Limits, graph: &NodeGraph, matrices: usize) -> Budget { /// 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 step = model.select_step(step, |p| matches!(p, Procedure::Static {}), "static")?; let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); @@ -337,7 +314,7 @@ const TOLERANCE: f64 = 1e-9; /// 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 step = model.select_step(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())); }; diff --git a/crates/solver/tests/couplings.rs b/crates/solver/tests/couplings.rs index eacdb59..b3d5368 100644 --- a/crates/solver/tests/couplings.rs +++ b/crates/solver/tests/couplings.rs @@ -108,7 +108,7 @@ fn clamping_through_a_coupled_reference_point_is_clamping() { assert!((b.reaction_total[2] - 1.0).abs() < 1e-9, "{:?}", b.reaction_total); let modes = Procedure::Frequency { modes: 3, min_hz: None, max_hz: None }; - let (direct, coupled) = (build(false, modes), build(true, modes)); + let (direct, coupled) = (build(false, modes.clone()), build(true, modes)); let f = |g: &Grid| solve_frequency(&g.model, Some("Step"), &FrequencyOptions::default(), &Limits::default()).unwrap(); let (fa, fb) = (f(&direct), f(&coupled)); for (x, y) in fa.modes.iter().zip(&fb.modes) { diff --git a/docs/architecture.md b/docs/architecture.md index e7a8ead..6512e09 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -21,10 +21,11 @@ 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-fem` | Element kernels shared by every structural solver: the element table (tet4, tet10, wedge6, 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-fem` | Element kernels shared by every structural solver: the element table (tet4, tet10, wedge6, 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, and the uniform-strain (one-point) hex8 and tet4 kinematics, hourglass shape vectors and element frequency bounds of the explicit solver. 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-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 | +| L2 | `modelcraft-solver` | Implicit solver: assembly from the resolved model, boundary conditions, sparse direct solution, linear static and modal analysis, stresses and reactions. The lattice Boltzmann solver (`modelcraft-lbm`) joins this layer. | started | +| L2 | `modelcraft-explicit` | Explicit dynamics: central-difference time integration with lumped mass, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, rigid walls, energy account, histories and field frames. | 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 | @@ -170,6 +171,39 @@ modes), tension through a distributing coupling against the equivalent pressure end, pure bending by a moment through a distributing coupling, torsion of a tube by a moment on a kinematically coupled point (twist and shear stress). +## Explicit solver + +`modelcraft-explicit` integrates a model's `explicit_dynamic` steps in time (`solve.explicit`; +`solve.history` and `solve.probe` read the stored histories and frames). Methods and papers are +cited in the module docs and in docs/solvers/literature.md (Part H). + +1. **System.** `Model::resolve()`; hex8 sections with reduced integration become one-point + uniform-strain hexahedra with Flanagan–Belytschko hourglass control, tet4 sections Bonet–Burton + average nodal pressure tetrahedra (`modelcraft_fem::uniform` holds their kinematics). Other + element kinds and formulations, couplings, reference points and moments return "not supported + yet". Masses are row-sum lumped. A zero-valued boundary condition holds a degree of freedom; a + displacement value is reached at constant velocity over the step. Loads act in full from the + start. The step lists initial velocities and fixed frictionless rigid planes by set. +2. **Time step.** Each element's highest frequency is bounded in closed form (exact for a cube); + the increment is the step's scale (0.9 by default) times the smallest element increment with + Wilkins' bulk-viscosity correction, recomputed every increment. The ten controlling elements + are reported; a step that would need more increments than its cap is refused before it runs. +3. **Integration.** Central differences in the half-step form, elements in label order, no + parallel reductions (bit-for-bit repeatable). Rigid walls and prescribed degrees of freedom are + kinematic velocity constraints whose impulses are the reactions. +4. **Energy account.** Kinetic, internal (stress work), hourglass, bulk viscosity, external and + constraint work, and their balance; a run warns when the balance is off by more than 1% of the + largest energy or the hourglass energy exceeds 10% of the internal energy. +5. **Results.** Histories (energies, momentum, constraint force, chosen nodes) at the history + interval and field frames (displacement, velocity, element stress and strain) at the output + interval, capped in size; kept in the session like the implicit results. + +The first version is geometrically linear with linear elastic materials. Verification +(`crates/explicit/tests/verification.rs`, analytic only): stress waves in hex8 and tet4 bars +(front at c·t, σ = ρ c v, work input), a bar striking a rigid wall (contact time 2L/c, rebound +speed, plateau force ρ c v₀ A), the energy balance, bit-for-bit repeatability; in the engine +(`solve.rs`), a cantilever's free-vibration period against `solve.frequency`. + ## Commands Every user-visible action is a `CommandSpec` in `crates/engine/src/cmd/.rs` with a JSON diff --git a/docs/features.md b/docs/features.md index 3468cc0..7b76de9 100644 --- a/docs/features.md +++ b/docs/features.md @@ -5,7 +5,7 @@ Generated from the capability catalogue (`crates/engine/src/catalogue.rs`) by Every entry is a toolbar button. **Ready** entries run their command; **coming soon** entries are shown greyed with the roadmap milestone that brings them (ROADMAP.md). -35 ready, 65 coming soon. +36 ready, 64 coming soon. ## Geometry @@ -121,7 +121,7 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md). | Decks | Read Deck — Read a deck, keeping keywords ModelCraft does not know | coming soon (M9) | | | Implicit | Static — Linear static solve of a step with ModelCraft's own solver | ready | `solve.static` | | Implicit | Frequency — Natural frequencies and mode shapes of a frequency step | ready | `solve.frequency` | -| Explicit | Explicit — Explicit dynamics: central differences, contact, erosion | coming soon (M10) | | +| Explicit | Explicit — Transient dynamics of an explicit step by central differences (contact and erosion coming in M10) | ready | `solve.explicit` | | Explicit | Blast & Impact — EOS-driven fluids with fluid–structure contact (ALE or SPH) | coming soon (M14) | | | Fluids | Lattice Boltzmann — Lattice Boltzmann flow on the CPU or the GPU | coming soon (M12) | | | Fluids | Fluid–Structure — Coupled fluid and structure solve | coming soon (M13) | | @@ -134,7 +134,7 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md). | Display | Contours — Contour plots of stress, strain and displacement | coming soon (M4) | | | Display | Mode Shapes — Mode shapes of a frequency step | coming soon (M4) | | | Display | Animate — Animate modes and time histories | coming soon (M9) | | -| Query | Probe — Displacement, reaction and stress at nodes, or a mode shape | ready | `solve.probe` | +| Query | Probe — Displacement, reaction and stress at nodes, a mode shape, or an explicit field frame | ready | `solve.probe`, `solve.history` | | Query | Reactions — Reaction force and moment sums | coming soon (M4) | | | Query | Free Body — Free-body forces across a cut | coming soon (M9) | | | Plots | XY Plot — Plot a quantity against time or another | coming soon (M9) | | diff --git a/docs/parity.md b/docs/parity.md index 00e5332..979fc82 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -61,6 +61,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | 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 | started (solid elements: tet4, tet10, wedge6, hex8, hex20; reference points, kinematic and distributing couplings, moments) | +| Solve | Own explicit dynamics solver (central differences, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, energy balance), verified against analytic wave and impact solutions | M10 | yes | started (X1 core: small-strain elastic, rigid walls; large deformation, plasticity and contact planned) | | 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 | diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index 01f3731..bc964e2 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -991,3 +991,34 @@ the module docs): element accuracy*, Finite Elem. Anal. Des. 1(1):3–20, 1985 (patch tests); Lamé, Kirsch, Timoshenko and Euler–Bernoulli closed forms; Saint-Venant's principle (tension through a coupling), pure bending σ = M z / I, torsion of a circular tube φ = T L / (G J), τ = T r / J. + +## Part H: Explicit solver core (X1), as implemented + +The papers and books `modelcraft-explicit` and `modelcraft_fem::uniform` are written from (cited +again in the module docs): + +- Time integration, energy balance, element eigenvalue inequality: T. Belytschko, W. K. Liu, + B. Moran, K. Elkhodary, *Nonlinear Finite Elements for Continua and Structures*, 2nd ed., Wiley + 2014 (ch. 6: central differences in the half-step form with variable increments, the energy + balance check, the critical time step). +- One-point hexahedron and hourglass control: D. P. Flanagan, T. Belytschko, IJNME + 17(5):679–706, 1981 (integrated gradients b, mean strain, the γ hourglass shape vectors + orthogonal to linear fields, viscous and stiffness resisting forces). Our scaling of the + coefficient with ρ, c and the element size is our own: κ is the fraction of critical damping + (viscous) or the squared frequency ratio (stiffness) of a cube's hourglass modes at 2c/h. +- Element frequency bound: in the spirit of D. P. Flanagan, T. Belytschko, J. Appl. Mech. + 51(1):35–40, 1984; we bound uᵀKu ≤ (k/V) λ_max(Σ b bᵀ) ‖u‖² with k = max(3λ + 2μ, 2μ) + (isotropic) or 2 λ_max(D), which is exact for a cube; the 3 × 3 eigenvalue in closed form after + O. K. Smith, Commun. ACM 4(4):168, 1961. +- Average nodal pressure tetrahedron: J. Bonet, A. J. Burton, Commun. Numer. Methods Eng. + 14(5):437–449, 1998. +- Bulk viscosity and its time-step correction: J. von Neumann, R. D. Richtmyer, J. Appl. Phys. + 21(3):232–237, 1950; M. L. Wilkins, J. Comput. Phys. 36(3):281–303, 1980. +- Rigid walls as kinematic velocity constraints: T. J. R. Hughes, R. L. Taylor, J. L. Sackman, + A. Curnier, W. Kanoknukulchai, CMAME 8(3):249–276, 1976. +- Verification: d'Alembert's solution of the 1D wave equation (front at c t, σ = ρ c v), the + elastic bar impact on a rigid wall (contact 2L/c, rebound at the impact speed; Hughes et al. + 1976), Euler–Bernoulli and the implicit modal solver for the cantilever period. The one-point + element keeps 1 − 1/n² of the bending energy with n elements through the depth, which predicts + its period with viscous hourglass control to 0.03%. + diff --git a/xtask/src/layers.rs b/xtask/src/layers.rs index 94927e3..0b0d910 100644 --- a/xtask/src/layers.rs +++ b/xtask/src/layers.rs @@ -43,8 +43,9 @@ pub const TABLE: &[(&str, Class)] = &[ ("fem", Class::Layer(1)), ("decks", Class::Layer(1)), ("render", Class::Layer(1)), - // L2: solvers (the explicit and lattice Boltzmann solvers join here). + // L2: solvers (the lattice Boltzmann solver joins here). ("solver", Class::Layer(2)), + ("explicit", Class::Layer(2)), // L3: coupling between solvers (`fsi`, later). // L4 and up: the document, files, commands, UI. ("doc", Class::Layer(4)),