From 29ae46ff84a2238e5c75a1225f10a5fb630944cf Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Tue, 13 Oct 2026 22:00:25 -0400 Subject: [PATCH] Shells and beams: MITC4/MITC3 shells, Timoshenko beams, six-DOF nodes (M4) (#47) * fem: MITC4/MITC3 flat shells and Timoshenko beams with profile properties (ElementDofs STRUCTURAL) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * solver: shells and beams on the DofLayout beside solids, 2D elements and trusses; analytic verification and tagged benchmarks Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm * engine: probe rotations, reaction moments and shell and beam element results; docs for shells and beams Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --------- Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 --- ROADMAP.md | 2 +- crates/engine/src/catalogue.rs | 2 +- crates/engine/src/cmd/solve.rs | 95 ++++- crates/fem/src/beam.rs | 429 +++++++++++++++++++++ crates/fem/src/frame.rs | 142 +++++++ crates/fem/src/lib.rs | 9 +- crates/fem/src/material.rs | 20 +- crates/fem/src/shell.rs | 593 ++++++++++++++++++++++++++++++ crates/solver/src/analysis.rs | 115 +++++- crates/solver/src/constraints.rs | 16 +- crates/solver/src/system.rs | 204 +++++++++- crates/solver/tests/benchmarks.rs | 182 +++++++++ crates/solver/tests/structural.rs | 591 +++++++++++++++++++++++++++++ docs/architecture.md | 33 +- docs/elements.md | 4 +- docs/features.md | 2 +- docs/parity.md | 2 +- docs/solvers/literature.md | 18 + 18 files changed, 2411 insertions(+), 48 deletions(-) create mode 100644 crates/fem/src/beam.rs create mode 100644 crates/fem/src/frame.rs create mode 100644 crates/fem/src/shell.rs create mode 100644 crates/solver/tests/benchmarks.rs create mode 100644 crates/solver/tests/structural.rs diff --git a/ROADMAP.md b/ROADMAP.md index 9404ce6..05a322d 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -22,7 +22,7 @@ hex meshes today and it gives the best elements where stresses are read. | M1 Geometry in and surface preparation | Open `.solidcraft` designs with persistent face names, assemblies placed (done); face-tagged 3MF (done); test corpus, `cargo xtask corpus` (done); exact predicates (done); quality metrics and certified validity gate (done); weld, closed/manifold/orientation and self-intersection checks, feature edges (`modelcraft-mesh::surface`, `mesh.surface.check`, done) | in progress | | M2 Tet meshing and size field | Size field (global, per face and edge seeds, curvature and proximity; the `modelcraft_mesh::size::SizeField` trait and edge-division helpers are in); Bowyer–Watson Delaunay (exact, symbolic perturbation, verified; `modelcraft-mesh::delaunay`, done); robust CDT boundary recovery (Diazzi et al. 2023: implicit LNC Steiner points on segments, gift-wrapping facet recovery, inside/outside labelling, B-rep face tags on boundary tet faces, few-ulp rounding repair checked by the validity gate; `modelcraft-mesh::cdt` and `mesh.tet.generate`, done: all 20 valid corpus parts); Delaunay refinement, ODT smoothing, sliver removal; tet4 and tet10 with mid-side nodes on the B-rep; mesh verify per set | in progress | | M3 First deck | Sets and surfaces bound to SolidCraft face, edge, vertex and body names, resolved after meshing; materials (isotropic and orthotropic elastic, density) with a unit-system label and consistency check; sections (solid, shell, beam profiles) assigned by set; element formulation per region mapped to each solver; BCs (fixed, displacement, symmetry); loads (force and moment on reference points, pressure, gravity); reference points with kinematic and distributing coupling; steps (linear static, frequency) with explicit per-step load and BC tables and output presets; `model.check` naming every missing or bad item; `model.*` and `mesh.tag` commands (model definition: done); writers for Abaqus input (done: `deck.write`), Nastran SOL 101/103 (done: `deck.write`) and LS-DYNA implicit (done); the `.modelcraft` file (text definition, binary mesh, geometry references reopened on load, journal; atomic saves, autosave and recovery; undo/redo on the journal: done) | in progress | -| M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 1D elements I: rod2/rod3 trusses, point masses, rotary inertia, springs and dashpots in both solvers and all three decks, verified by the method of joints, the three-bar truss, spring–mass, two-mass chain and torsional pendulum frequencies, and damped decay; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started | +| M4 Own solver and first results | Element library per docs/elements.md (1D and 2D first); Linear static and modal solver (tet4/10, wedge6, hex8, hex20; sparse direct, Lanczos), verified against analytic solutions (solver: done — `modelcraft-fem` kernels, `modelcraft-solver` with nested dissection, sparse Cholesky and block Lanczos, `solve.static`/`solve.frequency`/`solve.probe`, patch, cantilever, Lamé, Kirsch and beam-frequency tests; reference points with kinematic and distributing couplings and moments by master–slave constraint elimination, verified by rigid motion, clamping, Saint-Venant tension, bending and tube torsion; 1D elements I: rod2/rod3 trusses, point masses, rotary inertia, springs and dashpots in both solvers and all three decks, verified by the method of joints, the three-bar truss, spring–mass, two-mass chain and torsional pendulum frequencies, and damped decay; 2D elements: plane stress, plane strain and axisymmetric tri3, tri6, quad4 (incompatible modes QM6), quad8 with two degrees of freedom per node, centrifugal loads, `mesh.planar` 2D meshes of polygons and planar faces, verified by patch tests, Kirsch, Lamé (plane strain and axisymmetric), rotating discs and cylinders, and a plane stress cantilever; deck mapping to Abaqus CPS/CPE/CAX, Nastran membranes, LS-DYNA ELFORM 12/13); shells and beams (done: quad4 MITC4 and tri3 MITC3 flat shells with a Hughes–Brezzi drilling rotation, two-node Timoshenko beams with profile properties and orientations, six degrees of freedom at shell and beam nodes shared with solids and trusses, section forces and surface stresses, verified by shell patch tests, Navier's plate, cantilever and portal-frame beams, and the Scordelis–Lo and pinched-cylinder benchmarks); results display (done: `results.*` commands, deformed shape to scale, banded contours of displacement, von Mises, principal and component stresses with a legend, mode shapes with their frequencies and animation, probing a node by clicking it, reaction force and moment sums per set, stale results marked after an edit; explicit steps: field frames by slider or animation with contours per frame, energy histories plotted) | started | | M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | | M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, connectors; amplitudes, multiple steps; nonlinear static steps in decks | planned | diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index 6af762d..4a8fb7a 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -366,7 +366,7 @@ pub static CATALOGUE: &[Capability] = &[ "Probe", Results, "Query", - "Displacement, reaction and stress at nodes, a mode shape, or an explicit field frame (or click the model)", + "Displacement, rotation, reaction and stress at nodes, element section forces, a mode shape, or an explicit field frame (or click the model)", &["solve.history"], ), ready("results.reactions", "Reactions", Results, "Query", "Reaction force and moment sums per set"), diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index bc542e5..d2afdd5 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -5,7 +5,7 @@ use modelcraft_explicit::{ExplicitResult, solve_explicit}; use modelcraft_lbm::{FluidResult, Obstacle}; use modelcraft_mesh::surface::{BoundarySurface, SurfaceOptions}; use modelcraft_model::FluidDomain; -use modelcraft_solver::{FrequencyOptions, Limits, Location, MassKind, StepResult as ImplicitResult, solve_frequency, solve_static}; +use modelcraft_solver::{FrequencyOptions, Limits, Location, MassKind, StaticResult, StepResult as ImplicitResult, solve_frequency, solve_static}; use serde_json::{Value, json}; use super::{labels, opt, points, req}; @@ -54,7 +54,7 @@ pub const FLUID: CommandSpec = CommandSpec { pub const PROBE: CommandSpec = CommandSpec { id: "solve.probe", label: "Probe Results", - 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) at `nodes` (or at the nodes nearest to points: "at": [[x, y, z], ...], lowest label on a tie, each entry then giving its `xyz` and `distance`) and, with `elements` (explicit), each element's stress, von Mises stress, strain, equivalent plastic strain, damage, temperature and whether it has eroded; a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity and pressure in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories"#, + params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1, "frame": 0} — the stored result at nodes: displacement, rotation, reaction, reaction moment, stress and von Mises (static, with `elements` each element's results: solid, 2D and truss nodal stresses, shell section forces (N, M, Q per length in the element's axes) and bottom, middle and top stresses, beam end forces (N, V_y, V_z, T, M_y, M_z in its axes)), the shape of mode `mode` with its rotations (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default) at `nodes` (or at the nodes nearest to points: "at": [[x, y, z], ...], lowest label on a tie, each entry then giving its `xyz` and `distance`) and, with `elements` (explicit), each element's stress, von Mises stress, strain, equivalent plastic strain, damage, temperature and whether it has eroded; a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity and pressure in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories"#, kind: CommandKind::Model, run: probe, }; @@ -355,7 +355,7 @@ fn probe(s: &mut Session, p: &Value) -> Result { let found = nearest_nodes(&s.model, result.nodes(), &at)?; (found.iter().map(|(l, _, _)| *l).collect(), Some(found)) } - (None, None) if matches!(result.as_ref(), StepResult::Explicit(_)) => (Vec::new(), None), + (None, None) if matches!(result.as_ref(), StepResult::Explicit(_)) || p.get("elements").is_some() => (Vec::new(), None), _ => return Err(EngineError::Params("give either `nodes` or `at`".into())), }; let mut out = match result.as_ref() { @@ -393,6 +393,7 @@ fn nearest_nodes(model: &modelcraft_model::Model, with_result: &[u32], at: &[[f6 /// The result of an implicit step at nodes. fn probe_implicit(result: &ImplicitResult, step: &str, nodes: &[u32], p: &Value) -> Result { + let elements = if p.get("elements").is_some() { labels(p, "elements", MAX_PROBE)? } else { Vec::new() }; let locate = |label: u32| { result .locate(label) @@ -410,7 +411,9 @@ fn probe_implicit(result: &ImplicitResult, step: &str, nodes: &[u32], p: &Value) Location::Node(i) => json!({ "node": l, "displacement": r.displacement.get(i), + "rotation": r.rotation.get(i), "reaction": r.reaction.get(i), + "reaction_moment": r.reaction_moment.get(i), "stress": r.stress.get(i), "von_mises": r.von_mises.get(i), }), @@ -418,6 +421,9 @@ fn probe_implicit(result: &ImplicitResult, step: &str, nodes: &[u32], p: &Value) }) .collect::>()?, ImplicitResult::Frequency(r) => { + if !elements.is_empty() { + return Err(EngineError::Params("element results are kept for static and explicit steps".into())); + } let mode = opt::(p, "mode")?.unwrap_or(1); let m = mode .checked_sub(1) @@ -430,13 +436,43 @@ fn probe_implicit(result: &ImplicitResult, step: &str, nodes: &[u32], p: &Value) Location::Point(k) => { json!({"node": l, "reference_point": r.reference_points.get(k).map(|p| &p.0), "shape": m.point_shape.get(k)}) } - Location::Node(i) => json!({"node": l, "shape": m.shape.get(i)}), + Location::Node(i) => json!({"node": l, "shape": m.shape.get(i), "rotation": m.rotation.get(i)}), }) }) .collect::>()? } }; - Ok(json!({"step": step, "nodes": out})) + let mut response = json!({"step": step, "nodes": out}); + if let ImplicitResult::Static(r) = result + && !elements.is_empty() + { + let found: Vec = elements + .iter() + .map(|&l| element_result(r, l).ok_or_else(|| EngineError::Params(format!("element {l} has no result in step `{step}`")))) + .collect::>()?; + response["elements"] = Value::Array(found); + } + Ok(response) +} + +/// The entry of a label-sorted result list. +fn by_label(list: &[T], label: u32, key: impl Fn(&T) -> u32) -> Option<&T> { + list.binary_search_by_key(&label, key).ok().and_then(|i| list.get(i)) +} + +/// The stored results of one element of a static step. +fn element_result(r: &StaticResult, label: u32) -> Option { + if let Some(e) = by_label(&r.elements, label, |e| e.label) { + return Some(json!({"element": label, "kind": e.kind.name(), "stress": e.nodal})); + } + if let Some(e) = by_label(&r.shells, label, |e| e.label) { + let s = &e.stress; + return Some( + json!({"element": label, "kind": e.kind.name(), "axes": s.axes, "forces": s.forces, "bottom": s.bottom, "middle": s.middle, "top": s.top}), + ); + } + let e = by_label(&r.beams, label, |e| e.label)?; + Some(json!({"element": label, "kind": e.kind.name(), "axes": e.forces.axes, "ends": e.forces.ends})) } /// Displacement and velocity at nodes, and element results, in one field frame of an explicit @@ -1026,4 +1062,53 @@ pub(crate) mod tests { s.execute("model.spring.delete", &json!({"name": "c"})).unwrap(); assert!(s.model.springs().len() == 1 && s.model.mass("m").is_some()); } + #[test] + fn solve_and_probe_shells_and_beams() { + // A cantilever strip of two quad4 shells clamped at x = 0, stiffened along its tip edge + // by a beam; a downward force at a tip corner. + let mut s = Session::default(); + let mut script = Vec::new(); + for (label, xyz) in [(1, [0, 0, 0]), (2, [10, 0, 0]), (3, [20, 0, 0]), (4, [0, 5, 0]), (5, [10, 5, 0]), (6, [20, 5, 0])] { + script.push(json!({"command": "mesh.node.add", "params": {"label": label, "xyz": xyz}})); + } + script.extend([ + json!({"command": "mesh.element.add", "params": {"label": 1, "kind": "quad4", "nodes": [1, 2, 5, 4]}}), + json!({"command": "mesh.element.add", "params": {"label": 2, "kind": "quad4", "nodes": [2, 3, 6, 5]}}), + json!({"command": "mesh.element.add", "params": {"label": 3, "kind": "beam2", "nodes": [3, 6]}}), + json!({"command": "model.set.create", "params": {"name": "plate", "elements": [1, 2]}}), + json!({"command": "model.set.create", "params": {"name": "stiffener", "elements": [3]}}), + json!({"command": "model.set.create", "params": {"name": "root", "nodes": [1, 4]}}), + json!({"command": "model.set.create", "params": {"name": "corner", "nodes": [3]}}), + json!({"command": "model.material.create", "params": {"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3}, "density": 7.85e-9}}), + json!({"command": "model.section.create", "params": {"name": "skin", "set": "plate", "material": "Steel", "kind": {"type": "shell", "thickness": 1}}}), + json!({"command": "model.section.create", "params": {"name": "rib", "set": "stiffener", "material": "Steel", "kind": {"type": "beam", "profile": {"type": "rectangle", "width": 1, "height": 2}, "orientation": [1, 0, 0]}}}), + json!({"command": "model.boundary_condition.create", "params": {"name": "clamp", "type": "fixed", "target": {"set": "root"}}}), + json!({"command": "model.load.create", "params": {"name": "push", "type": "force", "target": {"set": "corner"}, "vector": [0, 0, -1]}}), + json!({"command": "model.step.create", "params": {"name": "Bend", "procedure": {"type": "static"}, "boundary_conditions": ["clamp"], "loads": ["push"]}}), + json!({"command": "model.step.create", "params": {"name": "Modes", "procedure": {"type": "frequency", "modes": 3}, "boundary_conditions": ["clamp"]}}), + ]); + s.run_script(&Value::Array(script)).unwrap(); + let r = s.execute("solve.static", &json!({"step": "Bend"})).unwrap(); + assert!((r["reaction_total"][2].as_f64().unwrap() - 1.0).abs() < 1e-9, "{r}"); + let probe = s.execute("solve.probe", &json!({"step": "Bend", "nodes": [3, 1], "elements": [1, 3]})).unwrap(); + assert!(probe["nodes"][0]["displacement"][2].as_f64().unwrap() < 0.0, "{probe}"); + assert!(probe["nodes"][0]["rotation"][0].as_f64().unwrap() != 0.0, "{probe}"); + assert!(probe["nodes"][1]["reaction_moment"].is_array(), "{probe}"); + assert_eq!(probe["elements"][0]["forces"].as_array().unwrap().len(), 4, "{probe}"); + assert_eq!(probe["elements"][1]["ends"].as_array().unwrap().len(), 2, "{probe}"); + let only = s.execute("solve.probe", &json!({"step": "Bend", "elements": [2]})).unwrap(); + assert_eq!(only["nodes"].as_array().unwrap().len(), 0, "{only}"); + s.execute("solve.frequency", &json!({"step": "Modes"})).unwrap(); + let shape = s.execute("solve.probe", &json!({"step": "Modes", "nodes": [3], "mode": 1})).unwrap(); + assert!(shape["nodes"][0]["rotation"].is_array(), "{shape}"); + for p in [ + json!({"step": "Bend", "elements": [99]}), + json!({"step": "Bend", "elements": [0]}), + json!({"step": "Bend", "elements": "1"}), + json!({"step": "Bend", "elements": vec![1; 20_000]}), + json!({"step": "Modes", "nodes": [3], "elements": [1]}), + ] { + assert!(s.execute("solve.probe", &p).is_err(), "{p}"); + } + } } diff --git a/crates/fem/src/beam.rs b/crates/fem/src/beam.rs new file mode 100644 index 0000000..a4f40b0 --- /dev/null +++ b/crates/fem/src/beam.rs @@ -0,0 +1,429 @@ +//! The two-node Timoshenko beam (six degrees of freedom per node), with the section properties of +//! the model's beam profiles. (Axial-only members are the trusses of [`crate::truss`].) +//! +//! **Local axes.** x runs from node 1 to node 2; y is the section's orientation vector (the +//! profile's local 1 direction, its `width`) made normal to x; z = x × y (the profile's local 2 +//! direction, its `height`). Degrees of freedom per node: translations along x, y, z, then +//! rotations about them; element matrices are built in local axes and rotated to the global ones +//! ([`crate::frame`]). +//! +//! **Beam.** Axial and torsional stiffness are the exact two-node ones (EA/L, GJ/L). Each bending +//! plane uses the interdependent interpolation of the Timoshenko beam: cubic deflection and +//! quadratic rotation, coupled through Φ = 12 EI / (G A_s L²) so the shear strain is constant and +//! the element is exact for end loads (any number of elements gives the nodal displacements of +//! Timoshenko beam theory; J. S. Przemieniecki, *Theory of Matrix Structural Analysis*, +//! McGraw-Hill 1968, §5.6; Z. Friedman, J. B. Kosmatka, *An improved two-node Timoshenko beam +//! finite element*, Computers & Structures 47(3):473–481, 1993, whose shape functions are used +//! here). With these functions w(ξ) and θ(ξ), ξ = x/L, +//! K = ∫ (EI θ'ᵀ θ' + G A_s (w' − θ)ᵀ (w' − θ)) dx and the consistent mass +//! M = ∫ (ρA N_wᵀ N_w + ρI N_θᵀ N_θ) dx (translational and rotary inertia) are integrated with +//! four Gauss points, exactly; Φ → 0 recovers the Euler–Bernoulli (Hermite) element. The +//! rotation about y is −dw/dx in the x–z plane, so that plane's rotations enter with a minus +//! sign. Lumped mass: half the mass on each node's translations, half the polar inertia on the +//! twist, and the consistent bending rotary diagonal scaled like the translations (Hinton, Rock & +//! Zienkiewicz 1976). End forces are K u minus the consistent loads of the element's body force, +//! reported as the section forces at each end (axial force N, shears V_y, V_z, torque T, moments +//! M_y, M_z on the cut whose outward normal is +x). +//! +//! **Section properties.** Area, second moments I_y = ∫ z² dA and I_z = ∫ y² dA, the polar +//! moment I_y + I_z (rotary inertia of torsion), the torsion constant J and the shear areas: +//! - rectangle: J from Saint-Venant's series, J = a b³/3 (1 − 192 b/(π⁵ a) Σ_{n odd} +//! tanh(nπa/2b)/n⁵) for sides a ≥ b (S. P. Timoshenko, J. N. Goodier, *Theory of Elasticity*, +//! 3rd ed., McGraw-Hill 1970, §109); shear coefficient 10(1 + ν)/(12 + 11ν); +//! - circle and pipe: J = I_y + I_z; shear coefficients 6(1 + ν)/(7 + 6ν) and, for a tube of +//! radius ratio m, 6(1 + ν)(1 + m²)² / ((7 + 6ν)(1 + m²)² + (20 + 12ν) m²) (G. R. Cowper, +//! *The shear coefficient in Timoshenko's beam theory*, J. Appl. Mech. 33(2):335–340, 1966); +//! - I and box (thin-walled): J = Σ b t³/3 over the midline walls (open I) and Bredt's +//! 4 A_m² t / p (closed box); shear areas the walls parallel to the shear, between the midlines +//! of the walls across it (the web for shear along z, 5/6 of the flanges along y; the two +//! webs or the two flanges of a box) (W. D. Pilkey, *Analysis and Design of Elastic Beams*, +//! Wiley 2002, ch. 6). + +use modelcraft_model::{BeamProfile, Elastic}; + +use crate::dense::Matrix; +use crate::dof::{DofSet, ElementDofs}; +use crate::frame::{Axes, Vec3, cross, dot, norm, rotate_matrix, rotate_to_global, rotate_to_local, sub, to_local, translation, unit}; +use crate::quadrature::gauss_legendre; +use crate::{Error, Result}; + +/// Cross-section properties in the beam's local axes (y along the profile's width, z along its +/// height). +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Properties { + pub area: f64, + /// ∫ z² dA: bending with deflection along z. + pub iy: f64, + /// ∫ y² dA: bending with deflection along y. + pub iz: f64, + /// Saint-Venant torsion constant. + pub torsion: f64, + /// Shear areas for shear along y and along z. + pub shear_y: f64, + pub shear_z: f64, +} + +/// Saint-Venant's torsion constant of a solid rectangle. +fn rectangle_torsion(width: f64, height: f64) -> f64 { + let (a, b) = (width.max(height), width.min(height)); + let pi = std::f64::consts::PI; + let series: f64 = (0..50).map(|k| f64::from(2 * k + 1)).map(|n| (n * pi * a / (2.0 * b)).tanh() / n.powi(5)).sum(); + a * b * b * b / 3.0 * (1.0 - 192.0 * b / (pi.powi(5) * a) * series) +} + +/// The section properties of a profile (`poisson` enters the shear coefficients). +pub fn properties(profile: &BeamProfile, poisson: f64) -> Properties { + let pi = std::f64::consts::PI; + let area = profile.area(); + match *profile { + BeamProfile::Rectangle { width, height } => { + let k = 10.0 * (1.0 + poisson) / (12.0 + 11.0 * poisson); + Properties { + area, + iy: width * height.powi(3) / 12.0, + iz: height * width.powi(3) / 12.0, + torsion: rectangle_torsion(width, height), + shear_y: k * area, + shear_z: k * area, + } + } + BeamProfile::Circle { radius } => { + let i = pi * radius.powi(4) / 4.0; + let k = 6.0 * (1.0 + poisson) / (7.0 + 6.0 * poisson); + Properties { area, iy: i, iz: i, torsion: 2.0 * i, shear_y: k * area, shear_z: k * area } + } + BeamProfile::Pipe { radius, thickness } => { + let inner = radius - thickness; + let i = pi * (radius.powi(4) - inner.powi(4)) / 4.0; + let m2 = (inner / radius).powi(2); + let k = 6.0 * (1.0 + poisson) * (1.0 + m2).powi(2) / ((7.0 + 6.0 * poisson) * (1.0 + m2).powi(2) + (20.0 + 12.0 * poisson) * m2); + Properties { area, iy: i, iz: i, torsion: 2.0 * i, shear_y: k * area, shear_z: k * area } + } + BeamProfile::I { height, width, web, flange } => { + let inside = height - 2.0 * flange; + Properties { + area, + iy: (width * height.powi(3) - (width - web) * inside.powi(3)) / 12.0, + iz: (2.0 * flange * width.powi(3) + inside * web.powi(3)) / 12.0, + torsion: (2.0 * width * flange.powi(3) + (height - flange) * web.powi(3)) / 3.0, + shear_y: 5.0 / 6.0 * 2.0 * width * flange, + shear_z: (height - flange) * web, + } + } + BeamProfile::Box { width, height, thickness: t } => { + let (wi, hi) = (width - 2.0 * t, height - 2.0 * t); + let (wm, hm) = (width - t, height - t); + Properties { + area, + iy: (width * height.powi(3) - wi * hi.powi(3)) / 12.0, + iz: (height * width.powi(3) - hi * wi.powi(3)) / 12.0, + torsion: 4.0 * (wm * hm).powi(2) * t / (2.0 * (wm + hm)), + shear_y: 2.0 * wm * t, + shear_z: 2.0 * hm * t, + } + } + } +} + +/// The Young's and shear moduli and Poisson's ratio of an isotropic material. +/// (`what` names the section kind in errors.) +fn moduli(elastic: &Elastic, what: &str) -> Result<(f64, f64, f64)> { + let (young, poisson) = crate::material::isotropic_constants(elastic, what)?; + Ok((young, young / (2.0 * (1.0 + poisson)), poisson)) +} + +/// A beam section: profile properties and an isotropic material. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct BeamSection { + pub properties: Properties, + pub young: f64, + pub shear_modulus: f64, +} + +impl BeamSection { + pub fn new(profile: &BeamProfile, elastic: &Elastic) -> Result { + profile.validate().map_err(|e| Error::Invalid(e.to_string()))?; + let (young, shear_modulus, poisson) = moduli(elastic, "beam")?; + Ok(BeamSection { properties: properties(profile, poisson), young, shear_modulus }) + } +} + +/// The section forces at the two ends of a beam: N, V_y, V_z, T, M_y, M_z in its local axes. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct BeamForces { + pub axes: Axes, + pub ends: [[f64; 6]; 2], +} + +/// One bending plane's stiffness and mass (deflection, rotation at node 1, then node 2), with +/// θ = dw/dx in the plane's own convention. +fn bending_plane(ei: f64, gas: f64, rho_a: f64, rho_i: f64, l: f64) -> ([[f64; 4]; 4], [[f64; 4]; 4]) { + let phi = 12.0 * ei / (gas * l * l); + let mu = 1.0 / (1.0 + phi); + let (mut k, mut m) = ([[0.0; 4]; 4], [[0.0; 4]; 4]); + for &(x, w) in gauss_legendre(4) { + let xi = 0.5 * (x + 1.0); + let w = 0.5 * w * l; + let x2 = xi * xi; + let x3 = x2 * xi; + let nw = [ + mu * (2.0 * x3 - 3.0 * x2 - phi * xi + 1.0 + phi), + mu * l * (x3 - (2.0 + phi / 2.0) * x2 + (1.0 + phi / 2.0) * xi), + -mu * (2.0 * x3 - 3.0 * x2 - phi * xi), + mu * l * (x3 - (1.0 - phi / 2.0) * x2 - phi / 2.0 * xi), + ]; + let nt = + [mu * 6.0 / l * (x2 - xi), mu * (3.0 * x2 - (4.0 + phi) * xi + 1.0 + phi), -mu * 6.0 / l * (x2 - xi), mu * (3.0 * x2 - (2.0 - phi) * xi)]; + // θ' and γ = w' − θ. + let dt = [ + mu * 6.0 / (l * l) * (2.0 * xi - 1.0), + mu / l * (6.0 * xi - 4.0 - phi), + -mu * 6.0 / (l * l) * (2.0 * xi - 1.0), + mu / l * (6.0 * xi - 2.0 + phi), + ]; + let dw = [ + mu / l * (6.0 * x2 - 6.0 * xi - phi), + mu * (3.0 * x2 - (4.0 + phi) * xi + 1.0 + phi / 2.0), + -mu / l * (6.0 * x2 - 6.0 * xi - phi), + mu * (3.0 * x2 - (2.0 - phi) * xi - phi / 2.0), + ]; + let g: [f64; 4] = std::array::from_fn(|a| dw[a] - nt[a]); + for a in 0..4 { + for b in 0..4 { + k[a][b] += w * (ei * dt[a] * dt[b] + gas * g[a] * g[b]); + m[a][b] += w * (rho_a * nw[a] * nw[b] + rho_i * nt[a] * nt[b]); + } + } + } + (k, m) +} + +/// A two-node beam element. +#[derive(Debug, Clone, Copy)] +pub struct Beam { + axes: Axes, + length: f64, + section: BeamSection, +} + +/// A line element's unit axis (node 1 to node 2) and length. +fn line_axis(xyz: &[[f64; 3]]) -> Result<(Vec3, f64)> { + let [a, b] = xyz else { + return Err(Error::Invalid(format!("a two-node line element has 2 nodes, got {}", xyz.len()))); + }; + modelcraft_model::check_coordinates(xyz.as_flattened()).map_err(|e| Error::Invalid(e.to_string()))?; + let d = sub(*b, *a); + let length = norm(d); + let ex = unit(d).ok_or_else(|| Error::Degenerate("a line element has zero length".into()))?; + Ok((ex, length)) +} + +/// The beam's local axes from its end points and orientation vector. +fn beam_axes(xyz: &[[f64; 3]], orientation: Vec3) -> Result<(Axes, f64)> { + let (ex, length) = line_axis(xyz)?; + let o = unit(orientation).ok_or_else(|| Error::Invalid("a beam's orientation vector must be nonzero and finite".into()))?; + let proj = sub(o, ex.map(|c| c * dot(o, ex))); + let ey = unit(proj) + .filter(|_| norm(proj) > 1e-6) + .ok_or_else(|| Error::Invalid("the beam section's orientation vector is parallel to the beam".into()))?; + Ok(([ex, ey, cross(ex, ey)], length)) +} + +impl Beam { + pub fn new(xyz: &[[f64; 3]], orientation: Vec3, section: BeamSection) -> Result { + let (axes, length) = beam_axes(xyz, orientation)?; + Ok(Beam { axes, length, section }) + } + + pub fn axes(&self) -> Axes { + self.axes + } + + pub fn length(&self) -> f64 { + self.length + } + + /// The element's unknowns: three translations and three rotations per node. + pub fn dofs(&self) -> ElementDofs { + ElementDofs::new(DofSet::STRUCTURAL) + } + + /// Local stiffness and consistent mass, 12 × 12. + fn local(&self, density: f64) -> (Matrix, Matrix) { + let (p, e, g, l) = (self.section.properties, self.section.young, self.section.shear_modulus, self.length); + let (mut k, mut m) = (Matrix::zeros(12), Matrix::zeros(12)); + let pair = |k: &mut Matrix, m: &mut Matrix, dofs: [usize; 2], stiff: f64, mass: f64| { + for (a, &i) in dofs.iter().enumerate() { + for (b, &j) in dofs.iter().enumerate() { + k.add(i, j, if a == b { stiff } else { -stiff }); + m.add(i, j, mass / 6.0 * if a == b { 2.0 } else { 1.0 }); + } + } + }; + pair(&mut k, &mut m, [0, 6], e * p.area / l, density * p.area * l); + pair(&mut k, &mut m, [3, 9], g * p.torsion / l, density * (p.iy + p.iz) * l); + // x–y plane: v and θ_z = dv/dx; x–z plane: w and θ_y = −dw/dx. + let planes = [([1, 5, 7, 11], [1.0; 4], p.iz, p.shear_y), ([2, 4, 8, 10], [1.0, -1.0, 1.0, -1.0], p.iy, p.shear_z)]; + for (dofs, sign, i, shear) in planes { + let (kp, mp) = bending_plane(e * i, g * shear, density * p.area, density * i, l); + for a in 0..4 { + for b in 0..4 { + k.add(dofs[a], dofs[b], sign[a] * sign[b] * kp[a][b]); + m.add(dofs[a], dofs[b], sign[a] * sign[b] * mp[a][b]); + } + } + } + (k, m) + } + + pub fn stiffness(&self) -> Matrix { + rotate_matrix(&self.local(0.0).0, &self.axes) + } + + pub fn mass(&self, density: f64) -> Matrix { + rotate_matrix(&self.local(density).1, &self.axes) + } + + /// The lumped mass, block diagonal (3 × 3 rotation blocks) in global axes. + pub fn lumped_mass(&self, density: f64) -> Matrix { + let (_, m) = self.local(density); + let (p, l) = (self.section.properties, self.length); + let half = density * p.area * l / 2.0; + let mut d = Matrix::zeros(12); + // Bending rotations scaled like their plane's translations (equal in both planes). + for (rotation, deflection) in [(5, 1), (4, 2)] { + let scale = 2.0 * half / (m.get(deflection, deflection) + m.get(deflection + 6, deflection + 6)); + d.add(rotation, rotation, scale * m.get(rotation, rotation)); + d.add(rotation + 6, rotation + 6, scale * m.get(rotation + 6, rotation + 6)); + } + for node in [0, 6] { + for c in 0..3 { + d.add(node + c, node + c, half); + } + d.add(node + 3, node + 3, density * (p.iy + p.iz) * l / 2.0); + } + rotate_matrix(&d, &self.axes) + } + + /// Consistent nodal loads of a uniform force per volume b, in local axes: the consistent + /// mass of unit density times the uniform nodal translation. + fn local_body_force(&self, b: [f64; 3]) -> Vec { + self.local(1.0).1.mul_vec(&translation(to_local(&self.axes, b), 2)) + } + + /// Consistent nodal forces and moments (global) of a uniform force per volume b. + pub fn body_force(&self, b: [f64; 3]) -> Vec { + rotate_to_global(&self.local_body_force(b), &self.axes) + } + + /// Section forces at the ends for the nodal displacements and rotations `u` (global, 12) + /// under a uniform force per volume `body` acting on the element. + pub fn forces(&self, u: &[f64], body: [f64; 3]) -> Result { + if u.len() != 12 { + return Err(Error::Invalid(format!("{} values for a beam element with 12 degrees of freedom", u.len()))); + } + let ku = self.local(0.0).0.mul_vec(&rotate_to_local(u, &self.axes)); + let end: Vec = ku.iter().zip(self.local_body_force(body)).map(|(a, b)| a - b).collect(); + let at = |i: usize| end.get(i).copied().unwrap_or(0.0); + Ok(BeamForces { axes: self.axes, ends: [std::array::from_fn(|c| -at(c)), std::array::from_fn(|c| at(6 + c))] }) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + const STEEL: Elastic = Elastic::Isotropic { young: 210e3, poisson: 0.3 }; + + /// Rigid translations and rotations about a point give no force. + fn rigid_motions(xyz: &[[f64; 3]]) -> Vec> { + (0..6) + .map(|m| { + xyz.iter() + .flat_map(|x| { + let mut d = [0.0; 6]; + if m < 3 { + d[m] = 1.0; + } else { + let mut w = [0.0; 3]; + w[m - 3] = 1.0; + let t = cross(w, *x); + d[..3].copy_from_slice(&t); + d[3..].copy_from_slice(&w); + } + d + }) + .collect() + }) + .collect() + } + + #[test] + fn one_beam_element_is_exact_for_a_tip_load() { + // Cantilever along a skew axis; tip force along the local z axis: δ = PL³/3EI + PL/(G A_s). + let profile = BeamProfile::Rectangle { width: 3.0, height: 5.0 }; + let section = BeamSection::new(&profile, &STEEL).unwrap(); + let xyz = [[1.0, 2.0, 3.0], [41.0, 32.0, 3.0]]; + let beam = Beam::new(&xyz, [0.0, 0.0, 1.0], section).unwrap(); + let l = beam.length(); + assert!((l - 50.0).abs() < 1e-12); + let k = beam.stiffness(); + for r in rigid_motions(&xyz) { + let f = k.mul_vec(&r); + assert!(f.iter().all(|v| v.abs() < 1e-6), "{f:?}"); + } + // Clamp node 1: solve the 6 × 6 block of node 2 for a unit force along local z. + let axes = beam.axes(); + let mut k22 = Matrix::zeros(6); + for i in 0..6 { + for j in 0..6 { + k22.add(i, j, k.get(6 + i, 6 + j)); + } + } + let mut f = vec![0.0; 6]; + f[..3].copy_from_slice(&axes[2]); + crate::dense::Cholesky::new(&k22).unwrap().solve(&mut f); + let p = section.properties; + let want = l.powi(3) / (3.0 * section.young * p.iy) + l / (section.shear_modulus * p.shear_z); + let got = dot([f[0], f[1], f[2]], axes[2]); + assert!((got - want).abs() < 1e-12 * want, "{got} vs {want}"); + // Section forces at the clamp: V_z = 1, M_y = −L (the moment about y of a +z load at +x). + let mut u = vec![0.0; 12]; + u[6..].copy_from_slice(&f); + let forces = beam.forces(&u, [0.0; 3]).unwrap(); + assert!((forces.ends[0][2] - 1.0).abs() < 1e-9 && (forces.ends[0][4] + l).abs() < 1e-8, "{forces:?}"); + assert!((forces.ends[1][2] - 1.0).abs() < 1e-9 && forces.ends[1][4].abs() < 1e-8, "{forces:?}"); + // Mass: the translational entries sum to the mass along each direction. + let m = beam.mass(7.8e-9); + let total = 7.8e-9 * p.area * l; + for d in 0..3 { + let s: f64 = [d, 6 + d].iter().flat_map(|&i| [d, 6 + d].map(|j| m.get(i, j))).sum(); + assert!((s - total).abs() < 1e-12 * total, "{s} {total}"); + } + let lumped = beam.lumped_mass(7.8e-9); + assert!((lumped.get(0, 0) - total / 2.0).abs() < 1e-12 * total); + assert!(Beam::new(&xyz, [40.0, 30.0, 0.0], section).is_err(), "orientation along the axis"); + assert!(Beam::new(&[xyz[0], xyz[0]], [0.0, 0.0, 1.0], section).is_err(), "zero length"); + } + + #[test] + fn profile_properties() { + // A square: J = 0.1406 a⁴ (Timoshenko & Goodier's table, 0.1406); thin strip J → b t³/3. + let sq = properties(&BeamProfile::Rectangle { width: 2.0, height: 2.0 }, 0.3); + assert!((sq.torsion / 16.0 - 0.1406).abs() < 1e-4, "{}", sq.torsion / 16.0); + let strip = properties(&BeamProfile::Rectangle { width: 1000.0, height: 1.0 }, 0.3); + assert!((strip.torsion / (1000.0 / 3.0) - 1.0).abs() < 1e-3, "{}", strip.torsion); + // A thin pipe: the shear coefficient tends to 2(1 + ν)/(4 + 3ν), J = 2 π r³ t. + let (r, t) = (100.0, 0.01); + let pipe = properties(&BeamProfile::Pipe { radius: r, thickness: t }, 0.3); + let k = pipe.shear_y / pipe.area; + assert!((k - 2.6 / 4.9).abs() < 1e-3, "{k}"); + assert!((pipe.torsion / (2.0 * std::f64::consts::PI * r.powi(3) * t) - 1.0).abs() < 1e-3); + // A box is a pair of I sections in bending: the same I_y with webs of twice the thickness. + let b = properties(&BeamProfile::Box { width: 10.0, height: 20.0, thickness: 1.0 }, 0.3); + let i = properties(&BeamProfile::I { height: 20.0, width: 10.0, web: 2.0, flange: 1.0 }, 0.3); + assert!((b.iy - i.iy).abs() < 1e-9 && (b.area - i.area).abs() < 1e-12); + let circle = properties(&BeamProfile::Circle { radius: 2.0 }, 0.3); + assert!((circle.torsion - std::f64::consts::PI * 8.0).abs() < 1e-12); + } +} diff --git a/crates/fem/src/frame.rs b/crates/fem/src/frame.rs new file mode 100644 index 0000000..e48c252 --- /dev/null +++ b/crates/fem/src/frame.rs @@ -0,0 +1,142 @@ +//! Local frames of structural elements (shells, beams, rods) and the rotation of their matrices +//! between the element's local axes and the global ones. +//! +//! A frame is a rotation R whose rows are the local axes e₁, e₂, e₃ in global components, so +//! local components are a_l = R a_g. Element degrees of freedom come in groups of three (the +//! translations, then the rotations, of each node), each rotated by R, so the element matrix is +//! K_g = Tᵀ K_l T with T = diag(R, R, …) (Cook, Malkus, Plesha & Witt 2002, §2.10 and §7.4). + +use crate::dense::Matrix; + +pub type Vec3 = [f64; 3]; + +/// The rows e₁, e₂, e₃ of a rotation from global to local components. +pub type Axes = [Vec3; 3]; + +pub fn sub(a: Vec3, b: Vec3) -> Vec3 { + [a[0] - b[0], a[1] - b[1], a[2] - b[2]] +} + +pub fn dot(a: Vec3, b: Vec3) -> f64 { + a[0] * b[0] + a[1] * b[1] + a[2] * b[2] +} + +pub fn cross(a: Vec3, b: Vec3) -> Vec3 { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] +} + +pub fn norm(a: Vec3) -> f64 { + dot(a, a).sqrt() +} + +/// The unit vector along `a`, or `None` when it is zero or not finite. +pub fn unit(a: Vec3) -> Option { + let n = norm(a); + (n > 0.0 && n.is_finite()).then(|| a.map(|c| c / n)) +} + +/// R a: global components to local. +pub fn to_local(r: &Axes, a: Vec3) -> Vec3 { + [dot(r[0], a), dot(r[1], a), dot(r[2], a)] +} + +/// Rᵀ a: local components to global. +pub fn to_global(r: &Axes, a: Vec3) -> Vec3 { + let mut g = [0.0; 3]; + for (k, axis) in r.iter().enumerate() { + for c in 0..3 { + g[c] += axis[c] * a[k]; + } + } + g +} + +/// Tᵀ K T for a symmetric matrix whose degrees of freedom come in groups of three local +/// components; the result is exactly symmetric. +pub fn rotate_matrix(k: &Matrix, r: &Axes) -> Matrix { + let n = k.n; + let groups = n / 3; + let mut out = Matrix::zeros(n); + for a in 0..groups { + for b in a..groups { + // Rᵀ K_ab R, mirrored to (b, a): K is symmetric. + let mut kr = [[0.0; 3]; 3]; + for i in 0..3 { + for j in 0..3 { + kr[i][j] = (0..3).map(|m| k.get(3 * a + i, 3 * b + m) * r[m][j]).sum(); + } + } + for i in 0..3 { + for j in 0..3 { + let v: f64 = (0..3).map(|m| r[m][i] * kr[m][j]).sum(); + out.add(3 * a + i, 3 * b + j, v); + if a != b { + out.add(3 * b + j, 3 * a + i, v); + } + } + } + } + } + crate::element::symmetrize(&mut out); + out +} + +/// T v: a vector of global components in groups of three to local components. +pub fn rotate_to_local(v: &[f64], r: &Axes) -> Vec { + v.as_chunks::<3>().0.iter().flat_map(|g| to_local(r, *g)).collect() +} + +/// Tᵀ v: local components in groups of three to global ones. +pub fn rotate_to_global(v: &[f64], r: &Axes) -> Vec { + v.as_chunks::<3>().0.iter().flat_map(|g| to_global(r, *g)).collect() +} + +/// A uniform translation `b` of `nodes` nodes with six degrees of freedom each (rotations +/// zero): the nodal field whose product with a consistent mass gives body-force loads. +pub fn translation(b: Vec3, nodes: usize) -> Vec { + (0..nodes).flat_map(|_| [b[0], b[1], b[2], 0.0, 0.0, 0.0]).collect() +} + +/// A symmetric tensor in local axes, [[xx, xy, xz], [xy, yy, yz], [xz, yz, zz]], in global +/// components as the 6-vector xx, yy, zz, xy, xz, yz: Rᵀ S R. +pub fn tensor_to_global(r: &Axes, s: [[f64; 3]; 3]) -> [f64; 6] { + let mut g = [[0.0; 3]; 3]; + for i in 0..3 { + for j in 0..3 { + g[i][j] = (0..3).map(|a| (0..3).map(|b| r[a][i] * s[a][b] * r[b][j]).sum::()).sum(); + } + } + [g[0][0], g[1][1], g[2][2], g[0][1], g[0][2], g[1][2]] +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn rotations_round_trip() { + let e1 = unit([1.0, 1.0, 0.0]).unwrap(); + let e3 = [0.0, 0.0, 1.0]; + let r: Axes = [e1, cross(e3, e1), e3]; + let v = [0.3, -1.2, 2.0]; + let back = to_global(&r, to_local(&r, v)); + assert!(back.iter().zip(v).all(|(a, b)| (a - b).abs() < 1e-15)); + // A uniaxial stress along e1 is xx = yy = xy = σ/2 globally. + let s = tensor_to_global(&r, [[2.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]); + assert!((s[0] - 1.0).abs() < 1e-15 && (s[1] - 1.0).abs() < 1e-15 && (s[3] - 1.0).abs() < 1e-15 && s[2] == 0.0, "{s:?}"); + // Rotating a matrix keeps its quadratic form: (T v)ᵀ K (T v) = vᵀ (Tᵀ K T) v. + let mut k = Matrix::zeros(6); + for i in 0..6 { + for j in 0..6 { + k.data[i * 6 + j] = 1.0 / (1.0 + i as f64 + j as f64); + } + } + let g = rotate_matrix(&k, &r); + let x = [0.1, 0.2, -0.3, 0.7, 0.0, 1.0]; + let xl = rotate_to_local(&x, &r); + let ql: f64 = (0..6).map(|i| (0..6).map(|j| xl[i] * k.get(i, j) * xl[j]).sum::()).sum(); + let qg: f64 = (0..6).map(|i| (0..6).map(|j| x[i] * g.get(i, j) * x[j]).sum::()).sum(); + assert!((ql - qg).abs() < 1e-14, "{ql} {qg}"); + assert!(unit([0.0; 3]).is_none() && unit([f64::NAN, 0.0, 0.0]).is_none()); + } +} diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index 42f8a1d..d5352e3 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -17,6 +17,10 @@ //! and lumped mass, body forces, integration-point stresses and their extrapolation to nodes; //! - [`plane`]: 2D elements (plane stress, plane strain, axisymmetric) on tri3, tri6, quad4 and //! quad8, with two displacement components per node; +//! - [`shell`]: flat MITC4 quadrilateral and MITC3 triangle shells with a drilling rotation, six +//! degrees of freedom per node; +//! - [`beam`]: two-node Timoshenko beams with profile section properties; +//! - [`frame`]: local axes of shells and beams and the rotation of their matrices; //! - [`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; @@ -29,7 +33,7 @@ //! - [`dense`]: the small dense linear algebra inside an element. //! //! Conventions: element matrices are ordered node by node, each node's degrees of freedom in -//! [`DofKind`] order (three translations for solids; docs/fem.md); stress +//! [`DofKind`] order (three translations for solids, six for shells and beams; docs/fem.md); stress //! and strain 6-vectors in the order xx, yy, zz, xy, xz, yz with engineering shear strains. #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] #![forbid(unsafe_code)] @@ -37,10 +41,12 @@ #![allow(clippy::needless_range_loop)] pub mod basis; +pub mod beam; pub mod dense; pub mod dof; pub mod element; pub mod face; +pub mod frame; pub mod jet; pub mod material; pub mod plane; @@ -48,6 +54,7 @@ pub mod plasticity; pub mod quadrature; pub mod rotation; pub mod shape; +pub mod shell; pub mod truss; pub mod uniform; diff --git a/crates/fem/src/material.rs b/crates/fem/src/material.rs index 97ff75d..b31c161 100644 --- a/crates/fem/src/material.rs +++ b/crates/fem/src/material.rs @@ -18,11 +18,27 @@ use crate::{Error, Result}; /// A 6 × 6 elasticity matrix, row-major. pub type Elasticity = [[f64; 6]; 6]; -/// The isotropic elasticity matrix. -pub fn isotropic(young: f64, poisson: f64) -> Result { +/// Check isotropic constants: E > 0 and −1 < ν < 0.5. +fn check_isotropic(young: f64, poisson: f64) -> Result<()> { if !(young > 0.0 && young.is_finite() && poisson > -1.0 && poisson < 0.5) { return Err(Error::Invalid(format!("isotropic constants E = {young}, ν = {poisson} are not a stable material"))); } + Ok(()) +} + +/// The checked (E, ν) of an isotropic material, for the structural sections (shells, beams, +/// rods) that need one; `what` names the section kind in the error for other materials. +pub fn isotropic_constants(elastic: &Elastic, what: &str) -> Result<(f64, f64)> { + let Elastic::Isotropic { young, poisson } = *elastic else { + return Err(Error::NotSupported(format!("orthotropic {what} sections (they need material orientations)"))); + }; + check_isotropic(young, poisson)?; + Ok((young, poisson)) +} + +/// The isotropic elasticity matrix. +pub fn isotropic(young: f64, poisson: f64) -> Result { + check_isotropic(young, poisson)?; let lambda = young * poisson / ((1.0 + poisson) * (1.0 - 2.0 * poisson)); let mu = young / (2.0 * (1.0 + poisson)); let mut d = [[0.0; 6]; 6]; diff --git a/crates/fem/src/shell.rs b/crates/fem/src/shell.rs new file mode 100644 index 0000000..9aa2b28 --- /dev/null +++ b/crates/fem/src/shell.rs @@ -0,0 +1,593 @@ +//! Flat shell elements: the four-node MITC4 quadrilateral and the three-node MITC3 triangle, +//! with six degrees of freedom per node (three translations, three rotations). +//! +//! Each element is flat: its nodes are projected onto the element's mean plane, where the +//! membrane, the Reissner–Mindlin plate and a drilling stiffness are superposed in local axes and +//! the result rotated to the global ones (the facet approach; Zienkiewicz, Taylor & Fox 2014, +//! *The Finite Element Method for Solid and Structural Mechanics*, 7th ed., ch. 9–10). Curved +//! shells are approximated by facets; quadrilaterals on a cylinder whose edges follow the +//! generators and the circumference are exactly flat. +//! +//! - **Kinematics.** With the reference surface at the nodes and z along the element normal, +//! u = u₀ + z θ_y, v = v₀ − z θ_x, w = w₀: membrane strains ε = (u₀,x, v₀,y, u₀,y + v₀,x), +//! curvatures κ = (θ_y,x, −θ_x,y, θ_y,y − θ_x,x) and transverse shear strains +//! γ = (w,x + θ_y, w,y − θ_x). +//! - **Section.** Linear elastic and homogeneous through the thickness, so the integration +//! through the thickness is done exactly in closed form (any number of section points gives the +//! same stiffness): with the mid-surface at z = e from the reference surface (e = −offset · t), +//! N = A ε + B κ, M = B ε + D κ with A = Q t, B = Q t e, D = Q (t³/12 + t e²) for the plane +//! stress matrix Q, and Q_s = k G t γ with the shear correction k = 5/6. +//! - **Transverse shear, quadrilateral (MITC4).** The covariant transverse shear strains are +//! sampled at the mid-sides and interpolated linearly across the element: ẽ_rt from the points +//! (0, ±1), ẽ_st from (±1, 0), then turned into Cartesian strains with the Jacobian +//! (E. N. Dvorkin, K.-J. Bathe, *A continuum mechanics based four-node shell element for +//! general nonlinear analysis*, Engineering Computations 1(1):77–88, 1984; in the flat-plate +//! form of K.-J. Bathe, E. N. Dvorkin, *A four-node plate bending element based on +//! Mindlin/Reissner plate theory and a mixed interpolation*, IJNME 21(2):367–383, 1985). This +//! removes shear locking while passing the patch tests; membrane and bending use the bilinear +//! fields, all integrated with 2 × 2 Gauss points. +//! - **Transverse shear, triangle (MITC3).** The covariant shear strains are tied at the three +//! mid-sides: ẽ_rt = e_rt⁽¹⁾ + c s, ẽ_st = e_st⁽²⁾ − c r with c = e_rt⁽³⁾ − e_rt⁽¹⁾ − e_st⁽³⁾ + +//! e_st⁽²⁾, tying points (½, 0), (0, ½) and (½, ½), which keeps the tangential shear strain of +//! each edge constant along it (P.-S. Lee, K.-J. Bathe, *Development of MITC isotropic +//! triangular shell finite elements*, Computers & Structures 82(11–13):945–962, 2004). The +//! membrane is the constant-strain triangle. +//! - **Drilling rotation.** The rotation about the normal has no stiffness in shell theory; it +//! is tied to the membrane's in-plane rotation ω = ½ (v,x − u,y) by the penalty +//! ∫ γ t (θ_z − ω)² dA of the variational formulation of T. J. R. Hughes, F. Brezzi, *On drilling +//! degrees of freedom*, CMAME 72(1):105–121, 1989, with γ = G/100. Hughes and Brezzi take +//! γ = G with a drilling-enhanced (Allman-type) membrane; on the plain bilinear membrane the +//! in-plane rotation of a parallelogram is linear, which θ_z interpolates exactly, so the penalty +//! does not lock, and the smaller γ keeps distorted elements free of it too. Being a function of +//! θ_z − ω it leaves rigid motions free. Scaled with G t like the membrane (not with the bending +//! rigidity), it ties each node's rotation about the surface normal to the membrane as the mesh +//! is refined, so curved shells modelled with facets converge steadily: the Scordelis–Lo and +//! pinched-cylinder benchmarks change by less than 0.5 % for γ from G/1000 to G, whereas a +//! penalty that is a fixed fraction of the bending stiffness lets that rotation go soft on fine +//! meshes. +//! - **Mass.** Consistent: ρ t on the translations, ρ (t³/12 + t e²) on the rotations and +//! ρ t e coupling translations and bending rotations when the reference surface is offset. The +//! drilling rotation gets the same rotary inertia as the bending ones (shell theory gives it +//! none; it is tiny next to the in-plane inertia of the translations and keeps the mass matrix +//! positive definite, as the frequency solver's Lanczos iteration needs). Lumped: the diagonal +//! scaled to the element's mass (Hinton, Rock & Zienkiewicz 1976), the rotary inertia +//! isotropic about the node so the lumped matrix stays diagonal in global axes. +//! - **Stresses.** Section forces and the stresses at the bottom, middle and top of the section +//! at the integration points, extrapolated bilinearly to the nodes (quadrilateral) or taken at +//! the nodes (triangle, constant membrane and bending fields). Transverse shear stresses are +//! reported at the mid-surface as the peak of the parabolic distribution, 3/2 Q/t. + +use modelcraft_model::{Elastic, ElementKind}; + +use crate::dense::Matrix; +use crate::dof::{DofSet, ElementDofs}; +use crate::frame::{Axes, cross, dot, norm, rotate_matrix, rotate_to_global, rotate_to_local, sub, tensor_to_global, to_local, translation, unit}; +use crate::quadrature::Rule; +use crate::{Error, Result}; + +/// Transverse shear correction factor of a homogeneous section. +pub const SHEAR_CORRECTION: f64 = 5.0 / 6.0; + +/// The drilling penalty γ as a fraction of the shear modulus. +pub const DRILLING: f64 = 1e-2; + +/// Degrees of freedom per shell node: translations along local (or global) x, y, z, then +/// rotations about them. +pub const NODE_DOFS: usize = 6; + +/// A shell section with an isotropic material. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct ShellSection { + pub thickness: f64, + /// The mid-surface's position along the element normal from the nodes (the reference + /// surface): −offset · thickness. + pub eccentricity: f64, + pub young: f64, + pub poisson: f64, +} + +impl ShellSection { + /// A section of the model's thickness and offset (a fraction of the thickness, 0.5 puts the + /// nodes on the side the normal points to). + pub fn new(thickness: f64, offset: f64, elastic: &Elastic) -> Result { + let (young, poisson) = crate::material::isotropic_constants(elastic, "shell")?; + if !(thickness > 0.0 && thickness.is_finite() && offset.abs() <= 0.5) { + return Err(Error::Invalid(format!( + "a shell section needs a positive thickness and an offset in −0.5..=0.5, got {thickness} and {offset}" + ))); + } + Ok(ShellSection { thickness, eccentricity: -offset * thickness, young, poisson }) + } + + /// The plane stress matrix Q. + fn plane_stress(&self) -> [[f64; 3]; 3] { + let (e, nu) = (self.young, self.poisson); + let c = e / (1.0 - nu * nu); + [[c, c * nu, 0.0], [c * nu, c, 0.0], [0.0, 0.0, c * (1.0 - nu) / 2.0]] + } + + /// The membrane–bending stiffness [[A, B], [B, D]], 6 × 6. + fn abd(&self) -> [[f64; 6]; 6] { + let q = self.plane_stress(); + let (t, e) = (self.thickness, self.eccentricity); + let (a, b, d) = (t, t * e, t * t * t / 12.0 + t * e * e); + let mut out = [[0.0; 6]; 6]; + for i in 0..3 { + for j in 0..3 { + out[i][j] = a * q[i][j]; + out[i][j + 3] = b * q[i][j]; + out[i + 3][j] = b * q[i][j]; + out[i + 3][j + 3] = d * q[i][j]; + } + } + out + } + + pub fn shear_modulus(&self) -> f64 { + self.young / (2.0 * (1.0 + self.poisson)) + } + + /// The transverse shear stiffness k G t. + fn shear(&self) -> f64 { + SHEAR_CORRECTION * self.shear_modulus() * self.thickness + } +} + +/// Section forces and stresses of a shell element, per node in the element's node order. +#[derive(Debug, Clone, PartialEq)] +pub struct ShellStress { + /// The element's local axes: e₁, e₂ in its plane, e₃ its normal (the right-hand normal of + /// the node order). + pub axes: Axes, + /// N_x, N_y, N_xy, M_x, M_y, M_xy, Q_x, Q_y in the local axes, per unit length, about the + /// reference surface. + pub forces: Vec<[f64; 8]>, + /// Stress (global xx, yy, zz, xy, xz, yz) at the bottom (−e₃ side), the middle and the top + /// (+e₃ side) of the section. + pub bottom: Vec<[f64; 6]>, + pub middle: Vec<[f64; 6]>, + pub top: Vec<[f64; 6]>, +} + +/// The shape functions of the element's corners at (r, s) and their derivatives along r and s. +fn shape(kind: ElementKind, r: f64, s: f64) -> ([f64; 4], [[f64; 2]; 4]) { + if kind == ElementKind::Tri3 { + ([1.0 - r - s, r, s, 0.0], [[-1.0, -1.0], [1.0, 0.0], [0.0, 1.0], [0.0, 0.0]]) + } else { + let c = [[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]]; + let mut n = [0.0; 4]; + let mut d = [[0.0; 2]; 4]; + for (i, ci) in c.iter().enumerate() { + n[i] = 0.25 * (1.0 + ci[0] * r) * (1.0 + ci[1] * s); + d[i] = [0.25 * ci[0] * (1.0 + ci[1] * s), 0.25 * ci[1] * (1.0 + ci[0] * r)]; + } + (n, d) + } +} + +/// The map's Jacobian at a point: J = [[x,r, y,r], [x,s, y,s]], its inverse and determinant. +struct Map { + n: [f64; 4], + /// ∂Nᵢ/∂r, ∂Nᵢ/∂s. + dn: [[f64; 2]; 4], + /// ∂Nᵢ/∂x, ∂Nᵢ/∂y. + dx: [[f64; 2]; 4], + j: [[f64; 2]; 2], + inv: [[f64; 2]; 2], + det: f64, +} + +/// Most degrees of freedom of a shell element (quad4). +const MAX_DOFS: usize = 4 * NODE_DOFS; + +/// A row of a strain operator over the element's local degrees of freedom (zero past them). +type Row = [f64; MAX_DOFS]; + +/// The covariant shear rows at the element's tying points, which depend on its geometry only: +/// quad4 e_rt at (0, 1), (0, −1) and e_st at (1, 0), (−1, 0); tri3 e_rt⁽¹⁾, e_st⁽²⁾ and the +/// combination c. +struct Tying([Row; 4]); + +/// A flat shell element. +#[derive(Debug, Clone)] +pub struct Shell { + kind: ElementKind, + nodes: usize, + axes: Axes, + /// The nodes in the element plane (local x, y). + local: [[f64; 2]; 4], + section: ShellSection, + area: f64, +} + +impl Shell { + pub fn new(kind: ElementKind, xyz: &[[f64; 3]], section: ShellSection) -> Result { + let nodes = match kind { + ElementKind::Tri3 => 3, + ElementKind::Quad4 => 4, + ElementKind::Tri6 | ElementKind::Quad8 => { + return Err(Error::NotSupported(format!("{} shells (the solver has quad4 and tri3 shells so far)", kind.name()))); + } + _ => return Err(Error::Invalid(format!("{} is not a shell element", kind.name()))), + }; + if xyz.len() != nodes { + return Err(Error::Invalid(format!("a {} element has {nodes} nodes, got {}", kind.name(), xyz.len()))); + } + modelcraft_model::check_coordinates(xyz.as_flattened()).map_err(|e| Error::Invalid(e.to_string()))?; + let degenerate = || Error::Degenerate(format!("a {} shell element is degenerate (zero area, inverted or not convex)", kind.name())); + let p = |i: usize| xyz.get(i).copied().unwrap_or_default(); + let (normal, along) = if nodes == 3 { + (cross(sub(p(1), p(0)), sub(p(2), p(0))), sub(p(1), p(0))) + } else { + let along = sub([0, 1, 2].map(|k| p(1)[k] + p(2)[k]), [0, 1, 2].map(|k| p(0)[k] + p(3)[k])); + (cross(sub(p(2), p(0)), sub(p(3), p(1))), along) + }; + let e3 = unit(normal).ok_or_else(degenerate)?; + let e1 = unit(sub(along, e3.map(|c| c * dot(along, e3)))).ok_or_else(degenerate)?; + let e2 = cross(e3, e1); + let axes = [e1, e2, e3]; + let mut centre = [0.0; 3]; + for x in xyz { + for k in 0..3 { + centre[k] += x[k] / nodes as f64; + } + } + let mut local = [[0.0; 2]; 4]; + let mut size = 0.0f64; + for (l, x) in local.iter_mut().zip(xyz) { + let d = sub(*x, centre); + *l = [dot(d, e1), dot(d, e2)]; + size = size.max(norm(d)); + } + // Corners turn left (counter-clockwise about e₃) by a margin: positive area, convex. + for i in 0..nodes { + let (a, b, c) = (local[(i + nodes - 1) % nodes], local[i], local[(i + 1) % nodes]); + let turn = (b[0] - a[0]) * (c[1] - b[1]) - (b[1] - a[1]) * (c[0] - b[0]); + if !(turn > 1e-10 * size * size) { + return Err(degenerate()); + } + } + // A strongly warped quadrilateral is not a facet. + let warp = xyz.iter().map(|x| dot(sub(*x, centre), e3).abs()).fold(0.0, f64::max); + if warp > 0.1 * size { + return Err(Error::Degenerate(format!( + "a quad4 shell element is warped by {warp:.3e} out of its plane (more than a tenth of its size); split it into triangles" + ))); + } + let mut shell = Shell { kind, nodes, axes, local, section, area: 0.0 }; + shell.area = shell.rule().points().iter().map(|q| shell.map(q.xi[0], q.xi[1]).map_or(0.0, |m| m.det * q.weight)).sum(); + if !(shell.area > 0.0) { + return Err(degenerate()); + } + Ok(shell) + } + + pub fn axes(&self) -> Axes { + self.axes + } + + pub fn area(&self) -> f64 { + self.area + } + + /// The element's unknowns: three translations and three rotations per node. + pub fn dofs(&self) -> ElementDofs { + ElementDofs::new(DofSet::STRUCTURAL) + } + + /// The size of the element's matrices. + fn size(&self) -> usize { + NODE_DOFS * self.nodes + } + + fn rule(&self) -> Rule { + if self.kind == ElementKind::Tri3 { Rule::TriConical(2) } else { Rule::Quad(2) } + } + + fn map(&self, r: f64, s: f64) -> Result { + let (n, d) = shape(self.kind, r, s); + let mut j = [[0.0; 2]; 2]; + for (di, x) in d.iter().zip(&self.local).take(self.nodes) { + for a in 0..2 { + for b in 0..2 { + j[a][b] += di[a] * x[b]; + } + } + } + let det = j[0][0] * j[1][1] - j[0][1] * j[1][0]; + if !(det > 0.0 && det.is_finite()) { + return Err(Error::Degenerate(format!("a {} shell element is inverted or degenerate", self.kind.name()))); + } + let inv = [[j[1][1] / det, -j[0][1] / det], [-j[1][0] / det, j[0][0] / det]]; + let mut dx = [[0.0; 2]; 4]; + for (g, di) in dx.iter_mut().zip(&d) { + g[0] = inv[0][0] * di[0] + inv[0][1] * di[1]; + g[1] = inv[1][0] * di[0] + inv[1][1] * di[1]; + } + Ok(Map { n, dn: d, dx, j, inv, det }) + } + + /// The covariant transverse shear strain e_rt (`along` 0) or e_st (`along` 1) of the + /// displacement interpolation at (r, s): w,r + θ_y x,r − θ_x y,r (and with s). + fn covariant(&self, r: f64, s: f64, along: usize) -> Result { + let m = self.map(r, s)?; + let (xr, yr) = (m.j[along][0], m.j[along][1]); + let mut row = [0.0; MAX_DOFS]; + for i in 0..self.nodes { + row[6 * i + 2] = m.dn[i][along]; + row[6 * i + 4] = m.n[i] * xr; + row[6 * i + 3] = -m.n[i] * yr; + } + Ok(row) + } + + fn tying(&self) -> Result { + if self.kind == ElementKind::Tri3 { + let e1 = self.covariant(0.5, 0.0, 0)?; + let e2 = self.covariant(0.0, 0.5, 1)?; + let e3r = self.covariant(0.5, 0.5, 0)?; + let e3s = self.covariant(0.5, 0.5, 1)?; + let c: Row = std::array::from_fn(|k| e3r[k] - e1[k] - e3s[k] + e2[k]); + Ok(Tying([e1, e2, c, [0.0; MAX_DOFS]])) + } else { + Ok(Tying([self.covariant(0.0, 1.0, 0)?, self.covariant(0.0, -1.0, 0)?, self.covariant(1.0, 0.0, 1)?, self.covariant(-1.0, 0.0, 1)?])) + } + } + + /// The assumed covariant shear strains (ẽ_rt, ẽ_st) at (r, s). + fn assumed_shear(&self, tying: &Tying, r: f64, s: f64) -> [Row; 2] { + let lerp = |a: &Row, wa: f64, b: &Row, wb: f64| -> Row { std::array::from_fn(|k| wa * a[k] + wb * b[k]) }; + let [t0, t1, t2, t3] = &tying.0; + if self.kind == ElementKind::Tri3 { + [lerp(t0, 1.0, t2, s), lerp(t1, 1.0, t2, -r)] + } else { + [lerp(t0, 0.5 * (1.0 + s), t1, 0.5 * (1.0 - s)), lerp(t2, 0.5 * (1.0 + r), t3, 0.5 * (1.0 - r))] + } + } + + /// The generalized strain operator at (r, s): membrane strains, curvatures, transverse shear + /// strains (8 rows) and the drilling row θ_z − ω, with the map. + fn operator(&self, tying: &Tying, r: f64, s: f64) -> Result<([Row; 9], Map)> { + let m = self.map(r, s)?; + let mut b = [[0.0; MAX_DOFS]; 9]; + for i in 0..self.nodes { + let [gx, gy] = m.dx[i]; + let (u, v, tx, ty, tz) = (6 * i, 6 * i + 1, 6 * i + 3, 6 * i + 4, 6 * i + 5); + b[0][u] = gx; + b[1][v] = gy; + b[2][u] = gy; + b[2][v] = gx; + b[3][ty] = gx; + b[4][tx] = -gy; + b[5][ty] = gy; + b[5][tx] = -gx; + b[8][tz] = m.n[i]; + b[8][u] = 0.5 * gy; + b[8][v] = -0.5 * gx; + } + let [ert, est] = self.assumed_shear(tying, r, s); + for k in 0..MAX_DOFS { + b[6][k] = m.inv[0][0] * ert[k] + m.inv[0][1] * est[k]; + b[7][k] = m.inv[1][0] * ert[k] + m.inv[1][1] * est[k]; + } + Ok((b, m)) + } + + /// The stiffness matrix in global axes, 6n × 6n, degrees of freedom node by node. + pub fn stiffness(&self) -> Result { + let n = self.size(); + let abd = self.section.abd(); + let ks = self.section.shear(); + let drill = DRILLING * self.section.shear_modulus() * self.section.thickness; + let tying = self.tying()?; + let mut k = Matrix::zeros(n); + for p in self.rule().points() { + let (b, m) = self.operator(&tying, p.xi[0], p.xi[1])?; + let w = m.det * p.weight; + // C B for the membrane–bending rows. + let cb: [Row; 6] = std::array::from_fn(|a| std::array::from_fn(|c| (0..6).map(|d| abd[a][d] * b[d][c]).sum())); + for i in 0..n { + for j in i..n { + let mut v: f64 = (0..6).map(|a| b[a][i] * cb[a][j]).sum(); + v += ks * (b[6][i] * b[6][j] + b[7][i] * b[7][j]); + v += drill * b[8][i] * b[8][j]; + if v != 0.0 { + k.add(i, j, w * v); + if i != j { + k.add(j, i, w * v); + } + } + } + } + } + Ok(rotate_matrix(&k, &self.axes)) + } + + /// The consistent mass matrix in global axes, 6n × 6n. + pub fn mass(&self, density: f64) -> Result { + Ok(rotate_matrix(&self.local_mass(density)?, &self.axes)) + } + + /// The consistent mass matrix in the element's axes. + fn local_mass(&self, density: f64) -> Result { + let (t, e) = (self.section.thickness, self.section.eccentricity); + let (m0, m1, m2) = (density * t, density * t * e, density * (t * t * t / 12.0 + t * e * e)); + let mut mm = Matrix::zeros(self.size()); + for p in self.rule().points() { + let m = self.map(p.xi[0], p.xi[1])?; + let w = m.det * p.weight; + for i in 0..self.nodes { + for j in 0..self.nodes { + let nn = w * m.n[i] * m.n[j]; + for c in 0..3 { + mm.add(6 * i + c, 6 * j + c, m0 * nn); + } + for c in 3..6 { + mm.add(6 * i + c, 6 * j + c, m2 * nn); + } + if m1 != 0.0 { + // u θ_y and −v θ_x, both ways. + mm.add(6 * i, 6 * j + 4, m1 * nn); + mm.add(6 * i + 4, 6 * j, m1 * nn); + mm.add(6 * i + 1, 6 * j + 3, -m1 * nn); + mm.add(6 * i + 3, 6 * j + 1, -m1 * nn); + } + } + } + } + Ok(mm) + } + + /// The lumped mass per degree of freedom (the same in local and global axes). + pub fn lumped_mass(&self, density: f64) -> Result> { + let (t, e) = (self.section.thickness, self.section.eccentricity); + let mut diag = vec![0.0; self.size()]; + let mut sum = 0.0; + for p in self.rule().points() { + let m = self.map(p.xi[0], p.xi[1])?; + for i in 0..self.nodes { + let v = m.det * p.weight * m.n[i] * m.n[i]; + diag[6 * i] += v; + sum += v; + } + } + if !(sum > 0.0) { + return Err(Error::Degenerate("a shell element has no area".into())); + } + let scale = self.area / sum; + let (m0, m2) = (density * t, density * (t * t * t / 12.0 + t * e * e)); + for i in 0..self.nodes { + let a = diag[6 * i] * scale; + for c in 0..3 { + diag[6 * i + c] = m0 * a; + diag[6 * i + 3 + c] = m2 * a; + } + } + Ok(diag) + } + + /// Consistent nodal forces and moments (global) of a uniform force per volume b: the + /// consistent mass of unit density times the uniform nodal translation b. + pub fn body_force(&self, b: [f64; 3]) -> Result> { + let r = translation(to_local(&self.axes, b), self.nodes); + Ok(rotate_to_global(&self.local_mass(1.0)?.mul_vec(&r), &self.axes)) + } + + /// Section forces and surface stresses for the element's nodal displacements and rotations + /// (global, 6 per node). + pub fn stresses(&self, u: &[f64]) -> Result { + if u.len() != self.size() { + return Err(Error::Invalid(format!("{} values for a shell element with {} degrees of freedom", u.len(), self.size()))); + } + let ul = rotate_to_local(u, &self.axes); + let tying = self.tying()?; + let strains = |r: f64, s: f64| -> Result<[f64; 8]> { + let (b, _) = self.operator(&tying, r, s)?; + Ok(std::array::from_fn(|a| b[a].iter().zip(&ul).map(|(x, y)| x * y).sum())) + }; + let nodal: Vec<[f64; 8]> = if self.kind == ElementKind::Tri3 { + [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]].iter().map(|c| strains(c[0], c[1])).collect::>()? + } else { + // 2 × 2 Gauss values, extrapolated with the bilinear functions of the points. + let g = 1.0 / 3f64.sqrt(); + let corners = [[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]]; + let at: Vec<[f64; 8]> = corners.iter().map(|c| strains(c[0] * g, c[1] * g)).collect::>()?; + corners + .iter() + .map(|c| { + let (n, _) = shape(ElementKind::Quad4, c[0] / g, c[1] / g); + std::array::from_fn(|k| n.iter().zip(&at).map(|(w, v)| w * v[k]).sum()) + }) + .collect() + }; + let abd = self.section.abd(); + let q = self.section.plane_stress(); + let ks = self.section.shear(); + let (t, e) = (self.section.thickness, self.section.eccentricity); + let mut out = ShellStress { axes: self.axes, forces: Vec::new(), bottom: Vec::new(), middle: Vec::new(), top: Vec::new() }; + for g in &nodal { + let mut f = [0.0; 8]; + for a in 0..6 { + f[a] = (0..6).map(|b| abd[a][b] * g[b]).sum(); + } + f[6] = ks * g[6]; + f[7] = ks * g[7]; + out.forces.push(f); + let at = |z: f64, shear: f64| { + let strain = [g[0] + z * g[3], g[1] + z * g[4], g[2] + z * g[5]]; + let s: [f64; 3] = std::array::from_fn(|a| (0..3).map(|b| q[a][b] * strain[b]).sum()); + let (txz, tyz) = (shear * f[6] / t, shear * f[7] / t); + tensor_to_global(&self.axes, [[s[0], s[2], txz], [s[2], s[1], tyz], [txz, tyz, 0.0]]) + }; + out.bottom.push(at(e - t / 2.0, 0.0)); + out.middle.push(at(e, 1.5)); + out.top.push(at(e + t / 2.0, 0.0)); + } + Ok(out) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + const STEEL: Elastic = Elastic::Isotropic { young: 210e3, poisson: 0.3 }; + + fn rigid(xyz: &[[f64; 3]], m: usize) -> Vec { + xyz.iter() + .flat_map(|x| { + let mut d = [0.0; 6]; + if m < 3 { + d[m] = 1.0; + } else { + let mut w = [0.0; 3]; + w[m - 3] = 1.0; + d[..3].copy_from_slice(&cross(w, *x)); + d[3..].copy_from_slice(&w); + } + d + }) + .collect() + } + + #[test] + fn shells_are_free_of_rigid_motion_and_symmetric() { + let section = ShellSection::new(0.5, 0.25, &STEEL).unwrap(); + // A skew quadrilateral and a triangle in a tilted plane. + let tilt = |x: f64, y: f64| [x, 0.6 * y, 0.8 * y + 1.0]; + let quad = [tilt(0.0, 0.0), tilt(2.0, 0.2), tilt(2.3, 1.7), tilt(-0.2, 1.4)]; + let tri = [tilt(0.0, 0.0), tilt(2.0, 0.2), tilt(0.4, 1.5)]; + for (kind, xyz) in [(ElementKind::Quad4, &quad[..]), (ElementKind::Tri3, &tri[..])] { + let shell = Shell::new(kind, xyz, section).unwrap(); + let k = shell.stiffness().unwrap(); + let scale = (0..k.n).map(|i| k.get(i, i)).fold(0.0, f64::max); + for i in 0..k.n { + for j in 0..k.n { + assert_eq!(k.get(i, j), k.get(j, i)); + } + } + for m in 0..6 { + let f = k.mul_vec(&rigid(xyz, m)); + assert!(f.iter().all(|v| v.abs() < 1e-10 * scale), "{kind:?} rigid {m}: {f:?}"); + } + // Mass: translations sum to ρ t A per direction (the offset couples rotations only). + let m = shell.mass(2.0).unwrap(); + let total = 2.0 * 0.5 * shell.area(); + let s: f64 = + (0..k.n).filter(|i| i % 6 == 0).flat_map(|i| (0..k.n).filter(|j| j % 6 == 0).map(move |j| (i, j))).map(|(i, j)| m.get(i, j)).sum(); + assert!((s - total).abs() < 1e-12 * total, "{s} {total}"); + let lumped = shell.lumped_mass(2.0).unwrap(); + let s: f64 = lumped.iter().step_by(6).sum(); + assert!((s - total).abs() < 1e-12 * total); + let body = shell.body_force([0.0, 0.0, -1.0]).unwrap(); + let fz: f64 = body.iter().skip(2).step_by(6).sum(); + assert!((fz + 0.5 * shell.area()).abs() < 1e-12, "{fz}"); + } + let flat = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.5, 0.2, 0.0]]; + assert!(Shell::new(ElementKind::Quad4, &flat, section).is_err(), "not convex"); + let line = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]]; + assert!(Shell::new(ElementKind::Tri3, &line, section).is_err()); + assert!(Shell::new(ElementKind::Quad8, &flat, section).is_err()); + assert!(ShellSection::new(1.0, 0.0, &Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0; 3] }).is_err()); + } +} diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 91e7f46..6359c71 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -1,8 +1,10 @@ //! Linear static and frequency analysis of a model's step, and their results. use modelcraft_fem::ElementType; +use modelcraft_fem::beam::BeamForces; use modelcraft_fem::element::{Extrapolation, extrapolation}; use modelcraft_fem::quadrature::Rule; +use modelcraft_fem::shell::ShellStress; use modelcraft_fem::{DofKind, DofLayout, PointStress, von_mises}; use modelcraft_model::{ElementKind, Model, Procedure, Resolved, Step}; use serde::Serialize; @@ -12,7 +14,7 @@ use crate::constraints::{Budget, Reduction, coupling_rows}; use crate::eigen::{Operator, Wanted, lanczos}; use crate::ordering::nested_dissection; use crate::sparse::{Csr, NodeGraph, add_block, add_diagonal, add_scalar}; -use crate::system::{Mesh, inelastic_warnings, loads, prescribed, translation_block}; +use crate::system::{Mesh, Structural, StructuralKernel, inelastic_warnings, loads, prescribed, structural_body_forces, translation_block}; use crate::{Error, Limits, Result}; /// Size and effort of a solution. @@ -29,6 +31,22 @@ pub struct Stats { pub lanczos_vectors: usize, } +/// Section forces and surface stresses of one shell element. +#[derive(Debug, Clone, PartialEq)] +pub struct ShellResult { + pub label: u32, + pub kind: ElementKind, + pub stress: ShellStress, +} + +/// End section forces of one beam element. +#[derive(Debug, Clone, PartialEq)] +pub struct BeamResult { + pub label: u32, + pub kind: ElementKind, + pub forces: BeamForces, +} + /// Stresses of one element. #[derive(Debug, Clone, PartialEq)] pub struct ElementStress { @@ -57,13 +75,23 @@ pub struct StaticResult { /// Node labels with degrees of freedom, sorted; the per-node vectors follow this order. pub nodes: Vec, pub displacement: Vec<[f64; 3]>, + /// Rotations about X, Y, Z (zero at nodes without rotations: those of solids and trusses + /// only). + pub rotation: Vec<[f64; 3]>, /// Reaction forces (zero where nothing is prescribed). pub reaction: Vec<[f64; 3]>, + /// Reaction moments at held rotations. + pub reaction_moment: Vec<[f64; 3]>, /// Stress averaged at the nodes over the elements that share them (xx, yy, zz, xy, xz, yz). pub stress: Vec<[f64; 6]>, /// Von Mises stress of the averaged nodal stress. pub von_mises: Vec, pub elements: Vec, + /// Shell section forces and surface stresses, by label (their top surface, the side their + /// normal points to, joins the nodal stress average). + pub shells: Vec, + /// Beam end section forces, by label. + pub beams: Vec, /// The coupled reference points, by name. pub reference_points: Vec, /// Sum of the applied forces and of the reaction forces (moments are not included). @@ -82,6 +110,8 @@ pub struct Mode { pub frequency: f64, /// Mass-normalized (φᵀ M φ = 1), per node in [`FrequencyResult::nodes`] order. pub shape: Vec<[f64; 3]>, + /// The rotations of the shape (zero at nodes without rotations). + pub rotation: Vec<[f64; 3]>, /// Participation factors Γ = φᵀ M r for unit translations along X, Y, Z. pub participation: [f64; 3], /// Effective mass Γ² along X, Y, Z. @@ -161,6 +191,15 @@ fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option) -> Result(full: &[f64], layout: &DofLayout, slot: usize) -> std::array::from_fn(|c| DofKind::from_direction(c).and_then(|k| layout.dof(slot, k)).and_then(|d| full.get(d)).copied().unwrap_or(0.0)) } -/// The mesh nodes' translations in a full vector, per node, and the reference points' -/// translations and rotations, per point. -fn split(full: &[f64], mesh: &Mesh) -> (Vec<[f64; 3]>, Vec<[f64; 6]>) { +/// A full vector split into the mesh nodes' translations and rotations (zero where a node has +/// none), per node, and the reference points' translations and rotations, per point. +fn split(full: &[f64], mesh: &Mesh) -> (Vec<[f64; 3]>, Vec<[f64; 3]>, Vec<[f64; 6]>) { let l = &mesh.layout; - ((0..l.nodes()).map(|n| directions(full, l, n)).collect(), (0..mesh.points.len()).map(|k| directions(full, l, l.point_slot(k))).collect()) + let (translation, rotation) = (0..l.nodes()) + .map(|n| { + let [a, b, c, d, e, f] = directions::<6>(full, l, n); + ([a, b, c], [d, e, f]) + }) + .unzip(); + (translation, rotation, (0..mesh.points.len()).map(|k| directions(full, l, l.point_slot(k))).collect()) } /// The matrix pattern of the mesh (node blocks of its elements, trusses and discrete elements) @@ -252,7 +297,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul let resolved = model.resolve()?; 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 mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).chain(mesh.warnings.iter().cloned()).collect(); warnings.extend(inelastic_warnings(model)); let (graph, order) = pattern(&mesh); ignored_dashpots(&mesh, &mut warnings); @@ -272,9 +317,9 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul if u.iter().any(|x| !x.is_finite()) { return Err(Error::Singular("the solution is not finite; the model is ill-conditioned".into())); } - let (displacement, point_motion) = split(&u, &mesh); + let (displacement, rotation, point_motion) = split(&u, &mesh); let reactions = red.reactions(&k, &u, &f); - let (reaction, point_reaction) = split(&reactions, &mesh); + let (reaction, reaction_moment, point_reaction) = split(&reactions, &mesh); let (mut applied_total, mut reaction_total) = ([0.0; 3], [0.0; 3]); for (dof, (fv, rv)) in f.iter().zip(&reactions).enumerate() { let Some(d) = mesh.kind(dof).and_then(DofKind::direction) else { continue }; @@ -337,6 +382,31 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul elements.push(ElementStress { label: t.label, kind: t.kind, points, nodal }); } elements.sort_by_key(|e| e.label); + let (mut shells, mut beams) = (Vec::new(), Vec::new()); + let per_volume = structural_body_forces(model, &resolved, &mesh, step)?; + for (s, body) in mesh.structural.iter().zip(&per_volume) { + let at = |err| Error::element(s.label, err); + let dofs = mesh + .layout + .element_dofs(&s.nodes, Structural::KINDS) + .ok_or_else(|| Error::Internal(format!("element {} degrees of freedom", s.label)))?; + let ue: Vec = dofs.iter().map(|d| u.get(*d).copied().unwrap_or(0.0)).collect(); + match s.kernel()? { + StructuralKernel::Shell(shell) => { + let stress = shell.stresses(&ue).map_err(at)?; + for (n, v) in s.nodes.iter().zip(&stress.top) { + if let (Some(acc), Some(c)) = (sum.get_mut(*n), count.get_mut(*n)) { + for k in 0..6 { + acc[k] += v[k]; + } + *c += 1; + } + } + shells.push(ShellResult { label: s.label, kind: s.kind, stress }); + } + StructuralKernel::Beam(beam) => beams.push(BeamResult { label: s.label, kind: s.kind, forces: beam.forces(&ue, *body).map_err(at)? }), + } + } let stress: Vec<[f64; 6]> = sum.iter().zip(&count).map(|(s, c)| s.map(|v| if *c > 0 { v / f64::from(*c) } else { 0.0 })).collect(); let von_mises = stress.iter().map(von_mises).collect(); let scale = applied_total.iter().chain(&reaction_total).fold(0.0f64, |m, v| m.max(v.abs())); @@ -349,10 +419,14 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul step: step.name.clone(), nodes: mesh.labels.clone(), displacement, + rotation, reaction, + reaction_moment, stress, von_mises, elements, + shells, + beams, reference_points, applied_total, reaction_total, @@ -398,7 +472,7 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt return Err(Error::Invalid(format!("the shift must be a finite number, got {s}"))); } let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; - let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).collect(); + let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).chain(mesh.warnings.iter().cloned()).collect(); warnings.extend(inelastic_warnings(model)); let (graph, order) = pattern(&mesh); ignored_dashpots(&mesh, &mut warnings); @@ -406,14 +480,18 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt let red = reduction(&setup, false, &order, budget(limits, &mesh, &graph, 2), &mut warnings)?; let Assembled { k, m } = assemble(&mesh, &graph, Some(options.mass))?; let m = m.ok_or_else(|| Error::Internal("mass matrix".into()))?; - // Every element's mass entries sum to its mass, once per direction. - // (the mean over the directions for a point mass given per direction; rotary inertia is left - // out), added in storage order. - let translation = |d: usize| mesh.kind(d).and_then(DofKind::direction).is_some_and(|c| c < 3); + // Every element's mass entries between translations sum to its mass, once per direction + // (the mean over the directions for a point mass given per direction; rotary inertia, and + // the translation–rotation coupling of shells and beams, are left out), added in storage + // order. + let translation: Vec = (0..m.n).map(|d| mesh.kind(d).and_then(DofKind::direction).is_some_and(|c| c < 3)).collect(); let mut translational = 0.0; - for d in (0..m.n).filter(|d| translation(*d)) { - for v in m.row(d).1 { - translational += v; + for d in (0..m.n).filter(|d| translation.get(*d).copied().unwrap_or(false)) { + let (cols, vals) = m.row(d); + for (c, v) in cols.iter().zip(vals) { + if translation.get(*c).copied().unwrap_or(false) { + translational += v; + } } } let total_mass = translational / mesh.translation_kinds().len() as f64; @@ -468,11 +546,12 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt let mut out = Vec::with_capacity(pairs.values.len()); for (lambda, phi) in pairs.values.iter().zip(&pairs.vectors) { let participation = [0, 1, 2].map(|d| mr.get(d).map_or(0.0, |y| y.iter().zip(phi).map(|(a, b)| a * b).sum())); - let (shape, point_shape) = split(&red.expand(phi, false), &mesh); + let (shape, rotation, point_shape) = split(&red.expand(phi, false), &mesh); out.push(Mode { eigenvalue: *lambda, frequency: lambda.max(0.0).sqrt() / two_pi, shape, + rotation, participation, effective_mass: participation.map(|g| g * g), point_shape, @@ -501,6 +580,8 @@ pub enum Location { /// The result of one step, as a session keeps it. #[derive(Debug, Clone, PartialEq)] +// One per solved step; boxing the larger variant would buy nothing. +#[allow(clippy::large_enum_variant)] pub enum StepResult { Static(StaticResult), Frequency(FrequencyResult), diff --git a/crates/solver/src/constraints.rs b/crates/solver/src/constraints.rs index 202a295..0937b13 100644 --- a/crates/solver/src/constraints.rs +++ b/crates/solver/src/constraints.rs @@ -10,7 +10,8 @@ //! D. S. Malkus, M. E. Plesha, R. J. Witt, *Concepts and Applications of Finite Element //! Analysis*, 4th ed., Wiley 2002, §9.2). The reference point's rotations take part through the //! lever arms, so a moment on it reaches the face as a couple even when only translations are -//! coupled; coupled rotations have nothing to act on at solid nodes, which have none. +//! coupled. Coupled rotations tie the rotations of shell and beam nodes, θ_s = θ_r; solid +//! nodes have none. //! - **Distributing** (an interpolation constraint): the reference point moves with the rigid //! motion that fits the set's nodal translations best in the least-squares sense, with equal //! weights (as the model and the deck writers define it). About the nodes' centroid c, with @@ -50,8 +51,8 @@ use std::collections::BTreeMap; -use modelcraft_fem::DofKind; use modelcraft_fem::dense::inverse3; +use modelcraft_fem::{DofKind, DofSet}; use modelcraft_model::{CouplingKind, Model, Resolved}; use crate::sparse::Csr; @@ -113,7 +114,8 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> let coupled = |range: std::ops::Range| range.filter(|k| c.dofs.contains(*k as u8 + 1)); match c.kind { CouplingKind::Kinematic => { - if coupled(0..3).next().is_none() { + let rotations = |n: usize| mesh.layout.kinds(n).includes(DofSet::ROTATIONS); + if coupled(0..3).next().is_none() && !(nodes.iter().any(|n| rotations(*n)) && coupled(3..6).next().is_some()) { return Err(Error::Invalid(format!( "coupling `{name}` couples only rotations, which the nodes of solid elements do not have; couple translations too" ))); @@ -133,6 +135,14 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> } rows.push(Row { terms, pivot, coupling: name }); } + // Shell and beam nodes also turn with the point: θ_s = θ_r. + if rotations(n) { + for k in coupled(3..6) { + let pivot = + DofKind::from_direction(k).ok_or_else(|| Error::Internal("direction".into())).and_then(|kind| mesh.dof(n, kind))?; + rows.push(Row { terms: vec![(pivot, 1.0), (rp(k)?, -1.0)], pivot, coupling: name }); + } + } } } CouplingKind::Distributing => { diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 943eaab..48e3962 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -7,15 +7,19 @@ //! translations for solids, Ux and Uy for 2D elements), nodes by label, then six per coupled reference point (points by //! name): translations along X, Y, Z and rotations about them. +use modelcraft_fem::beam::{Beam, BeamSection}; +use modelcraft_fem::dense::Matrix; use modelcraft_fem::element::{Formulation, Solid}; use modelcraft_fem::plane::{Plane, Theory}; +use modelcraft_fem::shell::{Shell, ShellSection}; use modelcraft_fem::truss::{Truss, young_modulus}; use modelcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, ElementDofs, PointStress, elasticity, face, plane}; use modelcraft_model::nodal::edge_pressure; use std::collections::BTreeSet; -use modelcraft_model::{ElementKind, Load, Model, Resolved, SectionKind, Step, Target, centrifugal_acceleration}; +use modelcraft_model::{ElementKind, Integration, Load, Model, Resolved, SectionKind, Step, Target, centrifugal_acceleration}; +use crate::analysis::MassKind; use crate::{Error, Result}; /// An element ready for the kernels: a solid, or a 2D element (`theory`). @@ -138,6 +142,79 @@ impl Bar { } } +/// What a shell or beam element is made of. +#[derive(Debug, Clone, Copy)] +pub enum Section { + Shell(ShellSection), + Beam { section: BeamSection, orientation: [f64; 3] }, +} + +/// A shell or beam element ready for its kernel: six degrees of freedom per node. +#[derive(Debug, Clone)] +pub struct Structural { + pub label: u32, + pub kind: ElementKind, + /// Indices into [`Mesh::labels`] (layout slots). + pub nodes: Vec, + pub xyz: Vec<[f64; 3]>, + pub density: Option, + pub section: Section, +} + +/// The kernel of a shell or beam element. +#[derive(Debug, Clone)] +pub enum StructuralKernel { + Shell(Shell), + Beam(Beam), +} + +impl Structural { + /// The unknowns every node of a shell or beam carries. + pub const KINDS: DofSet = DofSet::STRUCTURAL; + + pub fn kernel(&self) -> Result { + let k = match self.section { + Section::Shell(s) => Shell::new(self.kind, &self.xyz, s).map(StructuralKernel::Shell), + Section::Beam { section, orientation } => Beam::new(&self.xyz, orientation, section).map(StructuralKernel::Beam), + }; + k.map_err(|e| Error::element(self.label, e)) + } +} + +impl StructuralKernel { + pub fn stiffness(&self) -> modelcraft_fem::Result { + match self { + StructuralKernel::Shell(s) => s.stiffness(), + StructuralKernel::Beam(b) => Ok(b.stiffness()), + } + } + + /// The consistent or lumped mass over every element degree of freedom. + pub fn mass(&self, density: f64, kind: MassKind) -> modelcraft_fem::Result { + Ok(match (self, kind) { + (StructuralKernel::Shell(s), MassKind::Consistent) => s.mass(density)?, + (StructuralKernel::Shell(s), MassKind::Lumped) => { + let diagonal = s.lumped_mass(density)?; + let mut m = Matrix::zeros(diagonal.len()); + for (i, v) in diagonal.iter().enumerate() { + m.add(i, i, *v); + } + m + } + (StructuralKernel::Beam(b), MassKind::Consistent) => b.mass(density), + (StructuralKernel::Beam(b), MassKind::Lumped) => b.lumped_mass(density), + }) + } + + /// Consistent nodal forces and moments of a uniform force per volume, six per node. + pub fn body_force(&self, b: [f64; 3]) -> modelcraft_fem::Result> { + match self { + StructuralKernel::Shell(s) => s.body_force(b), + StructuralKernel::Beam(beam) => Ok(beam.body_force(b)), + } + } +} + /// The translations of a model's nodes: X and Y in a 2D model, else X, Y and Z. Trusses and /// translational discrete elements declare these. pub fn translation_kinds(planar: bool) -> DofSet { @@ -177,6 +254,10 @@ pub struct Mesh { /// The 2D theory when this is a 2D model (a model is 2D as a whole). pub theory: Option, pub trusses: Vec, + /// Shells and beams, by label. + pub structural: Vec, + /// What the solver changed or ignored in the model's sections. + pub warnings: Vec, /// Point masses, springs and dashpots ([`crate::discrete`]). pub discrete: crate::discrete::Discrete, } @@ -190,6 +271,11 @@ impl Mesh { self.trusses.binary_search_by_key(&label, |e| e.label).ok().and_then(|i| self.trusses.get(i)) } + /// The index of a shell or beam element in [`Self::structural`]. + pub fn structural_index(&self, label: u32) -> Option { + self.structural.binary_search_by_key(&label, |e| e.label).ok() + } + /// The kinds trusses and translational discrete elements carry at their nodes. pub fn translation_kinds(&self) -> DofSet { translation_kinds(self.theory.is_some()) @@ -239,6 +325,8 @@ impl Mesh { pub fn new(model: &Model, resolved: &Resolved, max_dofs: usize) -> Result { let mut elements = Vec::new(); let mut trusses = Vec::new(); + let mut structural = Vec::new(); + let mut warnings = Vec::new(); // A model is 2D as a whole (the model check refuses a mix). let planar = model.sections().values().find_map(|s| Theory::from_section(&s.kind)); if matches!(planar, Some(Theory::Axisymmetric)) @@ -257,7 +345,11 @@ impl Mesh { SectionKind::Rod { area } => Some(area), _ => None, }; - if !matches!(section.kind, SectionKind::Solid {}) && theory.is_none() && rod.is_none() { + let shell_or_beam = matches!(section.kind, SectionKind::Shell { .. } | SectionKind::Beam { .. }); + if shell_or_beam && planar.is_some() { + return Err(Error::NotSupported(format!("{} sections in 2D models (section `{name}`)", section.kind.name()))); + } + if !matches!(section.kind, SectionKind::Solid {}) && theory.is_none() && rod.is_none() && !shell_or_beam { return Err(Error::NotSupported(format!( "{} sections (section `{name}`); the solver has solid, 2D (plane stress, plane strain, axisymmetric) and truss elements so far", section.kind.name() @@ -272,6 +364,31 @@ impl Mesh { let elastic = material.elastic.as_ref().ok_or_else(|| Error::Invalid(format!("material `{}` has no elastic constants", section.material)))?; let set = resolved.sets.get(§ion.set).ok_or_else(|| Error::Invalid(format!("set `{}` is not resolved", section.set)))?; + let bad = |e: modelcraft_fem::Error| match e { + modelcraft_fem::Error::NotSupported(m) => Error::NotSupported(format!("section `{name}`: {m}")), + e => Error::Invalid(format!("section `{name}`, material `{}`: {e}", section.material)), + }; + let shared = match §ion.kind { + SectionKind::Shell { thickness, offset, .. } => { + if section.formulation.integration == Integration::Reduced { + warnings.push(format!( + "section `{name}`: shells are MITC elements with full integration; reduced integration is solved the same way" + )); + } + Some(Section::Shell(ShellSection::new(*thickness, *offset, elastic).map_err(bad)?)) + } + SectionKind::Beam { profile, orientation } => { + Some(Section::Beam { section: BeamSection::new(profile, elastic).map_err(bad)?, orientation: *orientation }) + } + _ => None, + }; + if let Some(body) = shared { + for &label in &set.elements { + let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; + structural.push((label, el.kind, el.nodes.clone(), material.density, body)); + } + continue; + } if let Some(area) = rod { let young = young_modulus(elastic).map_err(|e| Error::Invalid(format!("material `{}` of section `{name}`: {e}", section.material)))?; @@ -297,12 +414,14 @@ impl Mesh { } elements.sort_by_key(|e| e.0); trusses.sort_by_key(|e| e.0); + structural.sort_by_key(|e| e.0); let translations = translation_kinds(planar.is_some()); // What discrete elements need: at mesh nodes (some only they use), and at reference // points, which then get their six degrees of freedom like coupled ones. let needs: Vec<(u32, DofSet)> = crate::discrete::needs(resolved, translations).collect(); let mut labels: Vec = elements.iter().flat_map(|e| e.2.iter().copied()).collect(); labels.extend(trusses.iter().flat_map(|e| e.2.iter().copied())); + labels.extend(structural.iter().flat_map(|e| e.2.iter().copied())); labels.extend(needs.iter().map(|n| n.0).filter(|l| model.node(*l).is_some())); labels.sort_unstable(); labels.dedup(); @@ -323,6 +442,9 @@ impl Mesh { for e in &trusses { layout.element(&e.2, modelcraft_fem::ElementDofs::new(translations)); } + for e in &structural { + layout.element(&e.2, ElementDofs::new(Structural::KINDS)); + } for (label, kinds) in needs.iter().filter(|n| model.node(n.0).is_some()) { layout.node(*label, *kinds); } @@ -346,6 +468,8 @@ impl Mesh { layout, theory: planar, trusses: Vec::with_capacity(trusses.len()), + structural: Vec::with_capacity(structural.len()), + warnings, discrete: crate::discrete::Discrete::default(), }; for (label, kind, nodes, formulation, d, density, theory, _) in elements { @@ -358,6 +482,11 @@ impl Mesh { let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); mesh.trusses.push(Bar { label, kind, nodes: idx, xyz: exyz, young, area, density }); } + for (label, kind, nodes, density, section) in structural { + let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); + let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); + mesh.structural.push(Structural { label, kind, nodes: idx, xyz: exyz, density, section }); + } mesh.discrete = crate::discrete::Discrete::new(&mesh, resolved)?; Ok(mesh) } @@ -365,7 +494,12 @@ impl Mesh { /// The layout slots (mesh nodes, reference points) every element and discrete element /// couples, for the matrix pattern. pub fn blocks(&self) -> impl Iterator> { - let elements = self.elements.iter().map(|e| e.nodes.as_slice()).chain(self.trusses.iter().map(|t| t.nodes.as_slice())); + let elements = self + .elements + .iter() + .map(|e| e.nodes.as_slice()) + .chain(self.trusses.iter().map(|t| t.nodes.as_slice())) + .chain(self.structural.iter().map(|s| s.nodes.as_slice())); elements.map(std::borrow::Cow::Borrowed).chain(self.discrete.couplings().map(|dofs| { let mut slots: Vec = dofs.iter().filter_map(|d| self.layout.locate(*d).map(|(slot, _)| slot)).collect(); slots.dedup(); @@ -491,13 +625,16 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ "load `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh and no spring or point mass uses" ))); } - // Moments act on the rotations, which only reference points have. + // Moments act on the rotations, which reference points and the nodes of shells + // and beams have. let rotation = matches!(load, Load::Moment { .. }); for (label, v) in resolved.nodal_loads.get(name).map(Vec::as_slice).unwrap_or_default() { let slot = mesh.layout.slot(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: node {label} belongs to no solid element")))?; let dofs = vector_dofs(mesh, slot, rotation).ok_or_else(|| { - Error::Invalid(format!("moment `{name}`: solid nodes have no rotations; apply it to a coupled reference point")) + Error::Invalid(format!( + "moment `{name}`: node {label} has no rotations (only shell and beam nodes do); apply it to a coupled reference point" + )) })?; for (&dof, x) in dofs.iter().zip(v) { if *x != 0.0 && !active.get(dof).copied().unwrap_or(false) { @@ -564,6 +701,13 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ let rho = t.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", t.label)))?; bar_body(&mut f, mesh, t, acceleration.map(|a| rho * a))?; } + for e in &mesh.structural { + if chosen.is_some_and(|c| c.binary_search(&e.label).is_err()) { + continue; + } + let rho = e.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", e.label)))?; + structural_body(&mut f, mesh, e, acceleration.map(|a| rho * a))?; + } // Point masses take gravity on the whole model. if set.is_none() { mesh.discrete.add_gravity(&mut f, *acceleration, mesh); @@ -575,6 +719,10 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ bar_body(&mut f, mesh, t, *vector)?; continue; } + if let Some(e) = mesh.structural_index(*label).and_then(|i| mesh.structural.get(i)) { + structural_body(&mut f, mesh, e, *vector)?; + continue; + } let e = mesh.element(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: element {label} is not a solid element")))?; body(&mut f, mesh, e, *vector)?; } @@ -584,6 +732,9 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ if let Some(t) = mesh.trusses.iter().find(|t| chosen.is_none_or(|c| c.binary_search(&t.label).is_ok())) { return Err(Error::NotSupported(format!("centrifugal loads on trusses (load `{name}`, element {})", t.label))); } + if let Some(e) = mesh.structural.iter().find(|e| chosen.is_none_or(|c| c.binary_search(&e.label).is_ok())) { + return Err(Error::NotSupported(format!("centrifugal loads on shells and beams (load `{name}`, element {})", e.label))); + } if set.is_none() && !mesh.discrete.masses.is_empty() { return Err(Error::NotSupported(format!("centrifugal loads on point masses (load `{name}`)"))); } @@ -615,6 +766,49 @@ fn bar_body(f: &mut [f64], mesh: &Mesh, t: &Bar, b: [f64; 3]) -> Result<()> { Ok(()) } +/// Consistent forces and moments of a uniform force per volume on a shell or beam. +fn structural_body(f: &mut [f64], mesh: &Mesh, e: &Structural, b: [f64; 3]) -> Result<()> { + let forces = e.kernel()?.body_force(b).map_err(|err| Error::element(e.label, err))?; + let dofs = + mesh.layout.element_dofs(&e.nodes, Structural::KINDS).ok_or_else(|| Error::Internal(format!("element {} degrees of freedom", e.label)))?; + for (d, v) in dofs.iter().zip(&forces) { + if let Some(x) = f.get_mut(*d) { + *x += v; + } + } + Ok(()) +} + +/// The uniform force per volume a step's gravity and body forces put on each shell and beam +/// element ([`Mesh::structural`] order): beam end forces are recovered with it. +pub fn structural_body_forces(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step) -> Result> { + let mut out = vec![[0.0; 3]; mesh.structural.len()]; + for name in &step.loads { + let load = model.load(name).ok_or_else(|| Error::Invalid(format!("load `{name}` is missing")))?; + match load { + Load::Gravity { acceleration, set } => { + let chosen: Option<&[u32]> = set.as_ref().map(|s| resolved.sets.get(s).map(|r| r.elements.as_slice()).unwrap_or_default()); + for (e, b) in mesh.structural.iter().zip(out.iter_mut()) { + if chosen.is_some_and(|c| c.binary_search(&e.label).is_err()) { + continue; + } + let rho = e.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", e.label)))?; + *b = [0, 1, 2].map(|k| b[k] + rho * acceleration[k]); + } + } + Load::BodyForce { set, vector } => { + for label in resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default() { + if let Some(b) = mesh.structural_index(*label).and_then(|i| out.get_mut(i)) { + *b = [0, 1, 2].map(|k| b[k] + vector[k]); + } + } + } + _ => {} + } + } + Ok(out) +} + fn body(f: &mut [f64], mesh: &Mesh, e: &Elem, b: [f64; 3]) -> Result<()> { let forces = e.kernel()?.body_force(b).map_err(|err| Error::element(e.label, err))?; for (node, v) in e.nodes.iter().zip(forces) { diff --git a/crates/solver/tests/benchmarks.rs b/crates/solver/tests/benchmarks.rs new file mode 100644 index 0000000..1c5e046 --- /dev/null +++ b/crates/solver/tests/benchmarks.rs @@ -0,0 +1,182 @@ +//! Benchmark suite: shell problems compared with reference values published in papers, each +//! tagged with its source. These are not closed forms, so they are kept apart from the analytic +//! verification suite (`structural.rs`, `verification.rs`), which is the gate; a benchmark checks +//! convergence towards the published value within a stated margin. +//! +//! - **Scordelis–Lo roof** [benchmark: MacNeal & Harder 1985]: a cylindrical roof, R = 25, +//! L = 50, t = 0.25, 80° arc, E = 4.32 · 10⁸, ν = 0, self-weight 90 per area, rigid +//! diaphragms at the curved ends, straight edges free. Vertical deflection at the middle of a +//! free edge 0.3024 (R. H. MacNeal, R. L. Harder, *A proposed standard set of problems to test +//! finite element accuracy*, Finite Elements in Analysis and Design 1(1):3–20, 1985). +//! - **Pinched cylinder** [benchmark: Belytschko et al. 1985]: R = 300, L = 600, t = 3, +//! E = 3 · 10⁶, ν = 0.3, rigid diaphragms at both ends, two opposite unit forces at the middle. +//! Deflection under a load 1.8248 · 10⁻⁵ (T. Belytschko, H. Stolarski, W. K. Liu, N. Carpenter, +//! J. S.-J. Ong, *Stress projection for membrane and shear locking in shell finite elements*, +//! CMAME 51(1–3):221–258, 1985). +//! +//! Both are modelled on one symmetric part (a quarter of the roof, an eighth of the cylinder) +//! with quad4 (MITC4) shells on a structured mesh whose quadrilaterals are exactly flat. + +use modelcraft_model::{ + BoundaryCondition, Distribution, Elastic, ElementKind, Formulation, Load, Material, Model, Procedure, Section, SectionKind, Set, Step, Symmetry, + Target, +}; +use modelcraft_solver::{Limits, solve_static}; + +/// A cylindrical shell patch: x from 0 to `length`, angle φ from 0 (top, z = R) to `arc`, +/// y = R sin φ, z = R cos φ; n × n quad4 cells. Returns the model and the node at (x, φ). +struct Cylinder { + model: Model, + n: usize, +} + +impl Cylinder { + fn new(n: usize, radius: f64, length: f64, arc: f64, thickness: f64, elastic: Elastic) -> Cylinder { + let mut model = Model::default(); + let label = |i: usize, j: usize| (1 + i + (n + 1) * j) as u32; + for j in 0..=n { + for i in 0..=n { + let (x, phi) = (length * i as f64 / n as f64, arc * j as f64 / n as f64); + model.set_node(label(i, j), [x, radius * phi.sin(), radius * phi.cos()]).unwrap(); + } + } + let mut elements = Vec::new(); + for j in 0..n { + for i in 0..n { + let e = (1 + i + n * j) as u32; + // Counter-clockwise seen from outside: the normal points away from the axis. + model.set_element(e, ElementKind::Quad4, vec![label(i, j), label(i + 1, j), label(i + 1, j + 1), label(i, j + 1)]).unwrap(); + elements.push(e); + } + } + model.add_set("shell", Set { elements, ..Set::default() }).unwrap(); + model.add_material("Mat", Material { elastic: Some(elastic), density: None, ..Material::default() }).unwrap(); + model + .add_section( + "shell", + Section { + set: "shell".into(), + material: "Mat".into(), + kind: SectionKind::Shell { thickness, offset: 0.0, points: 5 }, + formulation: Formulation::default(), + }, + ) + .unwrap(); + model.add_step(Step::new("Step", Procedure::Static {}), None).unwrap(); + Cylinder { model, n } + } + + fn label(&self, i: usize, j: usize) -> u32 { + (1 + i + (self.n + 1) * j) as u32 + } + + fn nodes(&self, f: impl Fn(usize, usize) -> bool) -> Vec { + let mut out = Vec::new(); + for j in 0..=self.n { + for i in 0..=self.n { + if f(i, j) { + out.push(self.label(i, j)); + } + } + } + out.sort_unstable(); + out + } + + fn bc(&mut self, name: &str, nodes: Vec, bc: impl Fn(Target) -> BoundaryCondition) { + self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap(); + self.model.add_boundary_condition(name, bc(Target::Set(name.into()))).unwrap(); + self.model.set_boundary_condition_active("Step", name, true).unwrap(); + } + + /// A rigid diaphragm at x = length: no motion in its plane (y, z) and no rotation about x. + fn diaphragm(&mut self) { + let n = self.n; + let values = [None, Some(0.0), Some(0.0), Some(0.0), None, None]; + self.bc("diaphragm", self.nodes(|i, _| i == n), |target| BoundaryCondition::Displacement { target, values }); + } + + fn symmetry(&mut self, name: &str, plane: Symmetry, f: impl Fn(usize, usize) -> bool) { + let nodes = self.nodes(f); + self.bc(name, nodes, |target| BoundaryCondition::Symmetry { target, plane }); + } + + fn load(&mut self, name: &str, load: Load) { + self.model.add_load(name, load).unwrap(); + self.model.set_load_active("Step", name, true).unwrap(); + } +} + +/// The Scordelis–Lo roof's free-edge mid-span deflection on an n × n quarter mesh. +fn scordelis_lo(n: usize) -> f64 { + // A quarter: x from the mid-span (x = 0) to the diaphragm (x = 25), φ from the crown to the + // free edge at 40°. + let mut c = Cylinder::new(n, 25.0, 25.0, 40f64.to_radians(), 0.25, Elastic::Isotropic { young: 4.32e8, poisson: 0.0 }); + c.diaphragm(); + c.symmetry("midspan", Symmetry::X, |i, _| i == 0); + c.symmetry("crown", Symmetry::Y, |_, j| j == 0); + // Self-weight 90 per area: 360 per volume over t = 0.25. + c.load("weight", Load::BodyForce { set: "shell".into(), vector: [0.0, 0.0, -360.0] }); + let r = solve_static(&c.model, None, &Limits::default()).unwrap(); + assert!((r.reaction_total[2] + r.applied_total[2]).abs() < 1e-9 * r.applied_total[2].abs()); + let k = r.nodes.binary_search(&c.label(0, n)).unwrap(); + -r.displacement[k][2] +} + +/// Scordelis–Lo reference deflection (MacNeal & Harder 1985). +const SCORDELIS_LO: f64 = 0.3024; + +/// [benchmark: MacNeal & Harder 1985] Scordelis–Lo roof on a 16 × 16 quarter mesh: within 2 %. +#[test] +fn benchmark_scordelis_lo_roof_coarse() { + let ratio = scordelis_lo(16) / SCORDELIS_LO; + assert!((ratio - 1.0).abs() < 0.02, "{ratio}"); +} + +/// [benchmark: MacNeal & Harder 1985] Scordelis–Lo roof, free-edge mid-span deflection 0.3024: +/// monotone convergence to within 1 % at 64 × 64. +#[test] +#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"] +fn benchmark_scordelis_lo_roof() { + let results: Vec<(usize, f64, f64)> = [8, 16, 32, 64] + .iter() + .map(|&n| { + let w = scordelis_lo(n); + (n, w, w / SCORDELIS_LO) + }) + .collect(); + eprintln!("Scordelis–Lo roof (reference {SCORDELIS_LO}): (cells per side, deflection, ratio) {results:.4?}"); + let (_, _, finest) = results[results.len() - 1]; + assert!((finest - 1.0).abs() < 0.01, "{results:?}"); + assert!(results.windows(2).all(|w| w[1].2 > w[0].2), "converges from below: {results:?}"); +} + +/// [benchmark: Belytschko et al. 1985] Pinched cylinder with diaphragms, deflection under the +/// load 1.8248 · 10⁻⁵. +#[test] +#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"] +fn benchmark_pinched_cylinder() { + const REFERENCE: f64 = 1.8248e-5; + let mut results = Vec::new(); + for n in [8, 16, 32] { + // An eighth: x from the loaded section (x = 0) to the diaphragm (x = 300), φ from the + // loaded top to the side (z = 0). + let mut c = Cylinder::new(n, 300.0, 300.0, std::f64::consts::FRAC_PI_2, 3.0, Elastic::Isotropic { young: 3e6, poisson: 0.3 }); + c.diaphragm(); + c.symmetry("middle", Symmetry::X, |i, _| i == 0); + c.symmetry("top", Symmetry::Y, |_, j| j == 0); + c.symmetry("side", Symmetry::Z, |_, j| j == n); + let top = c.label(0, 0); + c.model.add_set("pinch", Set { nodes: vec![top], ..Set::default() }).unwrap(); + // A quarter of the unit force acts on the eighth. + c.load("pinch", Load::Force { target: Target::Set("pinch".into()), vector: [0.0, 0.0, -0.25], distribution: Distribution::Total }); + let r = solve_static(&c.model, None, &Limits::default()).unwrap(); + let k = r.nodes.binary_search(&top).unwrap(); + let w = -r.displacement[k][2]; + results.push((n, w, w / REFERENCE)); + } + eprintln!("pinched cylinder (reference {REFERENCE:e}): (cells per side, deflection, ratio) {results:.4?}"); + let (_, _, finest) = results[results.len() - 1]; + assert!((finest - 1.0).abs() < 0.02, "{results:?}"); + assert!(results.windows(2).all(|w| w[1].2 > w[0].2), "converges from below: {results:?}"); +} diff --git a/crates/solver/tests/structural.rs b/crates/solver/tests/structural.rs new file mode 100644 index 0000000..82c9eef --- /dev/null +++ b/crates/solver/tests/structural.rs @@ -0,0 +1,591 @@ +//! Analytic verification of the structural elements: shell patch tests (membrane and bending, +//! quad4 and tri3, in a tilted plane), the simply supported plate under uniform pressure against +//! Navier's series (deflection and first frequency), cantilever beams against Timoshenko and +//! Euler–Bernoulli theory (deflection and frequencies), a portal frame against the +//! slope-deflection hand calculation, and mixed meshes of solids with shells and rods. +//! +//! Every expected value is a closed form or a converged series; nothing comes from another +//! program. + +use std::collections::BTreeMap; + +use modelcraft_model::{ + BeamProfile, BoundaryCondition, Distribution, Elastic, ElementFace, ElementKind, Formulation, Load, Material, Model, Procedure, Section, + SectionKind, Set, Step, Target, +}; +use modelcraft_solver::{FrequencyOptions, Limits, MassKind, StaticResult, solve_frequency, solve_static}; + +/// A model under construction: nodes and elements numbered from 1, one material, one step. +struct Build { + model: Model, + xyz: BTreeMap, + elements: u32, +} + +impl Build { + fn new(elastic: Elastic, density: f64, procedure: Procedure) -> Build { + let mut model = Model::default(); + model.add_material("Mat", Material { elastic: Some(elastic), density: Some(density), ..Material::default() }).unwrap(); + model.add_step(Step::new("Step", procedure), None).unwrap(); + Build { model, xyz: BTreeMap::new(), elements: 0 } + } + + fn node(&mut self, label: u32, xyz: [f64; 3]) { + self.model.set_node(label, xyz).unwrap(); + self.xyz.insert(label, xyz); + } + + fn element(&mut self, kind: ElementKind, nodes: Vec) -> u32 { + self.elements += 1; + self.model.set_element(self.elements, kind, nodes).unwrap(); + self.elements + } + + fn section(&mut self, name: &str, elements: Vec, kind: SectionKind) { + self.model.add_set(name, Set { elements, ..Set::default() }).unwrap(); + self.model.add_section(name, Section { set: name.into(), material: "Mat".into(), kind, formulation: Formulation::default() }).unwrap(); + } + + fn nodes_where(&self, f: impl Fn([f64; 3]) -> bool) -> Vec { + self.xyz.iter().filter(|(_, x)| f(**x)).map(|(l, _)| *l).collect() + } + + /// A boundary condition on the nodes where `f` holds: `values` per direction (None = free). + fn hold(&mut self, name: &str, values: [Option; 6], f: impl Fn([f64; 3]) -> bool) { + let nodes = self.nodes_where(f); + self.hold_nodes(name, values, nodes); + } + + fn hold_nodes(&mut self, name: &str, values: [Option; 6], nodes: Vec) { + self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap(); + self.model.add_boundary_condition(name, BoundaryCondition::Displacement { target: Target::Set(name.into()), values }).unwrap(); + self.model.set_boundary_condition_active("Step", name, true).unwrap(); + } + + fn load(&mut self, name: &str, load: Load) { + self.model.add_load(name, load).unwrap(); + self.model.set_load_active("Step", name, true).unwrap(); + } + + fn force(&mut self, name: &str, nodes: Vec, vector: [f64; 3]) { + self.model.add_set(name, Set { nodes, ..Set::default() }).unwrap(); + self.load(name, Load::Force { target: Target::Set(name.into()), vector, distribution: Distribution::Total }); + } + + /// A structured surface mesh of n[0] × n[1] cells on the map of [0, 1]², quadrilaterals or + /// two triangles per cell, starting at node label `base`; returns the elements. + fn surface(&mut self, n: [usize; 2], kind: ElementKind, base: u32, map: impl Fn(f64, f64) -> [f64; 3]) -> Vec { + let label = |i: usize, j: usize| base + (i + (n[0] + 1) * j) as u32; + for j in 0..=n[1] { + for i in 0..=n[0] { + self.node(label(i, j), map(i as f64 / n[0] as f64, j as f64 / n[1] as f64)); + } + } + let mut out = Vec::new(); + for j in 0..n[1] { + for i in 0..n[0] { + let (a, b, c, d) = (label(i, j), label(i + 1, j), label(i + 1, j + 1), label(i, j + 1)); + if kind == ElementKind::Quad4 { + out.push(self.element(kind, vec![a, b, c, d])); + } else { + out.push(self.element(kind, vec![a, b, c])); + out.push(self.element(kind, vec![a, c, d])); + } + } + } + out + } + + fn solve(&self) -> StaticResult { + solve_static(&self.model, None, &Limits::default()).unwrap() + } +} + +fn at(r: &StaticResult, label: u32) -> usize { + r.nodes.binary_search(&label).unwrap() +} + +fn rel(a: f64, b: f64) -> f64 { + (a - b).abs() / b.abs() +} + +const FREE: Option = None; + +// --------------------------------------------------------------------------------------------- +// Shell patch tests +// --------------------------------------------------------------------------------------------- + +/// The distorted patch of five quadrilaterals in a 0.24 × 0.12 rectangle (the irregular patch of +/// R. H. MacNeal, R. L. Harder, *A proposed standard set of problems to test finite element +/// accuracy*, Finite Elements in Analysis and Design 1(1):3–20, 1985), in local plane +/// coordinates; the four inner nodes are 5–8. +const PATCH: [[f64; 2]; 8] = [[0.0, 0.0], [0.24, 0.0], [0.24, 0.12], [0.0, 0.12], [0.04, 0.02], [0.18, 0.03], [0.16, 0.08], [0.08, 0.08]]; +/// Thick enough that membrane and bending stiffness are within a few orders of magnitude, so the +/// comparison measures the elements and not the round-off of an ill-conditioned system. +const PATCH_THICKNESS: f64 = 0.02; +const PATCH_QUADS: [[u32; 4]; 5] = [[1, 2, 6, 5], [2, 3, 7, 6], [3, 4, 8, 7], [4, 1, 5, 8], [5, 6, 7, 8]]; + +/// A tilted plane: local (x, y) → global, with its axes. +fn tilted(p: [f64; 2]) -> [f64; 3] { + let (e1, e2) = axes(); + [0, 1, 2].map(|k| 1.0 + p[0] * e1[k] + p[1] * e2[k]) +} + +fn axes() -> ([f64; 3], [f64; 3]) { + let e1 = [0.6, 0.8, 0.0]; + let (c, s) = (0.5f64.sqrt(), 0.5f64.sqrt()); + // e2 = cos · (−0.8, 0.6, 0) + sin · z. + let e2 = [-0.8 * c, 0.6 * c, s]; + (e1, e2) +} + +fn normal() -> [f64; 3] { + let (a, b) = axes(); + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] +} + +/// Local displacement (u, v, w) and rotations (θx, θy, θz) to global 6-vector. +fn to_global(local: [f64; 6]) -> [f64; 6] { + let (e1, e2) = axes(); + let e3 = normal(); + let g = |a: [f64; 3]| [0, 1, 2].map(|k| a[0] * e1[k] + a[1] * e2[k] + a[2] * e3[k]); + let t = g([local[0], local[1], local[2]]); + let r = g([local[3], local[4], local[5]]); + [t[0], t[1], t[2], r[0], r[1], r[2]] +} + +/// The patch with a prescribed exact field on its boundary nodes 1–4: returns the worst nodal +/// error of the inner nodes (relative to the field's size) and the solution. +fn patch(kind: ElementKind, field: impl Fn([f64; 2]) -> [f64; 6]) -> (f64, StaticResult, Build) { + let mut b = Build::new(Elastic::Isotropic { young: 1e6, poisson: 0.25 }, 1.0, Procedure::Static {}); + for (i, p) in PATCH.iter().enumerate() { + b.node(i as u32 + 1, tilted(*p)); + } + let mut elements = Vec::new(); + for q in PATCH_QUADS { + if kind == ElementKind::Quad4 { + elements.push(b.element(kind, q.to_vec())); + } else { + elements.push(b.element(kind, vec![q[0], q[1], q[2]])); + elements.push(b.element(kind, vec![q[0], q[2], q[3]])); + } + } + b.section("shell", elements, SectionKind::Shell { thickness: PATCH_THICKNESS, offset: 0.0, points: 5 }); + for (i, p) in PATCH.iter().enumerate().take(4) { + let g = to_global(field(*p)); + b.hold_nodes(&format!("edge{i}"), g.map(Some), vec![i as u32 + 1]); + } + let r = b.solve(); + let size = PATCH.iter().flat_map(|p| field(*p)).fold(0.0f64, |m, v| m.max(v.abs())); + let mut worst = 0.0f64; + for (i, p) in PATCH.iter().enumerate().skip(4) { + let want = to_global(field(*p)); + let k = at(&r, i as u32 + 1); + let got = [r.displacement[k], r.rotation[k]].concat(); + for c in 0..6 { + worst = worst.max((got[c] - want[c]).abs() / size); + } + } + (worst, r, b) +} + +#[test] +fn shell_membrane_patch_test() { + // u = 10⁻³ (x + y/2), v = 10⁻³ (y + x/2): constant strain, no in-plane rotation. + let field = |p: [f64; 2]| [1e-3 * (p[0] + p[1] / 2.0), 1e-3 * (p[1] + p[0] / 2.0), 0.0, 0.0, 0.0, 0.0]; + for kind in [ElementKind::Quad4, ElementKind::Tri3] { + let (err, r, _) = patch(kind, field); + // σ in local axes: E/(1 − ν²) [[1, ν], [ν, 1]] ε, τ = G γ. + let (e, nu) = (1e6, 0.25); + let c = e / (1.0 - nu * nu); + let (sx, sy, txy) = (c * (1e-3 + nu * 1e-3), c * (1e-3 + nu * 1e-3), e / (2.0 * (1.0 + nu)) * 1e-3); + let mut stress_err = 0.0f64; + for s in &r.shells { + let [a1, a2, _] = s.stress.axes; + for layer in [&s.stress.bottom, &s.stress.middle, &s.stress.top] { + for g in layer { + // Back to the patch's local axes: σ_ab = a·S·b. + let t = [[g[0], g[3], g[4]], [g[3], g[1], g[5]], [g[4], g[5], g[2]]]; + let (e1, e2) = axes(); + let p = |a: [f64; 3], b: [f64; 3]| (0..3).map(|i| (0..3).map(|j| a[i] * t[i][j] * b[j]).sum::()).sum::(); + stress_err = stress_err.max(((p(e1, e1) - sx).abs()).max((p(e2, e2) - sy).abs()).max((p(e1, e2) - txy).abs()) / sx); + } + } + assert!(a1.iter().chain(&a2).all(|v| v.is_finite())); + } + eprintln!("membrane patch {kind:?}: displacement error {err:.1e}, stress error {stress_err:.1e}"); + assert!(err < 1e-9 && stress_err < 1e-9, "{kind:?}: {err:e} {stress_err:e}"); + } +} + +#[test] +fn shell_bending_patch_test() { + // w = 10⁻³ (x² + x y + y²)/2 with θx = w,y and θy = −w,x: constant curvature, no shear. + let k = 1e-3; + let field = |p: [f64; 2]| { + let [x, y] = p; + [0.0, 0.0, k * (x * x + x * y + y * y) / 2.0, k * (x / 2.0 + y), -k * (x + y / 2.0), 0.0] + }; + for kind in [ElementKind::Quad4, ElementKind::Tri3] { + let (err, r, _) = patch(kind, field); + // κx = θy,x = −k, κy = −θx,y = −k, κxy = θy,y − θx,x = −k. + let (e, nu, t) = (1e6, 0.25, PATCH_THICKNESS); + let d = e * t * t * t / (12.0 * (1.0 - nu * nu)); + let want = [d * (-k - nu * k), d * (-k - nu * k), d * (1.0 - nu) / 2.0 * -k]; + let mut moment_err = 0.0f64; + for s in &r.shells { + // Local element axes may differ from the patch's in-plane axes by a rotation in the + // plane; the invariants Mx + My and Mx My − Mxy² do not. + for f in &s.stress.forces { + let (tr, det) = (f[3] + f[4], f[3] * f[4] - f[5] * f[5]); + let (wtr, wdet) = (want[0] + want[1], want[0] * want[1] - want[2] * want[2]); + moment_err = moment_err.max(rel(tr, wtr)).max(rel(det, wdet)).max(f[6].abs().max(f[7].abs()) / wtr.abs()); + } + } + eprintln!("bending patch {kind:?}: displacement error {err:.1e}, moment error {moment_err:.1e}"); + assert!(err < 1e-9 && moment_err < 1e-8, "{kind:?}: {err:e} {moment_err:e}"); + } +} + +// --------------------------------------------------------------------------------------------- +// Simply supported plate (Navier) +// --------------------------------------------------------------------------------------------- + +/// Navier's series for the centre deflection of a simply supported square plate under uniform +/// pressure: w = 16 q a⁴/(π⁶ D) Σ_{m, n odd} (−1)^((m+n)/2 − 1) / (m n (m² + n²)²) +/// (S. Timoshenko, S. Woinowsky-Krieger, *Theory of Plates and Shells*, 2nd ed., McGraw-Hill +/// 1959, §30). +fn navier(q: f64, a: f64, d: f64) -> f64 { + let mut s = 0.0; + for m in (1..400).step_by(2) { + for n in (1..400).step_by(2) { + let (mf, nf) = (f64::from(m), f64::from(n)); + let sign = if ((m + n) / 2 - 1) % 2 == 0 { 1.0 } else { -1.0 }; + s += sign / (mf * nf * (mf * mf + nf * nf).powi(2)); + } + } + 16.0 * q * a.powi(4) / (std::f64::consts::PI.powi(6) * d) * s +} + +/// A quarter of a simply supported square plate (side a, edges x = 0 and y = 0 held, symmetry +/// on x = y = a/2), n × n cells, under pressure q or for its modes. +fn plate(kind: ElementKind, n: usize, procedure: Procedure) -> Build { + let (a, t, e, nu, rho) = (1.0, 0.01, 1e7, 0.3, 1.0); + let mut b = Build::new(Elastic::Isotropic { young: e, poisson: nu }, rho, procedure); + let elements = b.surface([n, n], kind, 1, |u, v| [u * a / 2.0, v * a / 2.0, 0.0]); + b.section("plate", elements.clone(), SectionKind::Shell { thickness: t, offset: 0.0, points: 5 }); + let on = |x: f64, c: f64| (x - c).abs() < 1e-9; + // Hard simple support: w and the twist along the edge held. + b.hold("x0", [FREE, FREE, Some(0.0), Some(0.0), FREE, FREE], |p| on(p[0], 0.0)); + b.hold("y0", [FREE, FREE, Some(0.0), FREE, Some(0.0), FREE], |p| on(p[1], 0.0)); + b.hold("xs", [Some(0.0), FREE, FREE, FREE, Some(0.0), Some(0.0)], |p| on(p[0], a / 2.0)); + b.hold("ys", [FREE, Some(0.0), FREE, Some(0.0), FREE, Some(0.0)], |p| on(p[1], a / 2.0)); + if matches!(b.model.steps()[0].procedure, Procedure::Static {}) { + let faces = elements.iter().map(|e| ElementFace { element: *e, face: 1 }).collect(); + b.model.add_set("top", Set { faces, ..Set::default() }).unwrap(); + b.load("q", Load::Pressure { set: "top".into(), magnitude: 1.0 }); + } + b +} + +#[test] +fn simply_supported_plate_deflection_and_frequency() { + let (a, t, e, nu, rho) = (1.0, 0.01, 1e7, 0.3, 1.0); + let d = e * t * t * t / (12.0 * (1.0 - nu * nu)); + let w = navier(1.0, a, d); + let f11 = std::f64::consts::PI * 2.0 / (a * a) * (d / (rho * t)).sqrt() / 2.0; + for kind in [ElementKind::Quad4, ElementKind::Tri3] { + let mut errors = Vec::new(); + for n in [4, 8, 16] { + let b = plate(kind, n, Procedure::Static {}); + let r = b.solve(); + let centre = b.nodes_where(|p| (p[0] - 0.5).abs() < 1e-9 && (p[1] - 0.5).abs() < 1e-9)[0]; + let got = -r.displacement[at(&r, centre)][2]; + assert!((r.reaction_total[2] + r.applied_total[2]).abs() < 1e-9 * r.applied_total[2].abs()); + let b = plate(kind, n, Procedure::Frequency { modes: 1, min_hz: None, max_hz: None }); + let fr = solve_frequency(&b.model, None, &FrequencyOptions::default(), &Limits::default()).unwrap(); + errors.push((n, rel(got, w), rel(fr.modes[0].frequency, f11))); + } + eprintln!("plate {kind:?} vs Navier w = {w:.6e}, f11 = {f11:.6}: (cells, deflection error, frequency error) {errors:.5?}"); + let (_, dw, df) = errors[2]; + assert!(dw < 0.01 && df < 0.01, "{kind:?} {errors:?}"); + assert!(errors[0].1 > errors[2].1, "converges: {errors:?}"); + } +} + +// --------------------------------------------------------------------------------------------- +// Beams +// --------------------------------------------------------------------------------------------- + +/// A straight cantilever along a skew direction, `n` beam elements, clamped at node 1. +fn cantilever(n: usize, length: f64, profile: BeamProfile, procedure: Procedure) -> (Build, [f64; 3], [f64; 3]) { + let dir = [2.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0]; + let orientation = [0.0, 0.0, 1.0]; + let mut b = Build::new(Elastic::Isotropic { young: 210e3, poisson: 0.3 }, 7.85e-9, procedure); + for i in 0..=n { + let s = length * i as f64 / n as f64; + b.node(i as u32 + 1, dir.map(|c| c * s)); + } + let elements: Vec = (1..=n as u32).map(|i| b.element(ElementKind::Beam2, vec![i, i + 1])).collect(); + b.section("beam", elements, SectionKind::Beam { profile, orientation }); + b.hold_nodes("clamp", [Some(0.0); 6], vec![1]); + // Local y: the orientation made normal to the axis; local z = x × y. + let dot = dir[2]; + let y = [0, 1, 2].map(|k| orientation[k] - dot * dir[k]); + let ny = (y[0] * y[0] + y[1] * y[1] + y[2] * y[2]).sqrt(); + let y = y.map(|c| c / ny); + let z = [dir[1] * y[2] - dir[2] * y[1], dir[2] * y[0] - dir[0] * y[2], dir[0] * y[1] - dir[1] * y[0]]; + (b, y, z) +} + +#[test] +fn cantilever_beam_deflection_is_exact() { + // A stocky I beam, so shear deformation matters: δ = P L³/(3 E I) + P L/(G A_s) in each plane. + let profile = BeamProfile::I { height: 200.0, width: 100.0, web: 6.0, flange: 10.0 }; + let length: f64 = 1000.0; + let p = modelcraft_fem::beam::properties(&profile, 0.3); + let (e, g) = (210e3, 210e3 / 2.6); + for n in [1, 4] { + for (along_y, name) in [(true, "weak"), (false, "strong")] { + let (mut b, y, z) = cantilever(n, length, profile, Procedure::Static {}); + let d = if along_y { y } else { z }; + b.force("tip", vec![n as u32 + 1], d.map(|c| 1000.0 * c)); + let r = b.solve(); + let u = r.displacement[at(&r, n as u32 + 1)]; + let got = u[0] * d[0] + u[1] * d[1] + u[2] * d[2]; + let (i, a_s) = if along_y { (p.iz, p.shear_y) } else { (p.iy, p.shear_z) }; + let want = 1000.0 * length.powi(3) / (3.0 * e * i) + 1000.0 * length / (g * a_s); + eprintln!( + "cantilever {n} element(s), {name} axis: tip {got:.6e} vs {want:.6e} (shear share {:.1}%)", + 100.0 * 1000.0 * length / (g * a_s) / want + ); + assert!(rel(got, want) < 1e-9, "{got} {want}"); + // The clamp's moment and the root section force. + let m = r.reaction_moment[at(&r, 1)]; + let mm = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]).sqrt(); + assert!(rel(mm, 1000.0 * length) < 1e-9, "{m:?}"); + let root = &r.beams[0].forces.ends[0]; + let (shear, moment) = if along_y { (root[1], root[5]) } else { (root[2], root[4]) }; + assert!(rel(shear.abs(), 1000.0) < 1e-9 && rel(moment.abs(), 1000.0 * length) < 1e-9, "{root:?}"); + } + } +} + +#[test] +fn cantilever_beam_frequencies() { + // A slender pipe, L/r = 1000 (shear and rotary inertia below 10⁻⁴ in the third mode): + // Euler–Bernoulli ω_k = (β_k L)² √(E I / (ρ A L⁴)). + let profile = BeamProfile::Pipe { radius: 10.0, thickness: 1.0 }; + let length: f64 = 10000.0; + let p = modelcraft_fem::beam::properties(&profile, 0.3); + let base = (210e3 * p.iy / (7.85e-9 * p.area * length.powi(4))).sqrt() / (2.0 * std::f64::consts::PI); + let beta = [1.875_104_068_711_961, 4.694_091_132_974_175, 7.854_757_438_237_613]; + for mass in [MassKind::Consistent, MassKind::Lumped] { + let (b, _, _) = cantilever(40, length, profile, Procedure::Frequency { modes: 6, min_hz: None, max_hz: None }); + let r = solve_frequency(&b.model, None, &FrequencyOptions { mass, shift: None }, &Limits::default()).unwrap(); + // The pipe bends alike in both planes: modes come in pairs. + let errors: Vec = beta.iter().enumerate().map(|(k, bl)| r.modes[2 * k].frequency / (bl * bl * base) - 1.0).collect(); + eprintln!("cantilever frequencies ({mass:?}, 40 elements) vs Euler–Bernoulli: {errors:.5?}"); + let mass_total = 7.85e-9 * p.area * length; + assert!(rel(r.total_mass, mass_total) < 1e-12, "{} {mass_total}", r.total_mass); + for (k, e) in errors.iter().enumerate() { + // Lumped mass converges from below, a few 10⁻³ low at 40 elements. + let tolerance = if mass == MassKind::Consistent { 5e-4 } else { 5e-3 }; + assert!(e.abs() < tolerance, "mode {k}: {e}"); + assert!(rel(r.modes[2 * k].frequency, r.modes[2 * k + 1].frequency) < 1e-6); + } + } +} + +#[test] +fn portal_frame_matches_the_slope_deflection_solution() { + // Columns of height h, a beam of span l, fixed bases, a horizontal load H at the top. + // Slope-deflection with Timoshenko members (Przemieniecki 1968, §5.6): with + // a = E I_c / (h (1 + Φ_c)) and b = E I_b / (l (1 + Φ_b)), Φ = 12 E I / (G A_s L²), the end + // moments are M_ij = k [(4 + Φ) θ_i + (2 − Φ) θ_j − 6 Δ/L]; antisymmetric sway gives the joint + // rotation θ = 6 a Δ / (h (a (4 + Φ_c) + 6 b)), the storey shear + // a (6 θ − 12 Δ/h) = −H h / 2, and the base moment M = a ((2 − Φ_c) θ − 6 Δ/h). The columns' + // axial deformation is left out of the hand calculation: their shortening and lengthening + // under the overturning moment tilts the girder by about 2 H h² / (l² E A), about 2.6 · 10⁻³ + // of θ here; sway and base moment feel it far less. + let (h, l, load) = (8000.0, 10000.0, 1000.0); + let column = BeamProfile::Rectangle { width: 100.0, height: 100.0 }; + let girder = BeamProfile::Rectangle { width: 100.0, height: 200.0 }; + let e = 210e3; + let mut b = Build::new(Elastic::Isotropic { young: e, poisson: 0.3 }, 7.85e-9, Procedure::Static {}); + let n = 8; + // Nodes: left column 1..=n+1 (base to top), beam to the right column top, right column down. + let mut path = Vec::new(); + for i in 0..=n { + path.push([0.0, 0.0, h * i as f64 / n as f64]); + } + for i in 1..=n { + path.push([l * i as f64 / n as f64, 0.0, h]); + } + for i in 1..=n { + path.push([l, 0.0, h - h * i as f64 / n as f64]); + } + for (i, x) in path.iter().enumerate() { + b.node(i as u32 + 1, *x); + } + let count = path.len() as u32; + let mut columns = Vec::new(); + let mut beam = Vec::new(); + for i in 1..count { + let el = b.element(ElementKind::Beam2, vec![i, i + 1]); + if (n as u32 + 1..2 * n as u32 + 1).contains(&i) { beam.push(el) } else { columns.push(el) } + } + // Both members bend in the x–z plane: their profile height (local z) must lie in it, so + // the orientation (local y, the width) is along global y. + b.section("columns", columns, SectionKind::Beam { profile: column, orientation: [0.0, 1.0, 0.0] }); + b.section("girder", beam, SectionKind::Beam { profile: girder, orientation: [0.0, 1.0, 0.0] }); + b.hold_nodes("bases", [Some(0.0); 6], vec![1, count]); + b.force("push", vec![n as u32 + 1], [load, 0.0, 0.0]); + let r = b.solve(); + let g = e / 2.6; + let (pc, pb) = (modelcraft_fem::beam::properties(&column, 0.3), modelcraft_fem::beam::properties(&girder, 0.3)); + let (phi_c, phi_b) = (12.0 * e * pc.iy / (g * pc.shear_z * h * h), 12.0 * e * pb.iy / (g * pb.shear_z * l * l)); + let (a, bb) = (e * pc.iy / (h * (1.0 + phi_c)), e * pb.iy / (l * (1.0 + phi_b))); + let c = 6.0 * a / (h * (a * (4.0 + phi_c) + 6.0 * bb)); + let sway = -load * h / 2.0 / (a * (6.0 * c - 12.0 / h)); + let theta = c * sway; + let base = a * ((2.0 - phi_c) * theta - 6.0 * sway / h); + let top = r.displacement[at(&r, n as u32 + 1)][0]; + let m = r.reaction_moment[at(&r, 1)][1]; + let rot = r.rotation[at(&r, n as u32 + 1)][1]; + eprintln!( + "portal frame: sway {top:.6e} vs {sway:.6e} ({:.1e}), joint rotation {rot:.6e} vs {theta:.6e} ({:.1e}), base moment {m:.6e} vs {:.6e} ({:.1e})", + rel(top, sway), + rel(rot.abs(), theta.abs()), + base.abs(), + rel(m.abs(), base.abs()) + ); + assert!(rel(top, sway) < 5e-4 && rel(rot.abs(), theta.abs()) < 3e-3 && rel(m.abs(), base.abs()) < 5e-4); + assert!((r.reaction_total[0] + load).abs() < 1e-9 * load); +} + +// --------------------------------------------------------------------------------------------- +// Mixed meshes +// --------------------------------------------------------------------------------------------- + +/// A hex8 column 1 × 1 × 2 (two elements) clamped at z = 0, with four rods of area 1/4 from its +/// top corners up to four loaded nodes: the elongation is the sum of both parts' (ν = 0). +#[test] +fn rods_continue_a_solid_column() { + let e = 1000.0; + let mut b = Build::new(Elastic::Isotropic { young: e, poisson: 0.0 }, 1.0, Procedure::Static {}); + let corners = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; + for z in 0..3 { + for (i, c) in corners.iter().enumerate() { + b.node(1 + i as u32 + 4 * z, [c[0], c[1], f64::from(z)]); + } + } + for (i, c) in corners.iter().enumerate() { + b.node(13 + i as u32, [c[0], c[1], 5.0]); + } + let hexes: Vec = (0..2).map(|k| b.element(ElementKind::Hex8, (1..=8).map(|n| n + 4 * k).collect())).collect(); + let rods: Vec = (0..4).map(|i| b.element(ElementKind::Rod2, vec![9 + i, 13 + i])).collect(); + b.section("column", hexes, SectionKind::Solid {}); + b.section("rods", rods, SectionKind::Rod { area: 0.25 }); + b.hold_nodes("base", [Some(0.0), Some(0.0), Some(0.0), FREE, FREE, FREE], vec![1, 2, 3, 4]); + // Rods carry no lateral stiffness: hold their ends sideways. + b.hold_nodes("guides", [Some(0.0), Some(0.0), FREE, FREE, FREE, FREE], vec![13, 14, 15, 16]); + b.force("pull", vec![13, 14, 15, 16], [0.0, 0.0, 4.0]); + let r = b.solve(); + let want = 4.0 * 2.0 / (e * 1.0) + 4.0 * 3.0 / (e * 1.0); + for n in 13..=16 { + let got = r.displacement[at(&r, n)][2]; + assert!(rel(got, want) < 1e-12, "{got} {want}"); + } + let rods: Vec<_> = r.elements.iter().filter(|e| e.kind == ElementKind::Rod2).collect(); + assert_eq!(rods.len(), 4); + assert!(rods.iter().flat_map(|e| &e.nodal).all(|s| rel(modelcraft_fem::von_mises(s), 4.0) < 1e-12), "{rods:?}"); + // Solid nodes have no rotations. + assert!(r.rotation.iter().all(|x| *x == [0.0; 3])); +} + +/// A hex8 slab with a shell skin on its top face sharing the nodes, stretched by prescribed end +/// displacements: the reaction is E (A_solid + A_shell) ε (ν = 0, uniform strain in both). +#[test] +fn shell_skin_on_a_solid_slab() { + let (e, length, width, depth, skin) = (1000.0, 4.0, 2.0, 1.0, 0.1); + let mut b = Build::new(Elastic::Isotropic { young: e, poisson: 0.0 }, 1.0, Procedure::Static {}); + let (nx, ny) = (4usize, 2usize); + let label = |i: usize, j: usize, k: usize| (1 + i + (nx + 1) * (j + (ny + 1) * k)) as u32; + for k in 0..2 { + for j in 0..=ny { + for i in 0..=nx { + b.node(label(i, j, k), [length * i as f64 / nx as f64, width * j as f64 / ny as f64, depth * k as f64]); + } + } + } + let mut hexes = Vec::new(); + let mut quads = Vec::new(); + for j in 0..ny { + for i in 0..nx { + let q = |k| vec![label(i, j, k), label(i + 1, j, k), label(i + 1, j + 1, k), label(i, j + 1, k)]; + hexes.push(b.element(ElementKind::Hex8, [q(0), q(1)].concat())); + quads.push(b.element(ElementKind::Quad4, q(1))); + } + } + b.section("slab", hexes, SectionKind::Solid {}); + b.section("skin", quads, SectionKind::Shell { thickness: skin, offset: 0.0, points: 5 }); + let strain = 1e-3; + b.hold("left", [Some(0.0), Some(0.0), Some(0.0), FREE, FREE, FREE], |p| p[0] == 0.0); + b.hold("right", [Some(strain * length), Some(0.0), Some(0.0), FREE, FREE, FREE], |p| p[0] == length); + let r = b.solve(); + let want = e * (width * depth + width * skin) * strain; + let pull: f64 = b.nodes_where(|p| p[0] == length).iter().map(|n| r.reaction[at(&r, *n)][0]).sum(); + eprintln!("solid with shell skin: end reaction {pull:.9e} vs {want:.9e}"); + assert!(rel(pull, want) < 1e-10, "{pull}"); + let top = r.shells[0].stress.middle[0][0]; + assert!(rel(top, e * strain) < 1e-10, "{top}"); +} + +// --------------------------------------------------------------------------------------------- +// Errors, not crashes +// --------------------------------------------------------------------------------------------- + +#[test] +fn structural_models_that_cannot_be_solved_are_errors() { + // A beam whose orientation is along it. + let (mut b, _, _) = cantilever(2, 100.0, BeamProfile::Circle { radius: 1.0 }, Procedure::Static {}); + b.model + .replace_section( + "beam", + Section { + set: "beam".into(), + material: "Mat".into(), + kind: SectionKind::Beam { profile: BeamProfile::Circle { radius: 1.0 }, orientation: [2.0, 2.0, 1.0] }, + formulation: Formulation::default(), + }, + ) + .unwrap(); + b.force("tip", vec![3], [1.0, 0.0, 0.0]); + let err = solve_static(&b.model, None, &Limits::default()).unwrap_err(); + assert!(err.to_string().contains("parallel"), "{err}"); + + // A quad8 shell is not supported yet; a degenerate quad4 is an error naming the element. + for (kind, nodes, want) in [(ElementKind::Quad8, 8, "not supported"), (ElementKind::Quad4, 4, "degenerate")] { + let mut b = Build::new(Elastic::Isotropic { young: 1.0, poisson: 0.0 }, 1.0, Procedure::Static {}); + for i in 0..nodes { + b.node(i as u32 + 1, [f64::from(i as u8 % 2), 0.0, 0.0]); + } + let el = b.element(kind, (1..=nodes as u32).collect()); + b.section("s", vec![el], SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }); + let err = solve_static(&b.model, None, &Limits::default()).unwrap_err(); + assert!(err.to_string().contains(want), "{kind:?}: {err}"); + } + + // A moment on the nodes of a rod: they have no rotations. + let mut b = Build::new(Elastic::Isotropic { young: 1.0, poisson: 0.0 }, 1.0, Procedure::Static {}); + b.node(1, [0.0; 3]); + b.node(2, [1.0, 0.0, 0.0]); + let el = b.element(ElementKind::Rod2, vec![1, 2]); + b.section("r", vec![el], SectionKind::Rod { area: 1.0 }); + b.hold_nodes("all", [Some(0.0); 6], vec![1]); + b.model.add_set("end", Set { nodes: vec![2], ..Set::default() }).unwrap(); + b.load("twist", Load::Moment { target: Target::Set("end".into()), vector: [1.0, 0.0, 0.0], distribution: Distribution::Total }); + let err = solve_static(&b.model, None, &Limits::default()).unwrap_err(); + assert!(err.to_string().contains("rotations"), "{err}"); +} diff --git a/docs/architecture.md b/docs/architecture.md index a348e0c..613dcc6 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -150,13 +150,20 @@ module docs; ADR 0001 records why the linear algebra is our own. the unknowns their nodes are given (`discrete::needs`, `discrete::dof_of`); dashpots are ignored with a warning (no velocities in static and frequency steps). Every reference point a coupling or a discrete element uses gets six degrees of freedom (translations, rotations) after - the mesh's. Shells - and beams return "not supported yet". The step's boundary conditions prescribe translations - (rotations do not exist on solid nodes) and a coupled point's six directions; its loads become - a load vector: forces and tractions from `Resolved::nodal_loads`, forces and moments on coupled - reference points, pressure as consistent face loads, gravity and body forces through the - element shape functions. -2. **Assembly** on the node graph (3 × 3 blocks for solids, 2 × 2 in 2D models, elements in label order, so sums are + the mesh's. Shell sections make quad4 (MITC4) and tri3 (MITC3) flat shells with a drilling + rotation (`modelcraft_fem::shell`; quad8 and tri6 shells are "not supported yet"), beam + sections two-node Timoshenko beams with the profile's section properties and the section's + orientation vector (`modelcraft_fem::beam`); both declare `ElementDofs::new(DofSet::STRUCTURAL)` + (`system::Structural`), so their nodes carry translations and rotations and the solids and + trusses sharing those nodes add into the translations. Shells and beams need isotropic + materials until material orientations exist, and are 3D only. The step's boundary conditions + prescribe the kinds a node carries (rotations of solid-only nodes are ignored) and a coupled + point's six directions; its loads become a load vector: forces, moments and tractions from + `Resolved::nodal_loads`, forces and moments on coupled reference points, pressure as + consistent face loads (shell faces too), gravity and body forces through the element shape + functions (consistent forces and moments on shells and beams). +2. **Assembly** on the node graph (one dense block per node pair: 3 × 3 for solids, 2 × 2 in 2D + models, up to 6 × 6 at shell and beam nodes; elements in label order, so sums are reproducible bit for bit). 3. **Ordering and elimination** (`constraints`). Nested dissection of the nodes by coordinate bisection. Couplings become linear constraint rows: kinematic = rigid link of each coupled @@ -170,8 +177,10 @@ module docs; ADR 0001 records why the linear algebra is our own. numbers are those of plain elimination, bit for bit. 4. **Static solution.** Up-looking sparse Cholesky (the factor size is computed and capped before any value), reactions T_pᵀ (K u − f) (at nodes and at prescribed reference points, forces and - moments), integration-point stresses extrapolated to the nodes and averaged, von Mises. A zero - pivot names the node or reference point and direction that can move freely. + moments), integration-point stresses extrapolated to the nodes and averaged, von Mises. Shells + report section forces and bottom, middle and top stresses at their nodes (their top surface + joins the nodal average), beams their end section forces; nodes report rotations and reaction + moments. A zero pivot names the node or reference point and direction that can move freely. 5. **Frequency.** Shift-invert block Lanczos (block 6, full reorthogonalization) on K − σM, σ = 0 for a supported structure and slightly negative for a free one, so rigid-body modes come out at zero; mass-normalized shapes, participation factors and effective masses. @@ -196,6 +205,12 @@ a kinematically coupled point (twist and shear stress). `crates/solver/tests/one forces of determinate trusses (a plane Pratt truss, a space tripod) against the method of joints, the three-bar truss, a rod3 constant-stress patch and a bar hanging under its weight, spring–mass, skew springs, the two-mass chain, a bar with a tip mass and the torsional pendulum. +`crates/solver/tests/structural.rs`: membrane and bending patch tests of quad4 and tri3 shells +in a tilted plane, the simply supported plate under pressure against Navier's series (deflection +and first frequency), cantilever beams against Timoshenko (exact) and Euler–Bernoulli +frequencies, a portal frame against slope-deflection, solids continued by rods and skinned by +shells. `crates/solver/tests/benchmarks.rs` (published reference values, tagged with their +source, separate from the analytic gate): the Scordelis–Lo roof and the pinched cylinder. ## Explicit solver diff --git a/docs/elements.md b/docs/elements.md index 9a59649..f729ba0 100644 --- a/docs/elements.md +++ b/docs/elements.md @@ -32,7 +32,7 @@ hex8/20/27 with first and second derivatives). |---|---|---| | Truss, 2-node (2D and 3D) | Axial members, cables (compression allowed) | done (rod2: implicit, in plane models with their two unknowns per node; explicit, corotational under nonlinear geometry; Abaqus T3D2 / T2D2, Nastran CROD, LS-DYNA truss beam (3D models)) | | Truss, 3-node (2D and 3D) | Curved axial members | done (rod3: implicit, plane models too; explicit with small displacements; Abaqus T3D3 / T2D3; no Nastran or LS-DYNA counterpart) | -| Beam, 2-node Timoshenko, 3D | Frames, stiffeners; shear deformable | in progress (shells and beams) | +| Beam, 2-node Timoshenko, 3D | Frames, stiffeners; shear deformable | done (exact interdependent interpolation, profile properties, orientation vector; `modelcraft_fem::beam`) | | Beam, 2-node Euler–Bernoulli (cubic), 2D and 3D | Slender frames; exact for end-loaded prismatic members | planned (1D elements II) | | Beam, 3-node Timoshenko (quadratic), 2D and 3D | Curved and stocky beams | planned (1D elements II) | | Beam, 2D (planar frames) | In-plane frames with 3 DOFs per node | planned (1D elements II) | @@ -57,7 +57,7 @@ hex8/20/27 with first and second derivatives). | Element | Use | Status | |---|---|---| -| Shell quad4 (MITC4), triangle (locking-free) | Thin and thick shells | in progress (shells and beams) | +| Shell quad4 (MITC4), triangle (locking-free) | Thin and thick shells | done (flat MITC4 and MITC3 with a Hughes–Brezzi drilling rotation; `modelcraft_fem::shell`) | | Shell quad8 / tri6 | Curved thick shells | planned | | Continuum shells (8-node, on swept hex) | Thick-to-thin shells, composites through thickness | planned (M8) | | Membranes (tri3, quad4, quad8) | Fabrics, inflatables, skins without bending | planned | diff --git a/docs/features.md b/docs/features.md index cbafa94..c903009 100644 --- a/docs/features.md +++ b/docs/features.md @@ -135,7 +135,7 @@ shown greyed with the roadmap milestone that brings them (ROADMAP.md). | Display | Contours — Contours of displacement and stress, with a legend | ready | `results.contour` | | Display | Mode Shapes — Mode shapes and frequencies of a frequency step | ready | `results.modes` | | Display | Animate — Animate the deformed shape, a mode or an explicit step's frames | ready | `results.animate`, `results.frame` | -| Query | Probe — Displacement, reaction and stress at nodes, a mode shape, or an explicit field frame (or click the model) | ready | `solve.probe`, `solve.history` | +| Query | Probe — Displacement, rotation, reaction and stress at nodes, element section forces, a mode shape, or an explicit field frame (or click the model) | ready | `solve.probe`, `solve.history` | | Query | Reactions — Reaction force and moment sums per set | ready | `results.reactions` | | 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 54cd3a6..35e1edb 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -62,7 +62,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | Decks | Nastran bulk data writer (`deck.write`, SOL 101 or 103 per deck; reports what Nastran cannot express) | M3 | | done | | Decks | LS-DYNA keyword writer (implicit; explicit in M8) | M3 | | done (implicit static and eigenvalue; beam sections planned) | | Decks | Readers keeping unknown keywords | M9 | | planned | -| Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, wedge6, hex8, hex20; 2D plane stress, plane strain and axisymmetric tri3, tri6, quad4 (full, incompatible modes), quad8 (full, reduced); reference points, kinematic and distributing couplings, moments; centrifugal loads) | +| Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, wedge6, hex8, hex20; 2D plane stress, plane strain and axisymmetric tri3, tri6, quad4 (full, incompatible modes), quad8 (full, reduced); quad4/tri3 MITC shells and Timoshenko beams with profiles, mixed with solids; reference points, kinematic and distributing couplings, moments; centrifugal loads) | | 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 and X2: large deformation, objective rates, J2 and Johnson–Cook plasticity, damage, erosion; 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 | diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index f87fec7..476b189 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -1003,6 +1003,24 @@ the module docs): Cook et al. 2002, §2.2 and §2.10; HRZ lumping as above. Verified against the method of joints (statically determinate trusses), the three-bar truss, the hanging bar, √(k/m), the two-mass chain and the torsional pendulum √(k_t/J). +- Shells: E. N. Dvorkin, K.-J. Bathe, *A continuum mechanics based four-node shell element for + general nonlinear analysis*, Engineering Computations 1(1):77–88, 1984 (MITC4 shear tying); + K.-J. Bathe, E. N. Dvorkin, IJNME 21(2):367–383, 1985 (the flat-plate form); P.-S. Lee, + K.-J. Bathe, *Development of MITC isotropic triangular shell finite elements*, Computers & + Structures 82(11–13):945–962, 2004 (MITC3); T. J. R. Hughes, F. Brezzi, *On drilling degrees of + freedom*, CMAME 72(1):105–121, 1989 (drilling penalty); facet shells, O. C. Zienkiewicz, + R. L. Taylor, D. D. Fox, *The Finite Element Method for Solid and Structural Mechanics*, 7th + ed., 2014, ch. 9–10. +- Beams: J. S. Przemieniecki, *Theory of Matrix Structural Analysis*, McGraw-Hill 1968, §5.6; + Z. Friedman, J. B. Kosmatka, Computers & Structures 47(3):473–481, 1993 (interdependent + Timoshenko interpolation); G. R. Cowper, J. Appl. Mech. 33(2):335–340, 1966 (shear + coefficients); S. P. Timoshenko, J. N. Goodier, *Theory of Elasticity*, 3rd ed., 1970, §109 + (torsion of a rectangle); W. D. Pilkey, *Analysis and Design of Elastic Beams*, Wiley 2002 + (thin-walled shear areas and torsion). +- Shell benchmarks (published values, tagged, not the analytic gate): Scordelis–Lo roof 0.3024 + (MacNeal & Harder 1985); pinched cylinder 1.8248 · 10⁻⁵ (T. Belytschko, H. Stolarski, W. K. + Liu, N. Carpenter, J. S.-J. Ong, CMAME 51(1–3):221–258, 1985). Plate: S. Timoshenko, + S. Woinowsky-Krieger, *Theory of Plates and Shells*, 2nd ed., 1959, §30 (Navier). - Verification: R. H. MacNeal, R. L. Harder, *A proposed standard set of problems to test finite element accuracy*, Finite Elem. Anal. Des. 1(1):3–20, 1985 (patch tests); Lamé, Kirsch, Timoshenko and Euler–Bernoulli closed forms; Saint-Venant's principle (tension through a