diff --git a/crates/explicit/src/elements.rs b/crates/explicit/src/elements.rs index a762570..f718f40 100644 --- a/crates/explicit/src/elements.rs +++ b/crates/explicit/src/elements.rs @@ -38,7 +38,7 @@ use modelcraft_fem::material::stress; use modelcraft_fem::uniform::UniformStrain; use modelcraft_model::{BulkViscosity, HourglassKind}; -use crate::system::{Material, System}; +use crate::system::{Material, NODE_DOFS, System}; /// Scratch vectors and the element state of the latest evaluation. #[derive(Debug, Clone)] @@ -82,13 +82,13 @@ impl Workspace { } pub(crate) fn gather(x: &[f64], nodes: &[usize; N]) -> [[f64; 3]; N] { - nodes.map(|n| [0, 1, 2].map(|k| x.get(3 * n + k).copied().unwrap_or(0.0))) + nodes.map(|n| [0, 1, 2].map(|k| x.get(NODE_DOFS * n + k).copied().unwrap_or(0.0))) } /// Add a node's vector into a per-degree-of-freedom vector. pub(crate) fn add_at(f: &mut [f64], node: usize, v: &[f64; 3]) { for (k, c) in v.iter().enumerate() { - if let Some(x) = f.get_mut(3 * node + k) { + if let Some(x) = f.get_mut(NODE_DOFS * node + k) { *x += c; } } diff --git a/crates/explicit/src/integrate.rs b/crates/explicit/src/integrate.rs index fabe939..77d1b57 100644 --- a/crates/explicit/src/integrate.rs +++ b/crates/explicit/src/integrate.rs @@ -24,7 +24,7 @@ use modelcraft_model::{ExplicitDynamic, Model, OutputVariable, Procedure, Step}; use crate::elements::{Workspace, internal_forces}; use crate::result::{ControllingElement, Energies, ExplicitResult, Frame, History}; -use crate::system::System; +use crate::system::{NODE_DOFS, System}; use crate::{Error, Limits, Result}; /// How many elements the result names as controlling the time increment. @@ -76,7 +76,7 @@ fn constrain(sys: &System, u: &[f64], v: &mut [f64], dt: f64) { for w in &sys.walls { for &node in &w.nodes { let x = sys.xyz.get(node).copied().unwrap_or_default(); - let at = |k: usize| 3 * node + k; + let at = |k: usize| NODE_DOFS * node + k; let gap: f64 = (0..3).map(|k| (x[k] + u.get(at(k)).copied().unwrap_or(0.0) - w.point[k]) * w.normal[k]).sum(); let normal: f64 = (0..3).map(|k| v.get(at(k)).copied().unwrap_or(0.0) * w.normal[k]).sum(); if gap + dt * normal < 0.0 { @@ -113,7 +113,7 @@ fn impulses(mass: &[f64], free: &[f64], constrained: &[f64], dt: f64, r: &mut [f } fn per_node(x: &[f64]) -> Vec<[f64; 3]> { - x.as_chunks::<3>().0.to_vec() + x.as_chunks::().0.to_vec() } /// Solve an explicit dynamic step: the step named `step`, or the model's only one. @@ -171,7 +171,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res } let n = sys.dofs(); - let mass: Vec = (0..n).map(|d| sys.mass.get(d / 3).copied().unwrap_or(0.0)).collect(); + let mass: Vec = (0..n).map(|d| sys.mass.get(d / NODE_DOFS).copied().unwrap_or(0.0)).collect(); if mass.iter().any(|m| !(*m > 0.0 && m.is_finite())) { return Err(Error::Internal("a node has no positive mass".into())); } @@ -209,7 +209,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res history.energies.push(e); let sum = |x: &[f64], weight: &dyn Fn(usize) -> f64| { let mut s = [0.0; 3]; - for (node, c) in x.as_chunks::<3>().0.iter().enumerate() { + for (node, c) in x.as_chunks::().0.iter().enumerate() { for k in 0..3 { s[k] += weight(node) * c[k]; } @@ -218,7 +218,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res }; history.momentum.push(sum(v, &|node| sys.mass.get(node).copied().unwrap_or(0.0))); history.constraint_force.push(sum(r, &|_| 1.0)); - let pick = |x: &[f64]| sys.history_nodes.iter().map(|i| [0, 1, 2].map(|k| x.get(3 * i + k).copied().unwrap_or(0.0))).collect(); + let pick = |x: &[f64]| sys.history_nodes.iter().map(|i| [0, 1, 2].map(|k| x.get(NODE_DOFS * i + k).copied().unwrap_or(0.0))).collect(); history.displacement.push(pick(u)); history.velocity.push(pick(v)); }; diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index 9d2ed7e..d61dc9b 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -2,14 +2,17 @@ //! lengths, lumped masses, prescribed velocities, rigid walls, initial velocities and the load //! vector of an explicit dynamic step. -use modelcraft_fem::face; use modelcraft_fem::material::Elasticity; use modelcraft_fem::uniform::{Hex8, Tet4, modulus_bound}; +use modelcraft_fem::{DofKind, DofLayout, DofSet, face}; use modelcraft_model::{Elastic, ElementKind, ExplicitDynamic, HourglassKind, Integration, Load, Model, Resolved, SectionKind, Step, Target}; use crate::elements::add_at; use crate::{Error, Result}; +/// Degrees of freedom per node: the three translations, node-major (see [`System`]). +pub const NODE_DOFS: usize = 3; + /// Most nodes whose histories a step records. pub const MAX_HISTORY_NODES: usize = 1000; @@ -67,10 +70,15 @@ pub struct Wall { } /// The discrete system of an explicit dynamic step. +/// +/// Its elements are solids, so every node carries the three translations and the layout is +/// node-major: degree of freedom 3 n + axis, which the integrator's per-node loops rely on +/// (checked when the system is built). #[derive(Debug, Clone)] pub struct System { /// Node labels with degrees of freedom (those of the elements), sorted. pub labels: Vec, + pub layout: DofLayout, pub xyz: Vec<[f64; 3]>, /// Element labels, sorted. pub elements: Vec, @@ -188,25 +196,30 @@ impl System { let mut labels: Vec = chosen.iter().flat_map(|e| e.2.iter().copied()).collect(); labels.sort_unstable(); labels.dedup(); - if 3 * labels.len() > max_dofs { - return Err(Error::Limit(format!("the model has {} degrees of freedom; the limit is {max_dofs}", 3 * labels.len()))); + let layout = DofLayout::uniform(labels.iter().copied(), DofSet::TRANSLATIONS); + if layout.dofs() > max_dofs { + return Err(Error::Limit(format!("the model has {} degrees of freedom; the limit is {max_dofs}", layout.dofs()))); + } + if layout.stride() != Some(NODE_DOFS) { + return Err(Error::Internal("explicit dynamics needs three translations per node".into())); } if labels.is_empty() { return Err(invalid("the model has no elements to analyse")); } let xyz: Vec<[f64; 3]> = labels.iter().map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| invalid(format!("node {l} is missing")))).collect::>()?; - let n = labels.len(); + let (n, dofs) = (labels.len(), layout.dofs()); let mut sys = System { + layout, elements: chosen.iter().map(|e| e.0).collect(), materials, hexes: Vec::new(), tets: Vec::new(), mass: vec![0.0; n], - prescribed: vec![None; 3 * n], + prescribed: vec![None; dofs], walls: Vec::new(), - initial_velocity: vec![0.0; 3 * n], - external: vec![0.0; 3 * n], + initial_velocity: vec![0.0; dofs], + external: vec![0.0; dofs], history_nodes: Vec::new(), hourglass: settings.hourglass.kind, warnings: resolved.warnings.iter().map(ToString::to_string).collect(), @@ -270,7 +283,7 @@ impl System { } pub fn dofs(&self) -> usize { - 3 * self.labels.len() + self.layout.dofs() } /// Prescribed velocities from the step's boundary conditions: zero for held degrees of @@ -283,13 +296,15 @@ impl System { } for &label in resolved.target_nodes(bc.target()) { let Some(node) = self.node(label) else { continue }; - for dof in 1..=3u8 { + for kind in self.layout.kinds(node).iter() { + let dof = kind.number(); let Some(value) = bc.value(dof) else { continue }; let velocity = value / end_time; if !velocity.is_finite() { return Err(invalid(format!("boundary condition `{name}`: the velocity {value} / {end_time} is not finite"))); } - let slot = self.prescribed.get_mut(3 * node + usize::from(dof - 1)).ok_or_else(|| Error::Internal("dof index".into()))?; + let index = self.layout.dof(node, kind).ok_or_else(|| Error::Internal("dof index".into()))?; + let slot = self.prescribed.get_mut(index).ok_or_else(|| Error::Internal("dof index".into()))?; match *slot { Some(old) if old != velocity => { return Err(invalid(format!( @@ -381,7 +396,9 @@ impl System { for iv in &settings.initial_velocities { for node in self.nodes_of(resolved, &iv.set) { for (k, c) in iv.velocity.iter().enumerate() { - if let Some(v) = self.initial_velocity.get_mut(3 * node + k) { + if let Some(v) = + DofKind::translation(k).and_then(|kind| self.layout.dof(node, kind)).and_then(|d| self.initial_velocity.get_mut(d)) + { *v = *c; } } diff --git a/crates/fem/src/basis.rs b/crates/fem/src/basis.rs new file mode 100644 index 0000000..5ee7c47 --- /dev/null +++ b/crates/fem/src/basis.rs @@ -0,0 +1,762 @@ +//! The shape-function library: every element family's interpolation with its first and second +//! derivatives in natural coordinates, and their mapping to physical gradients, Hessians and +//! Laplacians through the isoparametric map. +//! +//! Families ([`Basis`]) and their reference elements (natural coordinates ξ, η, ζ; unused ones +//! are zero): +//! - **lines** on ξ ∈ [−1, 1]: line2 (ends at ξ = −1, 1), line3 (then the middle node at 0); +//! - **triangles** on {ξ, η ≥ 0, ξ + η ≤ 1}: tri3 at (0,0), (1,0), (0,1), area coordinates +//! L₁ = 1 − ξ − η, L₂ = ξ, L₃ = η; tri6 adds mid-side nodes on edges 1-2, 2-3, 3-1 with +//! N = 4 Lₐ L_b and corners N = L(2L − 1); +//! - **quadrilaterals** on [−1, 1]²: corners counter-clockwise from (−1, −1); quad4 bilinear, +//! quad8 serendipity (mid-side nodes on edges 1-2, 2-3, 3-4, 4-1), quad9 biquadratic Lagrange +//! (quad8 plus the centre); +//! - **tetrahedra**: tet4, tet10 (see [`crate::shape`]); +//! - **wedges**: the triangle times ζ ∈ [−1, 1]; wedge6, and wedge15 with mid-side nodes on the +//! bottom triangle's edges, the top's, then the vertical edges 1-4, 2-5, 3-6 (the order of +//! [`ElementKind::Wedge6`]'s edges): corners N = ½ L(2L − 1)(1 + ζζᵢ) − ½ L(1 − ζ²), triangle +//! mid-sides N = 2 Lₐ L_b (1 + ζζᵢ), vertical mid-sides N = L(1 − ζ²); +//! - **pyramid5**: square base [−1, 1]² at ζ = 0 (corners as quad4), apex (0, 0, 1); the rational +//! functions N = ¼[(1 − ζ) + ξξᵢ + ηηᵢ + ξᵢηᵢ ξη / (1 − ζ)] and N₅ = ζ, linearly complete and +//! conforming with neighbouring tetrahedra and hexahedra (L. Bedrosian, *Shape functions and +//! integration formulas for three-dimensional finite element analysis*, IJNME 35(1):95–108, +//! 1992). At the apex itself the rational term and its derivatives take their limit along the +//! axis, zero; quadrature points never lie there; +//! - **hexahedra** on [−1, 1]³: hex8 trilinear, hex20 serendipity (see [`crate::shape`]), hex27 +//! triquadratic Lagrange (hex20, then the face centres in the order of [`ElementKind::Hex8`]'s +//! faces, then the centre). +//! +//! Node orders are corners first, then mid-side nodes in the model's edge order (the model is the +//! one place they are written down), then face and body nodes, the Abaqus/VTK convention of +//! [`modelcraft_model::ElementKind`]. The functions are the textbook ones (O. C. Zienkiewicz, +//! R. L. Taylor, J. Z. Zhu, *The Finite Element Method: Its Basis and Fundamentals*, 7th ed., +//! 2013, ch. 6; T. J. R. Hughes, *The Finite Element Method*, Dover 2000, ch. 3). +//! +//! **Derivatives.** [`Basis::eval`] gives values and first derivatives; the families the solid +//! kernels used before the library existed (tet4, tet10, wedge6, hex8, hex20) keep their +//! hand-written functions there, so element results do not change by a bit. [`Basis::eval2`] +//! adds the Hessian ∂²N/∂ξₐ∂ξ_b, computed from the same formulas by second-order automatic +//! differentiation ([`crate::jet`]): exact to round-off for every family. +//! +//! **Physical derivatives** (Hughes 2000, §3.9 and §4.2 for the second derivatives): with +//! J_ab = ∂x_b/∂ξₐ, the gradient is ∇ₓN = J⁻¹ ∇_ξN and the Hessian +//! Hₓ = J⁻¹ (H_ξ − Σ_b (∂N/∂x_b) X_b″) J⁻ᵀ, where X_b″ = Σᵢ x_{i,b} H_ξ(Nᵢ) is the second +//! derivative of the map. On a curved (non-affine) element that last term is what makes a field +//! linear in x have zero Hessian; dropping it is a classic error in residual-based stabilization. +//! A family of dimension d maps onto the first d physical coordinates (plane elements in x–y, +//! lines along x); elements embedded in a higher dimension (shells, beams, trusses in 3D) pass +//! coordinates in their own local frame. + +use std::sync::LazyLock; + +use modelcraft_model::ElementKind; + +use crate::dense::positive_inverse; +use crate::jet::Jet; +use crate::shape::{self, HEX_CORNERS, MAX_NODES, Shape, TET_CORNERS, WEDGE_CORNERS}; +use crate::{Error, Result}; + +/// An element family's interpolation. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub enum Basis { + Line2, + Line3, + Tri3, + Tri6, + Quad4, + Quad8, + Quad9, + Tet4, + Tet10, + Wedge6, + Wedge15, + Pyramid5, + Hex8, + Hex20, + Hex27, +} + +/// Values, first and second natural derivatives at one point. +#[derive(Debug, Clone, Copy)] +pub struct Shape2 { + pub shape: Shape, + /// ∂²Nᵢ/∂ξₐ∂ξ_b. + pub d2n: [[[f64; 3]; 3]; MAX_NODES], +} + +impl Shape2 { + pub fn hessians(&self) -> &[[[f64; 3]; 3]] { + self.d2n.get(..self.shape.count).unwrap_or_default() + } +} + +/// Physical gradients at one point. +#[derive(Debug, Clone, Copy)] +pub struct Gradients { + pub count: usize, + /// ∂Nᵢ/∂x_b. + pub dx: [[f64; 3]; MAX_NODES], + /// det J (a length, area or volume ratio for dimension 1, 2, 3). + pub det: f64, +} + +impl Gradients { + pub fn gradients(&self) -> &[[f64; 3]] { + self.dx.get(..self.count).unwrap_or_default() + } +} + +/// Physical gradients and Hessians at one point. +#[derive(Debug, Clone, Copy)] +pub struct Hessians { + pub gradients: Gradients, + /// ∂²Nᵢ/∂x_a∂x_b. + pub d2x: [[[f64; 3]; 3]; MAX_NODES], +} + +impl Hessians { + pub fn hessians(&self) -> &[[[f64; 3]; 3]] { + self.d2x.get(..self.gradients.count).unwrap_or_default() + } + + /// ∇²Nᵢ, the trace of the Hessian. + pub fn laplacian(&self, i: usize) -> f64 { + self.d2x.get(i).map_or(0.0, |h| h[0][0] + h[1][1] + h[2][2]) + } +} + +const LINE_ENDS: [[f64; 3]; 2] = [[-1., 0., 0.], [1., 0., 0.]]; +const TRI_CORNERS: [[f64; 3]; 3] = [[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]]; +const QUAD_CORNERS: [[f64; 3]; 4] = [[-1., -1., 0.], [1., -1., 0.], [1., 1., 0.], [-1., 1., 0.]]; +const PYRAMID_CORNERS: [[f64; 3]; 5] = [[-1., -1., 0.], [1., -1., 0.], [1., 1., 0.], [-1., 1., 0.], [0., 0., 1.]]; + +fn midpoint(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [0, 1, 2].map(|k| (a[k] + b[k]) / 2.0) +} + +/// Corners, then the midpoints of `edges` (pairs of corner positions). +fn with_midsides(corners: &[[f64; 3]], edges: &[&[usize]]) -> Vec<[f64; 3]> { + let mut v = corners.to_vec(); + for e in edges { + if let (Some(a), Some(b)) = (e.first().and_then(|a| corners.get(*a)), e.get(1).and_then(|b| corners.get(*b))) { + v.push(midpoint(*a, *b)); + } + } + v +} + +/// One-dimensional quadratic Lagrange function of the node at t = c ∈ {−1, 0, 1}. +fn lagrange2(t: Jet, c: f64) -> Jet { + if c < 0.0 { + t * (t - 1.0) * 0.5 + } else if c > 0.0 { + t * (t + 1.0) * 0.5 + } else { + 1.0 - t * t + } +} + +impl Basis { + pub const ALL: [Basis; 15] = [ + Basis::Line2, + Basis::Line3, + Basis::Tri3, + Basis::Tri6, + Basis::Quad4, + Basis::Quad8, + Basis::Quad9, + Basis::Tet4, + Basis::Tet10, + Basis::Wedge6, + Basis::Wedge15, + Basis::Pyramid5, + Basis::Hex8, + Basis::Hex20, + Basis::Hex27, + ]; + + pub fn name(self) -> &'static str { + match self { + Basis::Line2 => "line2", + Basis::Line3 => "line3", + Basis::Tri3 => "tri3", + Basis::Tri6 => "tri6", + Basis::Quad4 => "quad4", + Basis::Quad8 => "quad8", + Basis::Quad9 => "quad9", + Basis::Tet4 => "tet4", + Basis::Tet10 => "tet10", + Basis::Wedge6 => "wedge6", + Basis::Wedge15 => "wedge15", + Basis::Pyramid5 => "pyramid5", + Basis::Hex8 => "hex8", + Basis::Hex20 => "hex20", + Basis::Hex27 => "hex27", + } + } + + /// The interpolation of a model element kind (line elements: along their axis). + pub fn of(kind: ElementKind) -> Basis { + match kind { + ElementKind::Rod2 | ElementKind::Beam2 => Basis::Line2, + ElementKind::Tri3 => Basis::Tri3, + ElementKind::Tri6 => Basis::Tri6, + ElementKind::Quad4 => Basis::Quad4, + ElementKind::Quad8 => Basis::Quad8, + ElementKind::Tet4 => Basis::Tet4, + ElementKind::Tet10 => Basis::Tet10, + ElementKind::Wedge6 => Basis::Wedge6, + ElementKind::Hex8 => Basis::Hex8, + ElementKind::Hex20 => Basis::Hex20, + } + } + + /// Dimension of the reference element. + pub fn dim(self) -> usize { + match self { + Basis::Line2 | Basis::Line3 => 1, + Basis::Tri3 | Basis::Tri6 | Basis::Quad4 | Basis::Quad8 | Basis::Quad9 => 2, + _ => 3, + } + } + + /// Number of nodes (and functions). + pub fn count(self) -> usize { + match self { + Basis::Line2 => 2, + Basis::Line3 | Basis::Tri3 => 3, + Basis::Quad4 | Basis::Tet4 => 4, + Basis::Pyramid5 => 5, + Basis::Tri6 | Basis::Wedge6 => 6, + Basis::Quad8 | Basis::Hex8 => 8, + Basis::Quad9 => 9, + Basis::Tet10 => 10, + Basis::Wedge15 => 15, + Basis::Hex20 => 20, + Basis::Hex27 => 27, + } + } + + /// The corner (vertex) nodes, first in the node list. + pub fn corners(self) -> usize { + match self { + Basis::Line2 | Basis::Line3 => 2, + Basis::Tri3 | Basis::Tri6 => 3, + Basis::Quad4 | Basis::Quad8 | Basis::Quad9 | Basis::Tet4 | Basis::Tet10 => 4, + Basis::Pyramid5 => 5, + Basis::Wedge6 | Basis::Wedge15 => 6, + Basis::Hex8 | Basis::Hex20 | Basis::Hex27 => 8, + } + } + + /// The degree of the complete polynomials the family reproduces exactly (in natural + /// coordinates, and in physical ones on an affine map). + pub fn degree(self) -> usize { + match self { + Basis::Line2 | Basis::Tri3 | Basis::Quad4 | Basis::Tet4 | Basis::Wedge6 | Basis::Pyramid5 | Basis::Hex8 => 1, + _ => 2, + } + } + + /// The natural coordinates of the nodes, in node order. + pub fn nodes(self) -> &'static [[f64; 3]] { + static TABLES: LazyLock>> = LazyLock::new(|| Basis::ALL.iter().map(|b| b.build_nodes()).collect()); + TABLES.get(self as usize).map(Vec::as_slice).unwrap_or_default() + } + + fn build_nodes(self) -> Vec<[f64; 3]> { + match self { + Basis::Line2 => LINE_ENDS.to_vec(), + Basis::Line3 => with_midsides(&LINE_ENDS, ElementKind::Rod2.edges()), + Basis::Tri3 => TRI_CORNERS.to_vec(), + Basis::Tri6 => with_midsides(&TRI_CORNERS, ElementKind::Tri6.edges()), + Basis::Quad4 => QUAD_CORNERS.to_vec(), + Basis::Quad8 => with_midsides(&QUAD_CORNERS, ElementKind::Quad8.edges()), + Basis::Quad9 => { + let mut v = with_midsides(&QUAD_CORNERS, ElementKind::Quad8.edges()); + v.push([0.0; 3]); + v + } + Basis::Tet4 => TET_CORNERS.to_vec(), + Basis::Tet10 => with_midsides(&TET_CORNERS, ElementKind::Tet10.edges()), + Basis::Wedge6 => WEDGE_CORNERS.to_vec(), + Basis::Wedge15 => with_midsides(&WEDGE_CORNERS, ElementKind::Wedge6.edges()), + Basis::Pyramid5 => PYRAMID_CORNERS.to_vec(), + Basis::Hex8 => HEX_CORNERS.to_vec(), + Basis::Hex20 => with_midsides(&HEX_CORNERS, ElementKind::Hex20.edges()), + Basis::Hex27 => { + let mut v = with_midsides(&HEX_CORNERS, ElementKind::Hex20.edges()); + for face in ElementKind::Hex8.faces() { + let c = face.iter().filter_map(|i| HEX_CORNERS.get(*i)).fold([0.0; 3], |s, x| [0, 1, 2].map(|k| s[k] + x[k] / 4.0)); + v.push(c); + } + v.push([0.0; 3]); + v + } + } + } + + /// Values and first natural derivatives at `xi`. + pub fn eval(self, xi: [f64; 3]) -> Shape { + match self { + Basis::Tet4 => shape::tet4(xi), + Basis::Tet10 => shape::tet10(xi), + Basis::Wedge6 => shape::wedge6(xi), + Basis::Hex8 => shape::hex8(xi), + Basis::Hex20 => shape::hex20(xi), + _ => self.eval2(xi).shape, + } + } + + /// Values, first and second natural derivatives at `xi`. + pub fn eval2(self, xi: [f64; 3]) -> Shape2 { + let jets = self.jets(xi); + let count = self.count(); + let mut out = Shape2 { shape: Shape { count, n: [0.0; MAX_NODES], dn: [[0.0; 3]; MAX_NODES] }, d2n: [[[0.0; 3]; 3]; MAX_NODES] }; + for (i, j) in jets.iter().enumerate().take(count) { + out.shape.n[i] = j.v; + out.shape.dn[i] = j.g; + out.d2n[i] = j.h; + } + out + } + + /// The functions as jets in the natural coordinates (only the first `count` are set). + fn jets(self, xi: [f64; 3]) -> [Jet; MAX_NODES] { + let x = Jet::point(xi); + let dim = self.dim(); + let mut n = [Jet::constant(0.0); MAX_NODES]; + let nodes = self.nodes(); + // Area or volume coordinates of the simplex families. + let tri = [1.0 - x[0] - x[1], x[0], x[1]]; + let tet = [1.0 - x[0] - x[1] - x[2], x[0], x[1], x[2]]; + let set = |n: &mut [Jet; MAX_NODES], i: usize, v: Jet| { + if let Some(slot) = n.get_mut(i) { + *slot = v; + } + }; + match self { + Basis::Line2 | Basis::Quad4 | Basis::Hex8 => { + for (i, c) in nodes.iter().enumerate() { + let v = (0..dim).fold(Jet::constant(1.0), |p, k| p * (1.0 + x[k] * c[k]) * 0.5); + set(&mut n, i, v); + } + } + Basis::Line3 | Basis::Quad9 | Basis::Hex27 => { + for (i, c) in nodes.iter().enumerate() { + let v = (0..dim).fold(Jet::constant(1.0), |p, k| p * lagrange2(x[k], c[k])); + set(&mut n, i, v); + } + } + Basis::Quad8 | Basis::Hex20 => { + // Corners ∏(1 + ξc)/2ᵈ (Σ ξc − (d − 1)); a mid-side node with c_k = 0: + // (1 − ξ_k²) ∏_{l≠k}(1 + ξ_l c_l)/2ᵈ⁻¹. + for (i, c) in nodes.iter().enumerate() { + let along = (0..dim).find(|k| c[*k] == 0.0); + let v = match along { + None => { + let p = (0..dim).fold(Jet::constant(1.0), |p, k| p * (1.0 + x[k] * c[k]) * 0.5); + let s = (0..dim).fold(Jet::constant(-(dim as f64 - 1.0)), |s, k| s + x[k] * c[k]); + p * s + } + Some(a) => (0..dim).filter(|k| *k != a).fold(1.0 - x[a] * x[a], |p, k| p * (1.0 + x[k] * c[k]) * 0.5), + }; + set(&mut n, i, v); + } + } + Basis::Tri3 => { + for (i, l) in tri.iter().enumerate() { + set(&mut n, i, *l); + } + } + Basis::Tet4 => { + for (i, l) in tet.iter().enumerate() { + set(&mut n, i, *l); + } + } + Basis::Tri6 | Basis::Tet10 => { + let (l, edges): (&[Jet], _) = + if self == Basis::Tri6 { (&tri, ElementKind::Tri6.edges()) } else { (&tet, ElementKind::Tet10.edges()) }; + for (i, li) in l.iter().enumerate() { + set(&mut n, i, *li * (*li * 2.0 - 1.0)); + } + for e in edges { + if let [a, b, m] = **e + && let (Some(la), Some(lb)) = (l.get(a), l.get(b)) + { + set(&mut n, m, *la * *lb * 4.0); + } + } + } + Basis::Wedge6 => { + for i in 0..3 { + set(&mut n, i, tri[i] * (1.0 - x[2]) * 0.5); + set(&mut n, i + 3, tri[i] * (1.0 + x[2]) * 0.5); + } + } + Basis::Wedge15 => { + let bubble = 1.0 - x[2] * x[2]; + for i in 0..6 { + let (l, z) = (tri[i % 3], if i < 3 { -1.0 } else { 1.0 }); + set(&mut n, i, l * (l * 2.0 - 1.0) * (1.0 + x[2] * z) * 0.5 - l * bubble * 0.5); + } + for (e, edge) in ElementKind::Wedge6.edges().iter().enumerate() { + let [a, b] = **edge else { continue }; + let v = if e < 6 { + let z = if a < 3 { -1.0 } else { 1.0 }; + tri[a % 3] * tri[b % 3] * (1.0 + x[2] * z) * 2.0 + } else { + tri[a % 3] * bubble + }; + set(&mut n, 6 + e, v); + } + } + Basis::Pyramid5 => { + let d = 1.0 - x[2]; + let r = if d.v.abs() > 0.0 { x[0] * x[1] / d } else { Jet::constant(0.0) }; + for (i, c) in PYRAMID_CORNERS.iter().take(4).enumerate() { + set(&mut n, i, (d + x[0] * c[0] + x[1] * c[1] + r * (c[0] * c[1])) * 0.25); + } + set(&mut n, 4, x[2]); + } + } + n + } + + /// The Jacobian J_ab = ∂x_b/∂ξₐ over the family's dimensions, the identity in the others. + fn jacobian(self, s: &Shape, xyz: &[[f64; 3]]) -> Result<[[f64; 3]; 3]> { + if xyz.len() < s.count { + return Err(Error::Invalid(format!("a {} element has {} nodes, got {} coordinates", self.name(), s.count, xyz.len()))); + } + let dim = self.dim(); + let mut j = [[0.0; 3]; 3]; + for (d, x) in s.derivatives().iter().zip(xyz) { + for a in 0..dim { + for b in 0..dim { + j[a][b] += d[a] * x[b]; + } + } + } + for a in dim..3 { + j[a][a] = 1.0; + } + Ok(j) + } + + fn degenerate(self) -> Error { + Error::Degenerate(format!("a {} element is inverted or degenerate (its Jacobian is not positive)", self.name())) + } + + /// Physical gradients at a point from its natural shape (`eval`) and the node coordinates. + pub fn gradients(self, s: &Shape, xyz: &[[f64; 3]]) -> Result { + self.mapped(s, xyz).map(|m| m.0) + } + + /// The physical gradients and J⁻¹. + fn mapped(self, s: &Shape, xyz: &[[f64; 3]]) -> Result<(Gradients, [[f64; 3]; 3])> { + let (inv, det) = positive_inverse(&self.jacobian(s, xyz)?).ok_or_else(|| self.degenerate())?; + let mut dx = [[0.0; 3]; MAX_NODES]; + for (g, d) in dx.iter_mut().zip(s.derivatives()) { + for b in 0..3 { + g[b] = inv[b][0] * d[0] + inv[b][1] * d[1] + inv[b][2] * d[2]; + } + } + Ok((Gradients { count: s.count, dx, det }, inv)) + } + + /// Physical gradients and Hessians from a point's `eval2` and the node coordinates, with the + /// second-derivative term of the map. + pub fn hessians(self, s: &Shape2, xyz: &[[f64; 3]]) -> Result { + let (gradients, inv) = self.mapped(&s.shape, xyz)?; + let dim = self.dim(); + // X_b″ = Σᵢ x_{i,b} H_ξ(Nᵢ) for the mapped coordinates b < dim. + let mut map2 = [[[0.0; 3]; 3]; 3]; + for (h, x) in s.hessians().iter().zip(xyz) { + for b in 0..dim { + for p in 0..3 { + for q in 0..3 { + map2[b][p][q] += x[b] * h[p][q]; + } + } + } + } + let mut d2x = [[[0.0; 3]; 3]; MAX_NODES]; + for ((out, h), g) in d2x.iter_mut().zip(s.hessians()).zip(gradients.gradients()) { + // A = H_ξ − Σ_b g_b X_b″, then Hₓ = J⁻¹ A J⁻ᵀ. + let mut a = *h; + for b in 0..dim { + for p in 0..3 { + for q in 0..3 { + a[p][q] -= g[b] * map2[b][p][q]; + } + } + } + for r in 0..3 { + for c in 0..3 { + let mut v = 0.0; + for p in 0..3 { + for q in 0..3 { + v += inv[r][p] * a[p][q] * inv[c][q]; + } + } + out[r][c] = v; + } + } + } + Ok(Hessians { gradients, d2x }) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// A point inside every family's reference element. + const P: [f64; 3] = [0.21, 0.17, 0.33]; + + fn point(b: Basis) -> [f64; 3] { + match b.dim() { + 1 => [P[0], 0.0, 0.0], + 2 => [P[0], P[1], 0.0], + _ => P, + } + } + + /// Monomials ξᵃ ηᵇ ζᶜ of total degree ≤ `degree` in the family's dimensions, with their exact + /// gradient and Hessian. + fn monomials(dim: usize, degree: usize) -> Vec<[usize; 3]> { + let mut out = Vec::new(); + for a in 0..=degree { + for b in 0..=degree { + for c in 0..=degree { + let e = [a, b, c]; + if a + b + c <= degree && (dim..3).all(|k| e[k] == 0) { + out.push(e); + } + } + } + } + out + } + + fn monomial(e: [usize; 3], x: [f64; 3]) -> (f64, [f64; 3], [[f64; 3]; 3]) { + let p = |k: usize, d: usize| -> f64 { + // d-th derivative of x_k^{e_k}. + let n = e[k]; + if d > n { + return 0.0; + } + let coef: f64 = (0..d).map(|i| (n - i) as f64).product(); + coef * x[k].powi((n - d) as i32) + }; + let v = p(0, 0) * p(1, 0) * p(2, 0); + let mut g = [0.0; 3]; + let mut h = [[0.0; 3]; 3]; + for a in 0..3 { + let mut d = [0usize; 3]; + d[a] += 1; + g[a] = p(0, d[0]) * p(1, d[1]) * p(2, d[2]); + for b in 0..3 { + let mut dd = d; + dd[b] += 1; + h[a][b] = p(0, dd[0]) * p(1, dd[1]) * p(2, dd[2]); + } + } + (v, g, h) + } + + #[test] + fn kronecker_delta_and_partition_of_unity() { + for (i, b) in Basis::ALL.into_iter().enumerate() { + assert_eq!(b as usize, i, "Basis::ALL is in declaration order (node tables)"); + let nodes = b.nodes(); + assert_eq!(nodes.len(), b.count(), "{}", b.name()); + for (i, p) in nodes.iter().enumerate() { + for s in [b.eval(*p), b.eval2(*p).shape] { + for (j, v) in s.values().iter().enumerate() { + assert!((v - if i == j { 1.0 } else { 0.0 }).abs() < 1e-14, "{} N{j} at node {i}: {v}", b.name()); + } + } + } + let s = b.eval2(point(b)); + assert!((s.shape.values().iter().sum::() - 1.0).abs() < 1e-14, "{}", b.name()); + for a in 0..3 { + assert!(s.shape.derivatives().iter().map(|d| d[a]).sum::().abs() < 1e-13, "{}", b.name()); + for c in 0..3 { + assert!(s.hessians().iter().map(|h| h[a][c]).sum::().abs() < 1e-12, "{}", b.name()); + } + } + } + } + + #[test] + fn complete_polynomials_and_their_derivatives_are_reproduced() { + for b in Basis::ALL { + let nodes = b.nodes(); + let xi = point(b); + let s = b.eval2(xi); + for e in monomials(b.dim(), b.degree()) { + let (v, g, h) = monomial(e, xi); + let at = |i: usize| monomial(e, nodes[i]).0; + let iv: f64 = (0..b.count()).map(|i| s.shape.n[i] * at(i)).sum(); + assert!((iv - v).abs() < 1e-13, "{} {e:?}: {iv} vs {v}", b.name()); + for a in 0..3 { + let ig: f64 = (0..b.count()).map(|i| s.shape.dn[i][a] * at(i)).sum(); + assert!((ig - g[a]).abs() < 1e-12, "{} {e:?} d{a}: {ig} vs {}", b.name(), g[a]); + for c in 0..3 { + let ih: f64 = (0..b.count()).map(|i| s.d2n[i][a][c] * at(i)).sum(); + assert!((ih - h[a][c]).abs() < 1e-11, "{} {e:?} d{a}{c}: {ih} vs {}", b.name(), h[a][c]); + } + } + } + } + } + + #[test] + fn derivatives_match_finite_differences() { + let step = 1e-5; + for b in Basis::ALL { + let xi = point(b); + let s = b.eval2(xi); + let fast = b.eval(xi); + for i in 0..b.count() { + for a in 0..3 { + assert!((fast.dn[i][a] - s.shape.dn[i][a]).abs() < 1e-14, "{} eval vs eval2", b.name()); + } + assert!((fast.n[i] - s.shape.n[i]).abs() < 1e-15); + } + for a in 0..b.dim() { + let (mut p, mut m) = (xi, xi); + p[a] += step; + m[a] -= step; + let (sp, sm) = (b.eval2(p), b.eval2(m)); + for i in 0..b.count() { + let fd = (sp.shape.n[i] - sm.shape.n[i]) / (2.0 * step); + assert!((fd - s.shape.dn[i][a]).abs() < 1e-8, "{} dN{i}/dξ{a}", b.name()); + for c in 0..3 { + let fd2 = (sp.shape.dn[i][c] - sm.shape.dn[i][c]) / (2.0 * step); + assert!((fd2 - s.d2n[i][a][c]).abs() < 1e-7, "{} d²N{i}/dξ{a}dξ{c}: {fd2} vs {}", b.name(), s.d2n[i][a][c]); + } + } + } + } + } + + /// An affine map and a curved one (the affine map plus a smooth bend), applied to the + /// reference nodes. + fn affine(x: [f64; 3]) -> [f64; 3] { + [2.0 * x[0] + 0.3 * x[1] - 0.2 * x[2] + 1.0, 0.4 * x[0] + 1.5 * x[1] + 0.1 * x[2] - 2.0, -0.3 * x[0] + 0.2 * x[1] + 1.2 * x[2] + 0.5] + } + + fn curved(x: [f64; 3]) -> [f64; 3] { + let a = affine(x); + [a[0] + 0.15 * x[1] * x[1], a[1] + 0.1 * x[0] * x[0], a[2] + 0.12 * x[0] * x[1]] + } + + /// A quadratic field in physical coordinates (in the family's dimensions) with its gradient + /// and Hessian. + fn field(dim: usize, x: [f64; 3], degree: usize) -> (f64, [f64; 3], [[f64; 3]; 3]) { + let c = [0.7, -1.3, 0.4]; + let mut v = 0.5; + let mut g = [0.0; 3]; + let mut h = [[0.0; 3]; 3]; + for k in 0..dim { + v += c[k] * x[k]; + g[k] = c[k]; + } + if degree >= 2 { + let q = [[1.1, 0.3, -0.4], [0.3, -0.6, 0.25], [-0.4, 0.25, 0.9]]; + for a in 0..dim { + for b in 0..dim { + v += 0.5 * q[a][b] * x[a] * x[b]; + g[a] += q[a][b] * x[b]; + h[a][b] = q[a][b]; + } + } + } + (v, g, h) + } + + fn mapped(b: Basis, map: fn([f64; 3]) -> [f64; 3]) -> Vec<[f64; 3]> { + b.nodes() + .iter() + .map(|p| { + let x = map(*p); + [0, 1, 2].map(|k| if k < b.dim() { x[k] } else { 0.0 }) + }) + .collect() + } + + #[test] + fn physical_gradients_hessians_and_laplacians() { + for b in Basis::ALL { + let xi = point(b); + let s = b.eval2(xi); + for (map, degree) in [(affine as fn([f64; 3]) -> [f64; 3], b.degree()), (curved, 1)] { + let xyz = mapped(b, map); + let h = b.hessians(&s, &xyz).unwrap(); + let x: [f64; 3] = (0..b.count()).fold([0.0; 3], |acc, i| [0, 1, 2].map(|k| acc[k] + s.shape.n[i] * xyz[i][k])); + let (_, g, hh) = field(b.dim(), x, degree); + let nodal: Vec = xyz.iter().map(|p| field(b.dim(), *p, degree).0).collect(); + let grad: [f64; 3] = [0, 1, 2].map(|a| (0..b.count()).map(|i| h.gradients.dx[i][a] * nodal[i]).sum()); + for a in 0..3 { + assert!((grad[a] - g[a]).abs() < 1e-11, "{} grad {a}: {} vs {}", b.name(), grad[a], g[a]); + for c in 0..3 { + let v: f64 = (0..b.count()).map(|i| h.d2x[i][a][c] * nodal[i]).sum(); + assert!((v - hh[a][c]).abs() < 1e-10, "{} hess {a}{c}: {v} vs {}", b.name(), hh[a][c]); + } + } + let lap: f64 = (0..b.count()).map(|i| h.laplacian(i) * nodal[i]).sum(); + assert!((lap - (hh[0][0] + hh[1][1] + hh[2][2])).abs() < 1e-10, "{} laplacian", b.name()); + } + } + } + + #[test] + fn curved_map_hessians_match_finite_differences_in_physical_space() { + // On a curved element a quadratic field is not reproduced, but the physical Hessian of + // the interpolant must still be the derivative of its physical gradient. Differentiate + // the gradient along the natural directions and map: ∂(∇ₓu)/∂ξₐ = Σ_b J_ab Hₓ[b]. + for b in Basis::ALL.into_iter().filter(|b| b.degree() == 2) { + let xyz = mapped(b, curved); + let nodal: Vec = xyz.iter().map(|p| field(b.dim(), *p, 2).0 + 0.3 * p[0] * p[0] * p[0]).collect(); + let grad_at = |xi: [f64; 3]| -> [f64; 3] { + let g = b.gradients(&b.eval(xi), &xyz).unwrap(); + [0, 1, 2].map(|a| (0..b.count()).map(|i| g.dx[i][a] * nodal[i]).sum()) + }; + let xi = point(b); + let s = b.eval2(xi); + let h = b.hessians(&s, &xyz).unwrap(); + let hu: [[f64; 3]; 3] = [0, 1, 2].map(|r| [0, 1, 2].map(|c| (0..b.count()).map(|i| h.d2x[i][r][c] * nodal[i]).sum())); + let j = b.jacobian(&s.shape, &xyz).unwrap(); + let step = 1e-5; + for a in 0..b.dim() { + let (mut p, mut m) = (xi, xi); + p[a] += step; + m[a] -= step; + let (gp, gm) = (grad_at(p), grad_at(m)); + for c in 0..b.dim() { + let fd = (gp[c] - gm[c]) / (2.0 * step); + let exact: f64 = (0..3).map(|k| j[a][k] * hu[k][c]).sum(); + assert!((fd - exact).abs() < 1e-6, "{} ∂(∂u/∂x{c})/∂ξ{a}: {fd} vs {exact}", b.name()); + } + } + } + } + + #[test] + fn inverted_and_short_inputs_are_errors() { + let xyz = mapped(Basis::Quad4, affine); + let s = Basis::Quad4.eval([0.0; 3]); + assert!(Basis::Quad4.gradients(&s, &xyz[..3]).is_err()); + let flipped: Vec<[f64; 3]> = xyz.iter().map(|p| [p[0], -p[1], 0.0]).collect(); + assert!(matches!(Basis::Quad4.gradients(&s, &flipped), Err(Error::Degenerate(_)))); + let collapsed = vec![[0.0; 3]; 4]; + assert!(Basis::Quad4.gradients(&s, &collapsed).is_err()); + // The pyramid's apex is finite. + let apex = Basis::Pyramid5.eval2([0.0, 0.0, 1.0]); + assert!(apex.shape.values().iter().chain(apex.shape.derivatives().iter().flatten()).all(|v| v.is_finite())); + assert_eq!(Basis::of(ElementKind::Hex20), Basis::Hex20); + assert_eq!(Basis::of(ElementKind::Rod2).dim(), 1); + } +} diff --git a/crates/fem/src/dense.rs b/crates/fem/src/dense.rs index 307161c..294aa6b 100644 --- a/crates/fem/src/dense.rs +++ b/crates/fem/src/dense.rs @@ -105,6 +105,14 @@ pub fn inverse3(m: &[[f64; 3]; 3]) -> Option<([[f64; 3]; 3], f64)> { Some((inv, det)) } +/// The inverse and determinant of a Jacobian that maps a reference element onto a valid one: +/// `None` unless det J is positive beyond round-off relative to the product of its row norms. +pub fn positive_inverse(j: &[[f64; 3]; 3]) -> Option<([[f64; 3]; 3], f64)> { + let (inv, det) = inverse3(j)?; + let norms: f64 = j.iter().map(|r| (r[0] * r[0] + r[1] * r[1] + r[2] * r[2]).sqrt()).product(); + (det > 1e-12 * norms).then_some((inv, det)) +} + pub fn det3(m: &[[f64; 3]; 3]) -> f64 { m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) diff --git a/crates/fem/src/dof.rs b/crates/fem/src/dof.rs new file mode 100644 index 0000000..6fc7aba --- /dev/null +++ b/crates/fem/src/dof.rs @@ -0,0 +1,458 @@ +//! Degrees of freedom: what a node carries and how the global unknowns are numbered. +//! +//! Every unknown of a discrete system is a ([`DofKind`], node) pair. An element declares the +//! kinds each of its nodes uses ([`ElementDofs`]): three translations for solids, six +//! (translations and rotations) for shells and beams, two for plane elements, a temperature for +//! heat transfer, a pressure for mixed u–p elements, and the internal unknowns it condenses out +//! itself (incompatible modes, a discontinuous pressure), which never reach the global system. +//! A node carries the union of what its elements, constraints and boundary conditions need. +//! +//! [`DofLayout`] numbers them: mesh nodes in label order, each node's kinds in [`DofKind`] order, +//! then extra points (reference points) in the order they are added, so the same model always +//! gives the same numbering (the convention of T. J. R. Hughes, *The Finite Element Method*, +//! Dover 2000, §2.8–2.10, with the ID array replaced by per-node offsets). Assembly reads a +//! node's unknowns as one contiguous block, so a node pair is a dense block of the global matrix. +//! +//! Boundary conditions and loads address a node and a [`DofKind`]; the model's 1-based degree +//! of freedom numbers (1–3 translations, 4–6 rotations, 11 temperature, 8 pressure, the common +//! convention of the deck formats) convert with [`DofKind::from_number`]. + +use std::collections::BTreeMap; +use std::ops::Range; + +/// One kind of nodal unknown. The order is the numbering order within a node. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub enum DofKind { + /// Translation along X. + Ux, + /// Translation along Y. + Uy, + /// Translation along Z. + Uz, + /// Rotation about X. + Rx, + /// Rotation about Y. + Ry, + /// Rotation about Z. + Rz, + /// Temperature. + T, + /// Pressure (hydrostatic, of a mixed u–p element). + P, +} + +impl DofKind { + pub const ALL: [DofKind; 8] = [DofKind::Ux, DofKind::Uy, DofKind::Uz, DofKind::Rx, DofKind::Ry, DofKind::Rz, DofKind::T, DofKind::P]; + + /// Position in [`DofKind::ALL`]. + pub const fn index(self) -> usize { + self as usize + } + + /// The translation along axis `axis` (0–2). + pub fn translation(axis: usize) -> Option { + [DofKind::Ux, DofKind::Uy, DofKind::Uz].get(axis).copied() + } + + /// The rotation about axis `axis` (0–2). + pub fn rotation(axis: usize) -> Option { + [DofKind::Rx, DofKind::Ry, DofKind::Rz].get(axis).copied() + } + + /// The structural direction 0–5 (translations, then rotations), `None` for T and P. + pub fn direction(self) -> Option { + (self.index() < 6).then_some(self.index()) + } + + /// The kind of a structural direction 0–5. + pub fn from_direction(direction: usize) -> Option { + DofKind::ALL.get(direction).copied().filter(|k| k.index() < 6) + } + + /// The model's and decks' 1-based number: 1–3 translations, 4–6 rotations, 8 pressure, + /// 11 temperature. + pub fn number(self) -> u8 { + match self { + DofKind::T => 11, + DofKind::P => 8, + k => k.index() as u8 + 1, + } + } + + /// The kind of a model degree of freedom number ([`DofKind::number`]). + pub fn from_number(number: u8) -> Option { + DofKind::ALL.into_iter().find(|k| k.number() == number) + } + + /// What it is, in the user's terms. + pub fn describe(self) -> &'static str { + match self { + DofKind::Ux => "X", + DofKind::Uy => "Y", + DofKind::Uz => "Z", + DofKind::Rx => "rotation about X", + DofKind::Ry => "rotation about Y", + DofKind::Rz => "rotation about Z", + DofKind::T => "temperature", + DofKind::P => "pressure", + } + } + + const fn bit(self) -> u8 { + 1 << self as u8 + } +} + +/// A set of [`DofKind`]s, iterated in numbering order. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq, Hash)] +pub struct DofSet(u8); + +impl DofSet { + pub const EMPTY: DofSet = DofSet(0); + /// Ux, Uy, Uz: solids, trusses, membranes in 3D. + pub const TRANSLATIONS: DofSet = DofSet(0b0000_0111); + /// Ux, Uy: plane elements. + pub const PLANE: DofSet = DofSet(0b0000_0011); + /// Rx, Ry, Rz. + pub const ROTATIONS: DofSet = DofSet(0b0011_1000); + /// Translations and rotations: shells, beams, reference points. + pub const STRUCTURAL: DofSet = DofSet(0b0011_1111); + pub const TEMPERATURE: DofSet = DofSet(DofKind::T.bit()); + pub const PRESSURE: DofSet = DofSet(DofKind::P.bit()); + + pub const fn of(kinds: &[DofKind]) -> DofSet { + let mut bits = 0; + let mut i = 0; + while i < kinds.len() { + bits |= kinds[i].bit(); + i += 1; + } + DofSet(bits) + } + + pub const fn contains(self, kind: DofKind) -> bool { + self.0 & kind.bit() != 0 + } + + /// Whether every kind of `other` is in this set. + pub const fn includes(self, other: DofSet) -> bool { + self.0 & other.0 == other.0 + } + + pub const fn union(self, other: DofSet) -> DofSet { + DofSet(self.0 | other.0) + } + + pub const fn len(self) -> usize { + self.0.count_ones() as usize + } + + pub const fn is_empty(self) -> bool { + self.0 == 0 + } + + /// The position of `kind` among this set's kinds (its offset in a node's block). + pub const fn position(self, kind: DofKind) -> Option { + if self.contains(kind) { Some((self.0 & (kind.bit() - 1)).count_ones() as usize) } else { None } + } + + /// The kinds, in numbering order. + pub fn iter(self) -> impl Iterator { + DofKind::ALL.into_iter().filter(move |k| self.contains(*k)) + } +} + +impl FromIterator for DofSet { + fn from_iter>(iter: I) -> DofSet { + DofSet(iter.into_iter().fold(0, |m, k| m | k.bit())) + } +} + +/// What an element type contributes to the unknowns: the kinds every one of its nodes carries +/// (its element matrices are ordered node by node, these kinds in order within a node) and the +/// internal unknowns it condenses out before assembly. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct ElementDofs { + pub node: DofSet, + pub internal: Internal, +} + +/// Element-internal unknowns, condensed at the element level (static condensation; Hughes +/// 2000, §4.7 for incompatible modes, §4.4 for discontinuous pressure). +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)] +pub struct Internal { + /// Incompatible (bubble) displacement modes. + pub incompatible_modes: usize, + /// Discontinuous pressure unknowns of a mixed element. + pub pressure: usize, +} + +impl Internal { + pub const NONE: Internal = Internal { incompatible_modes: 0, pressure: 0 }; + + pub const fn count(self) -> usize { + self.incompatible_modes + self.pressure + } +} + +impl ElementDofs { + /// Three translations per node, nothing internal. + pub const SOLID: ElementDofs = ElementDofs { node: DofSet::TRANSLATIONS, internal: Internal::NONE }; + + pub const fn new(node: DofSet) -> ElementDofs { + ElementDofs { node, internal: Internal::NONE } + } + + /// Unknowns of an element with `nodes` nodes, before condensation. + pub const fn total(self, nodes: usize) -> usize { + nodes * self.node.len() + self.internal.count() + } +} + +/// One numbered node: its label, its kinds and its first global unknown. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +struct Slot { + label: u32, + kinds: DofSet, + first: usize, +} + +/// The global numbering of the unknowns. +/// +/// Slots (positions) `0..nodes()` are the mesh nodes in label order; `nodes()..slots()` are the +/// extra points in the order they were added. A slot's unknowns are contiguous, in +/// [`DofKind`] order. +#[derive(Debug, Clone, Default, PartialEq, Eq)] +pub struct DofLayout { + slots: Vec, + nodes: usize, + total: usize, +} + +/// Collects what every node needs before numbering: the union of the kinds asked for. +#[derive(Debug, Clone, Default)] +pub struct LayoutBuilder { + nodes: BTreeMap, + points: Vec<(u32, DofSet)>, +} + +impl LayoutBuilder { + /// Node `label` needs `kinds` (added to what it needs already). + pub fn node(&mut self, label: u32, kinds: DofSet) -> &mut Self { + let entry = self.nodes.entry(label).or_default(); + *entry = entry.union(kinds); + self + } + + /// Every node of an element of this kind needs the element's kinds. + pub fn element(&mut self, labels: &[u32], dofs: ElementDofs) -> &mut Self { + for &l in labels { + self.node(l, dofs.node); + } + self + } + + /// An extra point (a reference point) numbered after every mesh node, in the order added; + /// each call adds one. + pub fn point(&mut self, label: u32, kinds: DofSet) -> &mut Self { + self.points.push((label, kinds)); + self + } + + pub fn build(&self) -> DofLayout { + let mut slots = Vec::with_capacity(self.nodes.len() + self.points.len()); + let mut total = 0; + for (&label, &kinds) in self.nodes.iter().chain(self.points.iter().map(|p| (&p.0, &p.1))) { + slots.push(Slot { label, kinds, first: total }); + total += kinds.len(); + } + DofLayout { slots, nodes: self.nodes.len(), total } + } +} + +impl DofLayout { + pub fn builder() -> LayoutBuilder { + LayoutBuilder::default() + } + + /// Every label in `labels` (sorted or not) with the same kinds. + pub fn uniform(labels: impl IntoIterator, kinds: DofSet) -> DofLayout { + let mut b = LayoutBuilder::default(); + for l in labels { + b.node(l, kinds); + } + b.build() + } + + /// Total unknowns. + pub fn dofs(&self) -> usize { + self.total + } + + /// Mesh nodes (slots `0..nodes()`). + pub fn nodes(&self) -> usize { + self.nodes + } + + /// Mesh nodes and extra points. + pub fn slots(&self) -> usize { + self.slots.len() + } + + /// Unknowns of the mesh nodes; the extra points' follow them. + pub fn node_dofs(&self) -> usize { + self.slots.get(self.nodes).map_or(self.total, |s| s.first) + } + + /// The slot of extra point `k`. + pub fn point_slot(&self, k: usize) -> usize { + self.nodes + k + } + + /// The slot of a mesh node label. + pub fn node(&self, label: u32) -> Option { + self.slots.get(..self.nodes)?.binary_search_by_key(&label, |s| s.label).ok() + } + + /// The slot of a label: a mesh node, else an extra point. + pub fn slot(&self, label: u32) -> Option { + self.node(label).or_else(|| self.slots.get(self.nodes..)?.iter().position(|s| s.label == label).map(|k| self.nodes + k)) + } + + pub fn label(&self, slot: usize) -> Option { + self.slots.get(slot).map(|s| s.label) + } + + /// The kinds a slot carries. + pub fn kinds(&self, slot: usize) -> DofSet { + self.slots.get(slot).map_or(DofSet::EMPTY, |s| s.kinds) + } + + /// A slot's unknowns (empty for a slot that does not exist). + pub fn range(&self, slot: usize) -> Range { + self.slots.get(slot).map_or(0..0, |s| s.first..s.first + s.kinds.len()) + } + + /// The global unknown of `kind` at `slot`, if the slot carries it. + pub fn dof(&self, slot: usize, kind: DofKind) -> Option { + let s = self.slots.get(slot)?; + s.kinds.position(kind).map(|p| s.first + p) + } + + /// The slot and kind of a global unknown. + pub fn locate(&self, dof: usize) -> Option<(usize, DofKind)> { + if dof >= self.total { + return None; + } + let slot = self.slots.partition_point(|s| s.first <= dof).checked_sub(1)?; + let s = self.slots.get(slot)?; + let kind = s.kinds.iter().nth(dof - s.first)?; + Some((slot, kind)) + } + + /// The global unknowns of an element: its slots' `kinds`, node by node, `None` when a slot + /// lacks one of them. + pub fn element_dofs(&self, slots: &[usize], kinds: DofSet) -> Option> { + let mut out = Vec::with_capacity(slots.len() * kinds.len()); + for &slot in slots { + for k in kinds.iter() { + out.push(self.dof(slot, k)?); + } + } + Some(out) + } + + /// `Some(s)` when every slot carries the same `s` kinds, so unknown `s · slot + offset` is + /// the `offset`-th kind of `slot` (node-major vectors of `s`-vectors). + pub fn stride(&self) -> Option { + let first = self.slots.first()?.kinds; + self.slots.iter().all(|s| s.kinds == first).then_some(first.len()) + } + + /// What a global unknown is, in the user's terms. + pub fn describe(&self, dof: usize) -> String { + match self.locate(dof) { + Some((slot, kind)) => { + let what = if slot < self.nodes { "node" } else { "point" }; + format!("{what} {} in {}", self.label(slot).unwrap_or(0), kind.describe()) + } + None => format!("unknown degree of freedom {dof}"), + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn kinds_numbers_and_sets() { + for k in DofKind::ALL { + assert_eq!(DofKind::from_number(k.number()), Some(k)); + } + assert_eq!(DofKind::from_number(1), Some(DofKind::Ux)); + assert_eq!(DofKind::from_number(6), Some(DofKind::Rz)); + assert_eq!(DofKind::from_number(11), Some(DofKind::T)); + assert_eq!(DofKind::from_number(0), None); + assert_eq!(DofKind::from_number(7), None); + assert_eq!(DofKind::from_direction(4), Some(DofKind::Ry)); + assert_eq!(DofKind::from_direction(6), None); + let s = DofSet::of(&[DofKind::T, DofKind::Ux, DofKind::Uy, DofKind::Uz]); + assert_eq!(s.len(), 4); + assert_eq!(s.iter().collect::>(), vec![DofKind::Ux, DofKind::Uy, DofKind::Uz, DofKind::T]); + assert_eq!(s.position(DofKind::T), Some(3)); + assert_eq!(s.position(DofKind::Rx), None); + assert!(DofSet::STRUCTURAL.includes(DofSet::TRANSLATIONS)); + assert!(!DofSet::TRANSLATIONS.includes(DofSet::STRUCTURAL)); + assert_eq!(DofSet::TRANSLATIONS.union(DofSet::ROTATIONS), DofSet::STRUCTURAL); + assert_eq!(ElementDofs::SOLID.total(8), 24); + let hybrid = ElementDofs { node: DofSet::TRANSLATIONS, internal: Internal { incompatible_modes: 0, pressure: 1 } }; + assert_eq!(hybrid.total(8), 25); + } + + #[test] + fn layout_numbers_nodes_by_label_then_kind_then_points() { + let mut b = DofLayout::builder(); + // A shell on 7 and 3, a solid on 3, 5 and 7, a heat element on 5; a reference point. + b.element(&[7, 3], ElementDofs::new(DofSet::STRUCTURAL)); + b.element(&[3, 5, 7], ElementDofs::SOLID); + b.node(5, DofSet::TEMPERATURE); + b.point(100, DofSet::STRUCTURAL); + let l = b.build(); + assert_eq!((l.nodes(), l.slots(), l.dofs(), l.node_dofs()), (3, 4, 22, 16)); + assert_eq!(l.node(3), Some(0)); + assert_eq!(l.node(5), Some(1)); + assert_eq!(l.node(100), None); + assert_eq!(l.slot(100), Some(3)); + assert_eq!(l.range(1), 6..10); + assert_eq!(l.dof(1, DofKind::T), Some(9)); + assert_eq!(l.dof(1, DofKind::Rx), None); + assert_eq!(l.dof(2, DofKind::Rz), Some(15)); + assert_eq!(l.dof(3, DofKind::Ux), Some(16)); + assert_eq!(l.locate(9), Some((1, DofKind::T))); + assert_eq!(l.locate(16), Some((3, DofKind::Ux))); + assert_eq!(l.locate(22), None); + for d in 0..l.dofs() { + let (slot, kind) = l.locate(d).unwrap(); + assert_eq!(l.dof(slot, kind), Some(d)); + } + assert_eq!(l.element_dofs(&[0, 1], DofSet::TRANSLATIONS), Some(vec![0, 1, 2, 6, 7, 8])); + assert_eq!(l.element_dofs(&[0, 1], DofSet::STRUCTURAL), None); + assert_eq!(l.describe(9), "node 5 in temperature"); + assert_eq!(l.describe(19), "point 100 in rotation about X"); + // The same needs in another order give the same layout. + let mut c = DofLayout::builder(); + c.point(100, DofSet::STRUCTURAL).node(5, DofSet::TEMPERATURE).element(&[5, 7, 3], ElementDofs::SOLID); + c.element(&[3, 7], ElementDofs::new(DofSet::STRUCTURAL)); + assert_eq!(c.build(), l); + // A uniform translational layout is the classic 3 × node numbering. + let u = DofLayout::uniform([4, 2, 9], DofSet::TRANSLATIONS); + assert_eq!(u.dof(2, DofKind::Uy), Some(7)); + assert_eq!(u.dofs(), 9); + assert_eq!(u.range(7), 0..0); + assert_eq!(u.locate(usize::MAX), None); + assert_eq!(u.stride(), Some(3)); + assert_eq!(l.stride(), None); + assert_eq!(DofLayout::default().stride(), None); + } +} diff --git a/crates/fem/src/element.rs b/crates/fem/src/element.rs index 728554a..2ef905a 100644 --- a/crates/fem/src/element.rs +++ b/crates/fem/src/element.rs @@ -34,7 +34,8 @@ use modelcraft_model::ElementKind; -use crate::dense::{Cholesky, Matrix, inverse3}; +use crate::dense::{Cholesky, Matrix, inverse3, positive_inverse}; +use crate::dof::{DofSet, ElementDofs, Internal}; use crate::material::{Elasticity, stress}; use crate::quadrature::{Point, Rule}; use crate::shape::{self, MAX_NODES, Shape}; @@ -144,6 +145,13 @@ pub enum Formulation { } impl Formulation { + /// The unknowns of a solid element so formulated: three translations per node, and the + /// incompatible modes it condenses out. + pub fn dofs(self) -> ElementDofs { + let incompatible_modes = if self == Formulation::IncompatibleModes { INCOMPATIBLE_DOFS } else { 0 }; + ElementDofs { node: DofSet::TRANSLATIONS, internal: Internal { incompatible_modes, pressure: 0 } } + } + /// The formulation a model section asks for, if this element type has it. pub fn from_section(kind: ElementKind, f: &modelcraft_model::Formulation) -> Result { let ty = element_type(kind)?; @@ -220,6 +228,11 @@ impl<'a> Solid<'a> { Ok(Solid { ty, xyz, formulation }) } + /// The element's unknowns, in the order of its matrices. + pub fn dofs(&self) -> ElementDofs { + self.formulation.dofs() + } + /// The stiffness integration rule. pub fn rule(&self) -> Rule { match (self.formulation, self.ty.reduced) { @@ -249,11 +262,7 @@ impl<'a> Solid<'a> { fn gradients(&self, xi: [f64; 3]) -> Result { let shape = (self.ty.shape)(xi); let j = self.jacobian(&shape); - let (inv, det) = inverse3(&j).ok_or_else(|| self.degenerate())?; - let norms: f64 = j.iter().map(|r| (r[0] * r[0] + r[1] * r[1] + r[2] * r[2]).sqrt()).product(); - if !(det > 1e-12 * norms) { - return Err(self.degenerate()); - } + let (inv, det) = positive_inverse(&j).ok_or_else(|| self.degenerate())?; let mut dx = [[0.0; 3]; MAX_NODES]; for (g, d) in dx.iter_mut().zip(shape.derivatives()) { for b in 0..3 { diff --git a/crates/fem/src/jet.rs b/crates/fem/src/jet.rs new file mode 100644 index 0000000..1c0f3fe --- /dev/null +++ b/crates/fem/src/jet.rs @@ -0,0 +1,183 @@ +//! Second-order forward-mode automatic differentiation in three variables. +//! +//! A [`Jet`] carries a value with its gradient and Hessian with respect to the natural +//! coordinates (ξ, η, ζ). Arithmetic propagates them by the product and quotient rules, so a +//! shape function written once as a formula yields its exact first and second derivatives (to +//! round-off, no finite differences): the hyper-dual / truncated Taylor arithmetic of A. Griewank, +//! A. Walther, *Evaluating Derivatives*, 2nd ed., SIAM 2008, ch. 13. + +use std::ops::{Add, Div, Mul, Neg, Sub}; + +/// A value, its gradient and its (symmetric) Hessian. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Jet { + pub v: f64, + pub g: [f64; 3], + pub h: [[f64; 3]; 3], +} + +impl Jet { + pub const fn constant(v: f64) -> Jet { + Jet { v, g: [0.0; 3], h: [[0.0; 3]; 3] } + } + + /// The independent variable `k` (0–2) at value `v`. + pub fn var(v: f64, k: usize) -> Jet { + let mut j = Jet::constant(v); + if let Some(g) = j.g.get_mut(k) { + *g = 1.0; + } + j + } + + /// The three natural coordinates as variables. + pub fn point(xi: [f64; 3]) -> [Jet; 3] { + [Jet::var(xi[0], 0), Jet::var(xi[1], 1), Jet::var(xi[2], 2)] + } + + /// 1 / self; the Hessian term follows from (1/b)'' = 2 b'b'ᵀ / b³ − b'' / b². + pub fn recip(self) -> Jet { + let r = 1.0 / self.v; + let r2 = r * r; + let mut out = Jet::constant(r); + for a in 0..3 { + out.g[a] = -self.g[a] * r2; + for b in 0..3 { + out.h[a][b] = 2.0 * self.g[a] * self.g[b] * r2 * r - self.h[a][b] * r2; + } + } + out + } + + fn map(self, other: Jet, f: impl Fn(f64, f64) -> f64) -> Jet { + let mut out = Jet::constant(f(self.v, other.v)); + for a in 0..3 { + out.g[a] = f(self.g[a], other.g[a]); + for b in 0..3 { + out.h[a][b] = f(self.h[a][b], other.h[a][b]); + } + } + out + } + + fn scale(self, s: f64) -> Jet { + Jet { v: self.v * s, g: self.g.map(|x| x * s), h: self.h.map(|r| r.map(|x| x * s)) } + } +} + +impl Add for Jet { + type Output = Jet; + fn add(self, o: Jet) -> Jet { + self.map(o, |a, b| a + b) + } +} + +impl Sub for Jet { + type Output = Jet; + fn sub(self, o: Jet) -> Jet { + self.map(o, |a, b| a - b) + } +} + +impl Neg for Jet { + type Output = Jet; + fn neg(self) -> Jet { + self.scale(-1.0) + } +} + +impl Mul for Jet { + type Output = Jet; + fn mul(self, o: Jet) -> Jet { + let mut out = Jet::constant(self.v * o.v); + for a in 0..3 { + out.g[a] = self.v * o.g[a] + o.v * self.g[a]; + for b in 0..3 { + out.h[a][b] = self.v * o.h[a][b] + o.v * self.h[a][b] + self.g[a] * o.g[b] + o.g[a] * self.g[b]; + } + } + out + } +} + +impl Div for Jet { + type Output = Jet; + #[allow(clippy::suspicious_arithmetic_impl)] + fn div(self, o: Jet) -> Jet { + self * o.recip() + } +} + +impl Add for Jet { + type Output = Jet; + fn add(mut self, c: f64) -> Jet { + self.v += c; + self + } +} + +impl Sub for Jet { + type Output = Jet; + fn sub(mut self, c: f64) -> Jet { + self.v -= c; + self + } +} + +impl Mul for Jet { + type Output = Jet; + fn mul(self, c: f64) -> Jet { + self.scale(c) + } +} + +impl Add for f64 { + type Output = Jet; + fn add(self, j: Jet) -> Jet { + j + self + } +} + +impl Sub for f64 { + type Output = Jet; + fn sub(self, j: Jet) -> Jet { + -j + self + } +} + +impl Mul for f64 { + type Output = Jet; + fn mul(self, j: Jet) -> Jet { + j.scale(self) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn products_and_quotients_differentiate_exactly() { + let [x, y, z] = Jet::point([0.3, -0.7, 1.9]); + // f = x² y / (1 + z) + 3 x z. + let f = x * x * y / (1.0 + z) + 3.0 * x * z; + let (xv, yv, zv) = (0.3, -0.7, 1.9); + let d = 1.0 + zv; + assert!((f.v - (xv * xv * yv / d + 3.0 * xv * zv)).abs() < 1e-15); + let g = [2.0 * xv * yv / d + 3.0 * zv, xv * xv / d, -xv * xv * yv / (d * d) + 3.0 * xv]; + let h = [ + [2.0 * yv / d, 2.0 * xv / d, -2.0 * xv * yv / (d * d) + 3.0], + [2.0 * xv / d, 0.0, -xv * xv / (d * d)], + [-2.0 * xv * yv / (d * d) + 3.0, -xv * xv / (d * d), 2.0 * xv * xv * yv / (d * d * d)], + ]; + for a in 0..3 { + assert!((f.g[a] - g[a]).abs() < 1e-14, "g{a}"); + for b in 0..3 { + assert!((f.h[a][b] - h[a][b]).abs() < 1e-14, "h{a}{b}: {} vs {}", f.h[a][b], h[a][b]); + } + } + let s = 2.0 - x * 0.5 - (y - 1.0); + assert_eq!(s.g, [-0.5, -1.0, 0.0]); + assert_eq!(s.v, 2.0 - 0.15 + 1.7); + } +} diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index 99fbece..53e0ec8 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -3,8 +3,15 @@ //! Pure functions over element geometry and material constants, with no solver loop and no //! global state, so the implicit solver (`modelcraft-solver`) and the explicit solver that follows //! it use the same element technology: +//! - [`dof`]: degree-of-freedom kinds, what elements declare, and the global numbering +//! ([`DofLayout`]) every solver assembles on; //! - [`quadrature`]: Gauss–Legendre and tetrahedron rules; -//! - [`shape`]: shape functions of tet4, tet10, wedge6, hex8, hex20; +//! - [`basis`]: the shape-function library of every family (lines, triangles, quadrilaterals, +//! tetrahedra, wedges, pyramids, hexahedra) with first and second derivatives and their +//! physical gradients, Hessians and Laplacians; [`jet`] is the automatic differentiation behind +//! the second derivatives; +//! - [`shape`]: the hand-written first-order functions of tet4, tet10, wedge6, hex8, hex20 that +//! the solid kernels evaluate; //! - [`material`]: isotropic and orthotropic elasticity matrices, von Mises stress; //! - [`element`]: the element table, stiffness (full, reduced, incompatible modes), consistent //! and lumped mass, body forces, integration-point stresses and their extrapolation to nodes; @@ -13,21 +20,27 @@ //! shape vectors and element frequency bounds for explicit dynamics; //! - [`dense`]: the small dense linear algebra inside an element. //! -//! Conventions: three translational degrees of freedom per node, ordered node by node; stress +//! 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 //! and strain 6-vectors in the order xx, yy, zz, xy, xz, yz with engineering shear strains. #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] #![forbid(unsafe_code)] // Matrix kernels over fixed small dimensions read as the textbook formulas with index loops. #![allow(clippy::needless_range_loop)] +pub mod basis; pub mod dense; +pub mod dof; pub mod element; pub mod face; +pub mod jet; pub mod material; pub mod quadrature; pub mod shape; pub mod uniform; +pub use basis::Basis; +pub use dof::{DofKind, DofLayout, DofSet, ElementDofs}; pub use element::{ELEMENT_TYPES, ElementType, Family, Formulation, PointStress, Solid, element_type}; pub use material::{Elasticity, elasticity, principal, von_mises}; diff --git a/crates/fem/src/shape.rs b/crates/fem/src/shape.rs index 76da784..76f2de1 100644 --- a/crates/fem/src/shape.rs +++ b/crates/fem/src/shape.rs @@ -21,8 +21,8 @@ use modelcraft_model::ElementKind; -/// Most nodes an element in this crate has (hex20). -pub const MAX_NODES: usize = 20; +/// Most nodes an element family has (hex27, [`crate::basis::Basis::Hex27`]). +pub const MAX_NODES: usize = 27; /// Shape function values and natural derivatives at one point. #[derive(Debug, Clone, Copy)] @@ -157,33 +157,12 @@ pub fn hex20(xi: [f64; 3]) -> Shape { mod tests { use super::*; - fn nodes_of(kind: &str) -> Vec<[f64; 3]> { - match kind { - "tet4" => TET_CORNERS.to_vec(), - "tet10" => { - let mut v = TET_CORNERS.to_vec(); - for e in ElementKind::Tet10.edges() { - v.push([0, 1, 2].map(|k| (TET_CORNERS[e[0]][k] + TET_CORNERS[e[1]][k]) / 2.0)); - } - v - } - "hex8" => HEX_CORNERS.to_vec(), - "wedge6" => WEDGE_CORNERS.to_vec(), - _ => { - let mut v = HEX_CORNERS.to_vec(); - for e in ElementKind::Hex20.edges() { - v.push([0, 1, 2].map(|k| (HEX_CORNERS[e[0]][k] + HEX_CORNERS[e[1]][k]) / 2.0)); - } - v - } - } - } - #[test] fn kronecker_delta_partition_of_unity_and_derivatives() { let kinds: [(&str, fn([f64; 3]) -> Shape); 5] = [("tet4", tet4), ("tet10", tet10), ("wedge6", wedge6), ("hex8", hex8), ("hex20", hex20)]; for (name, f) in kinds { - let nodes = nodes_of(name); + let nodes = crate::Basis::ALL.iter().find(|b| b.name() == name).map(|b| b.nodes()).unwrap_or_default(); + assert!(!nodes.is_empty()); for (i, p) in nodes.iter().enumerate() { let s = f(*p); for (j, v) in s.values().iter().enumerate() { @@ -197,7 +176,7 @@ mod tests { assert!((s.values().iter().sum::() - 1.0).abs() < 1e-14); for k in 0..3 { assert!(s.derivatives().iter().map(|d| d[k]).sum::().abs() < 1e-13); - let x: f64 = s.values().iter().zip(&nodes).map(|(n, x)| n * x[k]).sum(); + let x: f64 = s.values().iter().zip(nodes).map(|(n, x)| n * x[k]).sum(); assert!((x - p[k]).abs() < 1e-14, "{name}"); let h = 1e-6; let (mut a, mut b) = (p, p); diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 01a619b..a143e59 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -2,7 +2,7 @@ use modelcraft_fem::element::{Extrapolation, extrapolation}; use modelcraft_fem::quadrature::Rule; -use modelcraft_fem::{PointStress, von_mises}; +use modelcraft_fem::{DofKind, DofLayout, PointStress, von_mises}; use modelcraft_model::{ElementKind, Model, Procedure, Resolved, Step}; use serde::Serialize; @@ -17,7 +17,8 @@ use crate::{Error, Limits, Result}; /// Size and effort of a solution. #[derive(Debug, Clone, Default, PartialEq, Serialize)] pub struct Stats { - /// Degrees of freedom: 3 per mesh node and 6 per coupled reference point. + /// Degrees of freedom: those of the mesh nodes (3 per solid node) and 6 per coupled + /// reference point. pub dofs: usize, /// Free degrees of freedom solved for. pub equations: usize, @@ -129,17 +130,18 @@ struct Assembled { /// Assemble the stiffness (and the mass when `mass` is set) on the mesh's node graph. fn assemble(mesh: &Mesh, graph: &NodeGraph, mass: Option) -> Result { - let mut k = graph.zero_matrix(); - let mut m = mass.map(|_| graph.zero_matrix()); + let layout = &mesh.layout; + let mut k = graph.zero_matrix(layout); + let mut m = mass.map(|_| graph.zero_matrix(layout)); for e in &mesh.elements { let solid = e.solid()?; let at = |err| Error::element(e.label, err); - add_block(&mut k, &e.nodes, &solid.stiffness(&e.d).map_err(at)?.data)?; + add_block(&mut k, layout, &e.nodes, solid.dofs().node, &solid.stiffness(&e.d).map_err(at)?.data)?; if let (Some(mm), Some(kind)) = (m.as_mut(), mass) { let rho = e.density.ok_or_else(|| Error::Invalid(format!("element {} has no density", e.label)))?; match kind { - MassKind::Consistent => add_scalar(mm, &e.nodes, &solid.mass(rho).map_err(at)?.data)?, - MassKind::Lumped => add_diagonal(mm, &e.nodes, &solid.lumped_mass(rho).map_err(at)?)?, + MassKind::Consistent => add_scalar(mm, layout, &e.nodes, solid.dofs().node, &solid.mass(rho).map_err(at)?.data)?, + MassKind::Lumped => add_diagonal(mm, layout, &e.nodes, solid.dofs().node, &solid.lumped_mass(rho).map_err(at)?)?, } } } @@ -167,10 +169,17 @@ fn factor_error(e: FactorError, mesh: &Mesh, reduction: &Reduction) -> Error { } } -/// The mesh nodes' part of a full vector, per node, and the reference points' part, per point. +/// The structural directions (0..N) of layout slot `slot` in a full vector, zero where the slot +/// has none. +fn directions(full: &[f64], layout: &DofLayout, slot: usize) -> [f64; N] { + 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]>) { - let (nodes, points) = full.split_at(mesh.mesh_dofs().min(full.len())); - (nodes.as_chunks::<3>().0.to_vec(), points.as_chunks::<6>().0.to_vec()) + 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()) } /// Constraint rows, prescribed values and their elimination for a step. @@ -188,8 +197,8 @@ fn reduction(s: &Setup, values: bool, order: &[usize], budget: Budget, warnings: } /// The reduced matrices' memory budget: the factorization's. -fn budget(limits: &Limits, graph: &NodeGraph, matrices: usize) -> Budget { - Budget { max_bytes: limits.max_factor_entries.saturating_mul(8), pattern_entries: graph.entries(), matrices } +fn budget(limits: &Limits, mesh: &Mesh, graph: &NodeGraph, matrices: usize) -> Budget { + Budget { max_bytes: limits.max_factor_entries.saturating_mul(8), pattern_entries: graph.entries(&mesh.layout), matrices } } /// Solve a linear static step: the step named `step`, or the model's only static step. @@ -201,7 +210,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); let order = nested_dissection(&mesh.xyz, &graph.adjacency); let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; - let red = reduction(&setup, true, &order, budget(limits, &graph, 1), &mut warnings)?; + let red = reduction(&setup, true, &order, budget(limits, &mesh, &graph, 1), &mut warnings)?; let f = loads(model, &resolved, &mesh, step, &red.active)?; let k = assemble(&mesh, &graph, None)?.k; if red.neq() == 0 { @@ -221,7 +230,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul let (reaction, 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 d = mesh.direction(dof); + let Some(d) = mesh.kind(dof).and_then(DofKind::direction) else { continue }; if let (Some(a), Some(r)) = (applied_total.get_mut(d), reaction_total.get_mut(d)) { *a += fv; *r += rv; @@ -328,7 +337,7 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt let graph = NodeGraph::new(mesh.labels.len(), mesh.elements.iter().map(|e| e.nodes.as_slice())); let order = nested_dissection(&mesh.xyz, &graph.adjacency); let setup = Setup { model, resolved: &resolved, mesh: &mesh, step }; - let red = reduction(&setup, false, &order, budget(limits, &graph, 2), &mut warnings)?; + 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. @@ -372,9 +381,10 @@ pub fn solve_frequency(model: &Model, step: Option<&str>, options: &FrequencyOpt let pairs = lanczos(&op, shift, wanted, BLOCK, max_basis, TOLERANCE)?; // Participation: Γ_d = φᵀ M r_d with r_d the unit translation along d on the free equations // (mesh and reference point translations). + let directions: Vec> = red.eq_dof.iter().map(|dof| mesh.kind(*dof).and_then(DofKind::direction)).collect(); let mr: Vec> = (0..3) .map(|d| { - let r: Vec = red.eq_dof.iter().map(|dof| if mesh.direction(*dof) == d { 1.0 } else { 0.0 }).collect(); + let r: Vec = directions.iter().map(|dir| if *dir == Some(d) { 1.0 } else { 0.0 }).collect(); let mut y = vec![0.0; neq]; mff.mul(&r, &mut y); y diff --git a/crates/solver/src/constraints.rs b/crates/solver/src/constraints.rs index a548347..fd5ee5a 100644 --- a/crates/solver/src/constraints.rs +++ b/crates/solver/src/constraints.rs @@ -50,6 +50,7 @@ use std::collections::BTreeMap; +use modelcraft_fem::DofKind; use modelcraft_fem::dense::inverse3; use modelcraft_model::{CouplingKind, Model, Resolved}; @@ -94,7 +95,10 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> for (name, c) in model.couplings() { let point = mesh.point(&c.reference_point).ok_or_else(|| Error::Internal(format!("reference point `{}`", c.reference_point)))?; let xr = mesh.points.get(point).map(|p| p.xyz).ok_or_else(|| Error::Internal("reference point index".into()))?; - let rp = mesh.point_dof(point, 0); + let slot = mesh.layout.point_slot(point); + // Reference point direction c (0–5), and translation k of mesh node n. + let rp = |c: usize| DofKind::from_direction(c).ok_or_else(|| Error::Internal("direction".into())).and_then(|kind| mesh.dof(slot, kind)); + let u = |n: usize, k: usize| DofKind::translation(k).ok_or_else(|| Error::Internal("axis".into())).and_then(|kind| mesh.dof(n, kind)); let set = resolved.sets.get(&c.set).ok_or_else(|| Error::Invalid(format!("coupling `{name}`: set `{}` is not resolved", c.set)))?; if set.nodes.is_empty() { return Err(Error::Invalid(format!("coupling `{name}`: set `{}` has no nodes", c.set))); @@ -119,14 +123,15 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> for k in coupled(0..3) { // u_s,k − u_r,k − (θ × r)_k = 0, (θ × r)_k = θ_a r_b − θ_b r_a. let (a, b) = ((k + 1) % 3, (k + 2) % 3); - let mut terms = vec![(3 * n + k, 1.0), (rp + k, -1.0)]; + let pivot = u(n, k)?; + let mut terms = vec![(pivot, 1.0), (rp(k)?, -1.0)]; if r[b] != 0.0 { - terms.push((rp + 3 + a, -r[b])); + terms.push((rp(3 + a)?, -r[b])); } if r[a] != 0.0 { - terms.push((rp + 3 + b, r[a])); + terms.push((rp(3 + b)?, r[a])); } - rows.push(Row { terms, pivot: 3 * n + k, coupling: name }); + rows.push(Row { terms, pivot, coupling: name }); } } } @@ -162,7 +167,8 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> let g: Vec<[[f64; 3]; 3]> = arms.iter().map(|r| inverse.map_or([[0.0; 3]; 3], |inv| mat_mul(&inv, &cross_matrix(*r)))).collect(); for k in coupled(0..6) { let mut terms = Vec::with_capacity(1 + 3 * nodes.len()); - terms.push((rp + k, 1.0)); + let pivot = rp(k)?; + terms.push((pivot, 1.0)); for (&n, gi) in nodes.iter().zip(&g) { for comp in 0..3 { let a = if k < 3 { @@ -172,11 +178,11 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> gi[k - 3][comp] }; if a != 0.0 { - terms.push((3 * n + comp, -a)); + terms.push((u(n, comp)?, -a)); } } } - rows.push(Row { terms, pivot: rp + k, coupling: name }); + rows.push(Row { terms, pivot, coupling: name }); } } } @@ -184,8 +190,8 @@ pub fn coupling_rows<'m>(model: &'m Model, resolved: &Resolved, mesh: &Mesh) -> Ok(rows) } -/// The degrees of freedom that exist: every mesh translation, and the reference point -/// directions some coupling row involves. +/// The degrees of freedom that exist: every mesh node's, and the reference point directions +/// some coupling row involves. fn active_dofs(mesh: &Mesh, rows: &[Row]) -> Vec { let mut active = vec![false; mesh.dofs()]; active.iter_mut().take(mesh.mesh_dofs()).for_each(|a| *a = true); @@ -288,8 +294,8 @@ fn combine(mut terms: Vec<(usize, f64)>) -> Vec<(usize, f64)> { pub struct Reduction { /// The degree of freedom of each free equation. pub eq_dof: Vec, - /// Whether a degree of freedom exists (mesh translations, and the reference point - /// directions some coupling ties). + /// Whether a degree of freedom exists (the mesh nodes', and the reference point directions + /// some coupling ties). pub active: Vec, /// T_f: per degree of freedom, (equation, coefficient). free: Terms, @@ -376,7 +382,7 @@ impl Reduction { for (name, mut dofs) in coupled { dofs.sort_unstable(); dofs.dedup(); - let mut nodes: Vec = dofs.iter().map(|d| d / 3).collect(); + let mut nodes: Vec = dofs.iter().filter_map(|d| mesh.layout.locate(*d).map(|(slot, _)| slot)).collect(); nodes.dedup(); blocks.push((name, nodes.len(), dofs.len().saturating_mul(dofs.len()))); } @@ -456,7 +462,7 @@ impl Reduction { // Equations: free mesh degrees of freedom in node order, then reference point ones. let mut eq = vec![usize::MAX; n]; let mut eq_dof = Vec::new(); - let mesh_order = order.iter().flat_map(|node| (0..3).map(move |c| 3 * node + c)); + let mesh_order = order.iter().flat_map(|node| mesh.layout.range(*node)); for dof in mesh_order.chain(mesh.mesh_dofs()..n) { if state.get(dof) == Some(&State::Free) && let Some(e) = eq.get_mut(dof) diff --git a/crates/solver/src/sparse.rs b/crates/solver/src/sparse.rs index 9bee509..c078d32 100644 --- a/crates/solver/src/sparse.rs +++ b/crates/solver/src/sparse.rs @@ -1,11 +1,14 @@ //! Sparse symmetric matrices in compressed rows, both triangles stored. //! //! Assembly works on the node graph: a node's row holds every node it shares an element with, -//! and each node pair is a dense 3 × 3 block, so the pattern is built once from the element -//! connectivity and the stiffness and mass share it. Sums run in a fixed order (elements by +//! and each node pair is a dense block of their degrees of freedom ([`DofLayout`]: a node's are +//! contiguous), so the pattern is built once from the element connectivity and the stiffness and +//! mass share it. Sums run in a fixed order (elements by //! label, then the element's own node order), so assembled values are the same bit for bit on //! every run and platform. +use modelcraft_fem::{DofLayout, DofSet}; + use crate::{Error, Result}; /// A square sparse matrix in compressed rows with sorted column indices. @@ -95,21 +98,22 @@ impl NodeGraph { NodeGraph { adjacency } } - /// Nonzeros of the dof-level matrix (3 × 3 per node pair). - pub fn entries(&self) -> usize { - self.adjacency.iter().map(|r| 9 * r.len()).sum() + /// Nonzeros of the dof-level matrix: a dense block per node pair. + pub fn entries(&self, layout: &DofLayout) -> usize { + self.adjacency.iter().enumerate().map(|(a, row)| layout.range(a).len() * row.iter().map(|b| layout.range(*b).len()).sum::()).sum() } - /// A zero matrix over every degree of freedom (3 per node) on this pattern. - pub fn zero_matrix(&self) -> Csr { - let n = 3 * self.adjacency.len(); + /// A zero matrix over the mesh nodes' degrees of freedom (layout slots `0..nodes`) on this + /// pattern. + pub fn zero_matrix(&self, layout: &DofLayout) -> Csr { + let n = layout.node_dofs(); let mut ptr = Vec::with_capacity(n + 1); ptr.push(0); - let mut idx = Vec::with_capacity(self.entries()); - for row in &self.adjacency { - for _ in 0..3 { + let mut idx = Vec::with_capacity(self.entries(layout)); + for (a, row) in self.adjacency.iter().enumerate() { + for _ in layout.range(a) { for &b in row { - idx.extend([3 * b, 3 * b + 1, 3 * b + 2]); + idx.extend(layout.range(b)); } ptr.push(idx.len()); } @@ -119,30 +123,56 @@ impl NodeGraph { } } -/// Where node b's 3 × 3 block starts in each of node a's three rows (the rows of a node share -/// one pattern, built by [`NodeGraph::zero_matrix`]): one search per node pair. -fn block_offsets(k: &Csr, a: usize, b: usize) -> Result<[usize; 3]> { - let (start, end) = (k.ptr.get(3 * a).copied().unwrap_or(0), k.ptr.get(3 * a + 1).copied().unwrap_or(0)); - let cols = k.idx.get(start..end).ok_or_else(|| Error::Internal("pattern".into()))?; - let at = cols.binary_search(&(3 * b)).map_err(|_| Error::Internal("entry outside the pattern".into()))?; - Ok([0, 1, 2].map(|r| k.ptr.get(3 * a + r).copied().unwrap_or(0) + at)) +/// Most kinds a node carries ([`modelcraft_fem::DofKind::ALL`]). +const MAX_KINDS: usize = 8; + +/// Where node b's block lies in node a's rows of a matrix built by [`NodeGraph::zero_matrix`], +/// for the degrees of freedom `kinds` of both: the value index of (row r, column c) is +/// `rows[r] + cols[c]`. The rows of a node share one pattern, so this is one search per node +/// pair. +struct Block { + rows: [usize; MAX_KINDS], + cols: [usize; MAX_KINDS], } -/// Add an element matrix (3n × 3n, row-major, degrees of freedom node by node) into a matrix -/// built by [`NodeGraph::zero_matrix`]; `nodes` are the element's global node indices. -pub fn add_block(k: &mut Csr, nodes: &[usize], data: &[f64]) -> Result<()> { - let width = 3 * nodes.len(); +impl Block { + fn new(k: &Csr, layout: &DofLayout, a: usize, b: usize, kinds: DofSet) -> Result { + let first = layout.range(a).start; + let (start, end) = (k.ptr.get(first).copied().unwrap_or(0), k.ptr.get(first + 1).copied().unwrap_or(0)); + let pattern = k.idx.get(start..end).ok_or_else(|| Error::Internal("pattern".into()))?; + let at = pattern.binary_search(&layout.range(b).start).map_err(|_| Error::Internal("entry outside the pattern".into()))?; + let missing = |kind| Error::Internal(format!("an element degree of freedom ({kind:?}) its node does not have")); + let mut block = Block { rows: [0; MAX_KINDS], cols: [0; MAX_KINDS] }; + for ((row, col), kind) in block.rows.iter_mut().zip(&mut block.cols).zip(kinds.iter()) { + let r = layout.dof(a, kind).ok_or_else(|| missing(kind))?; + *row = k.ptr.get(r).copied().unwrap_or(0) + at; + *col = layout.kinds(b).position(kind).ok_or_else(|| missing(kind))?; + } + Ok(block) + } +} + +fn add_at(k: &mut Csr, at: usize, v: f64) { + if let Some(x) = k.val.get_mut(at) { + *x += v; + } +} + +/// Add an element matrix into a matrix built by [`NodeGraph::zero_matrix`]: row-major over the +/// element's degrees of freedom, node by node and `kinds` in order within a node; `nodes` are +/// the element's layout slots. +pub fn add_block(k: &mut Csr, layout: &DofLayout, nodes: &[usize], kinds: DofSet, data: &[f64]) -> Result<()> { + let per = kinds.len(); + let width = per * nodes.len(); if data.len() != width * width { return Err(Error::Internal("element matrix size".into())); } for (i, &a) in nodes.iter().enumerate() { for (j, &b) in nodes.iter().enumerate() { - let rows = block_offsets(k, a, b)?; - for (r, start) in rows.iter().enumerate() { - for c in 0..3 { - if let (Some(x), Some(v)) = (k.val.get_mut(start + c), data.get((3 * i + r) * width + 3 * j + c)) { - *x += v; - } + let block = Block::new(k, layout, a, b, kinds)?; + for r in 0..per { + for c in 0..per { + add_at(k, block.rows[r] + block.cols[c], data.get((per * i + r) * width + per * j + c).copied().unwrap_or(0.0)); } } } @@ -150,34 +180,31 @@ pub fn add_block(k: &mut Csr, nodes: &[usize], data: &[f64]) -> Result<()> { Ok(()) } -/// Add an n × n element matrix to each displacement component's block (consistent mass). -pub fn add_scalar(k: &mut Csr, nodes: &[usize], data: &[f64]) -> Result<()> { +/// Add an n × n element matrix to the block of each of `kinds` (a solid's consistent mass on its +/// translations). +pub fn add_scalar(k: &mut Csr, layout: &DofLayout, nodes: &[usize], kinds: DofSet, data: &[f64]) -> Result<()> { let n = nodes.len(); if data.len() != n * n { return Err(Error::Internal("element matrix size".into())); } for (i, &a) in nodes.iter().enumerate() { for (j, &b) in nodes.iter().enumerate() { - let rows = block_offsets(k, a, b)?; + let block = Block::new(k, layout, a, b, kinds)?; let v = data.get(i * n + j).copied().unwrap_or(0.0); - for (c, start) in rows.iter().enumerate() { - if let Some(x) = k.val.get_mut(start + c) { - *x += v; - } + for c in 0..kinds.len() { + add_at(k, block.rows[c] + block.cols[c], v); } } } Ok(()) } -/// Add a per-node diagonal to each displacement component (lumped mass). -pub fn add_diagonal(k: &mut Csr, nodes: &[usize], diagonal: &[f64]) -> Result<()> { +/// Add a per-node diagonal to each of `kinds` (a solid's lumped mass on its translations). +pub fn add_diagonal(k: &mut Csr, layout: &DofLayout, nodes: &[usize], kinds: DofSet, diagonal: &[f64]) -> Result<()> { for (&a, v) in nodes.iter().zip(diagonal) { - let rows = block_offsets(k, a, a)?; - for (c, start) in rows.iter().enumerate() { - if let Some(x) = k.val.get_mut(start + c) { - *x += v; - } + let block = Block::new(k, layout, a, a, kinds)?; + for c in 0..kinds.len() { + add_at(k, block.rows[c] + block.cols[c], *v); } } Ok(()) @@ -192,14 +219,15 @@ mod tests { // Two "elements" sharing node 1. let g = NodeGraph::new(3, [vec![0usize, 1], vec![1, 2]].into_iter()); assert_eq!(g.adjacency, vec![vec![0, 1], vec![0, 1, 2], vec![1, 2]]); - let mut k = g.zero_matrix(); + let layout = DofLayout::uniform(0..3, DofSet::TRANSLATIONS); + let mut k = g.zero_matrix(&layout); let mut e = vec![0.0; 36]; for i in 0..6 { e[i * 6 + i] = 2.0; e[i * 6 + (i + 3) % 6] = -1.0; } - add_block(&mut k, &[0, 1], &e).unwrap(); - add_block(&mut k, &[1, 2], &e).unwrap(); + add_block(&mut k, &layout, &[0, 1], DofSet::TRANSLATIONS, &e).unwrap(); + add_block(&mut k, &layout, &[1, 2], DofSet::TRANSLATIONS, &e).unwrap(); assert_eq!(k.diagonal(), vec![2.0, 2.0, 2.0, 4.0, 4.0, 4.0, 2.0, 2.0, 2.0]); let x = vec![1.0; 9]; let mut y = vec![0.0; 9]; @@ -209,14 +237,38 @@ mod tests { let mut mirrored = k.clone(); mirrored.mirror_lower(); assert_eq!(mirrored, k); - assert!(add_block(&mut k, &[0, 2], &e).is_err(), "nodes 0 and 2 share no element"); + assert!(add_block(&mut k, &layout, &[0, 2], DofSet::TRANSLATIONS, &e).is_err(), "nodes 0 and 2 share no element"); // Scalar (per component) and diagonal additions land on each component. - let mut m = g.zero_matrix(); - add_scalar(&mut m, &[0, 1], &[2.0, 1.0, 1.0, 2.0]).unwrap(); - add_diagonal(&mut m, &[2], &[5.0]).unwrap(); + let mut m = g.zero_matrix(&layout); + add_scalar(&mut m, &layout, &[0, 1], DofSet::TRANSLATIONS, &[2.0, 1.0, 1.0, 2.0]).unwrap(); + add_diagonal(&mut m, &layout, &[2], DofSet::TRANSLATIONS, &[5.0]).unwrap(); assert_eq!(m.diagonal(), vec![2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 5.0, 5.0, 5.0]); assert_eq!(m.row(0).1, &[2.0, 0.0, 0.0, 1.0, 0.0, 0.0]); let s = k.add_scaled(-1.0, &k).unwrap(); assert!(s.val.iter().all(|v| *v == 0.0)); } + + #[test] + fn mixed_kinds_per_node_assemble_into_their_blocks() { + // A solid on nodes 0 and 1 (translations) and a beam-like element on 1 and 2 (six each): + // node 1 carries six, and the solid's block lands on its translations. + let mut b = DofLayout::builder(); + b.node(10, DofSet::TRANSLATIONS).node(11, DofSet::STRUCTURAL).node(12, DofSet::STRUCTURAL); + let layout = b.build(); + let g = NodeGraph::new(3, [vec![0usize, 1], vec![1, 2]].into_iter()); + let mut k = g.zero_matrix(&layout); + assert_eq!((k.n, k.idx.len(), g.entries(&layout)), (15, 3 * 9 + 6 * 15 + 6 * 12, 189)); + let solid: Vec = (0..36).map(|i| i as f64).collect(); + add_block(&mut k, &layout, &[0, 1], DofSet::TRANSLATIONS, &solid).unwrap(); + let beam = vec![1.0; 144]; + add_block(&mut k, &layout, &[1, 2], DofSet::STRUCTURAL, &beam).unwrap(); + // Row Uy of node 1 (global 4): the solid's row 4 on node 0's translations and node 1's. + let (cols, vals) = k.row(4); + let at = |c: usize| vals[cols.binary_search(&c).unwrap()]; + assert_eq!([at(0), at(1), at(2)], [24.0, 25.0, 26.0]); + assert_eq!([at(3), at(4), at(5)], [27.0 + 1.0, 28.0 + 1.0, 29.0 + 1.0]); + assert_eq!(at(6), 1.0, "rotation Rx of node 1 from the beam only"); + // A translations element cannot claim rotations its node lacks. + assert!(add_block(&mut k, &layout, &[0, 1], DofSet::STRUCTURAL, &beam).is_err()); + } } diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 45b16d7..3d8d828 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -2,11 +2,12 @@ //! formulations, the coupled reference points, the degrees of freedom, prescribed values and the //! load vector of a step. //! -//! Degrees of freedom: three translations per mesh node (nodes by label), then six per coupled -//! reference point (points by name): translations along X, Y, Z and rotations about them. +//! Degrees of freedom ([`DofLayout`]): the kinds each mesh node's elements declare (three +//! translations for solids), nodes by label, then six per coupled reference point (points by +//! name): translations along X, Y, Z and rotations about them. use modelcraft_fem::element::{Formulation, Solid}; -use modelcraft_fem::{Elasticity, elasticity, face}; +use modelcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, elasticity, face}; use std::collections::BTreeSet; use modelcraft_model::{ElementKind, Load, Model, Resolved, SectionKind, Step, Target}; @@ -41,18 +42,18 @@ pub struct Point { } /// The analysed mesh: the nodes that carry degrees of freedom (those of the solid elements), -/// sorted by label, the elements sorted by label, and the coupled reference points by name. +/// sorted by label, the elements sorted by label, the coupled reference points by name, and the +/// numbering of their degrees of freedom (mesh node `n` is layout slot `n`, reference point `k` +/// is [`DofLayout::point_slot`]`(k)`). #[derive(Debug, Clone)] pub struct Mesh { pub labels: Vec, pub xyz: Vec<[f64; 3]>, pub elements: Vec, pub points: Vec, + pub layout: DofLayout, } -/// The directions of a node's degrees of freedom, 0–5. -const DIRECTIONS: [&str; 6] = ["X", "Y", "Z", "rotation about X", "rotation about Y", "rotation about Z"]; - impl Mesh { pub fn node(&self, label: u32) -> Option { self.labels.binary_search(&label).ok() @@ -66,43 +67,35 @@ impl Mesh { self.points.binary_search_by(|p| p.name.as_str().cmp(name)).ok() } - /// The translational degrees of freedom of the mesh nodes. + /// The degrees of freedom of the mesh nodes; the reference points' follow them. pub fn mesh_dofs(&self) -> usize { - 3 * self.labels.len() + self.layout.node_dofs() } - /// Every degree of freedom: the mesh's, then six per reference point. + /// Every degree of freedom. pub fn dofs(&self) -> usize { - self.mesh_dofs() + 6 * self.points.len() + self.layout.dofs() } - /// Degree of freedom `c` (0–5) of reference point `k`. - pub fn point_dof(&self, k: usize, c: usize) -> usize { - self.mesh_dofs() + 6 * k + c + /// The degree of freedom of `kind` at layout slot `slot` (a mesh node or reference point). + pub fn dof(&self, slot: usize, kind: DofKind) -> Result { + self.layout.dof(slot, kind).ok_or_else(|| Error::Internal(format!("slot {slot} has no {kind:?} degree of freedom"))) } - /// The direction of a degree of freedom: 0–2 translations along X, Y, Z, 3–5 rotations. - pub fn direction(&self, dof: usize) -> usize { - if dof < self.mesh_dofs() { dof % 3 } else { (dof - self.mesh_dofs()) % 6 } - } - - /// The first degree of freedom and their count (3 or 6) of a node or coupled reference - /// point label. - pub fn label_dofs(&self, label: u32) -> Option<(usize, usize)> { - match self.node(label) { - Some(n) => Some((3 * n, 3)), - None => self.points.iter().position(|p| p.label == label).map(|k| (self.point_dof(k, 0), 6)), - } + /// The kind of a degree of freedom. + pub fn kind(&self, dof: usize) -> Option { + self.layout.locate(dof).map(|(_, kind)| kind) } /// What a degree of freedom is, in the user's terms. pub fn describe_dof(&self, dof: usize) -> String { - let direction = DIRECTIONS.get(self.direction(dof)).copied().unwrap_or("?"); - if dof < self.mesh_dofs() { - format!("node {} in {direction}", self.labels.get(dof / 3).copied().unwrap_or(0)) - } else { - let name = self.points.get((dof - self.mesh_dofs()) / 6).map_or("?", |p| p.name.as_str()); - format!("reference point `{name}` in {direction}") + let Some((slot, kind)) = self.layout.locate(dof) else { return format!("degree of freedom {dof}") }; + match self.labels.get(slot) { + Some(label) => format!("node {label} in {}", kind.describe()), + None => { + let name = self.points.get(slot - self.labels.len()).map_or("?", |p| p.name.as_str()); + format!("reference point `{name}` in {}", kind.describe()) + } } } @@ -141,7 +134,15 @@ impl Mesh { Ok(Point { name: name.to_string(), label: p.label, xyz: p.xyz }) }) .collect::>()?; - let dofs = 3 * labels.len() + 6 * points.len(); + let mut layout = DofLayout::builder(); + for e in &elements { + layout.element(&e.2, e.3.dofs()); + } + for p in &points { + layout.point(p.label, DofSet::STRUCTURAL); + } + let layout = layout.build(); + let dofs = layout.dofs(); if dofs > max_dofs { return Err(Error::Limit(format!("the model has {dofs} degrees of freedom; the limit is {max_dofs}"))); } @@ -149,7 +150,7 @@ impl Mesh { .iter() .map(|l| model.node(*l).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("node {l} is missing")))) .collect::>()?; - let mut mesh = Mesh { labels, xyz, elements: Vec::with_capacity(elements.len()), points }; + let mut mesh = Mesh { labels, xyz, elements: Vec::with_capacity(elements.len()), points, layout }; for (label, kind, nodes, formulation, d, density) in elements { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); @@ -183,10 +184,10 @@ pub fn prescribed( .push(format!("boundary condition `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh; it holds nothing")); continue; } - for (first, count) in resolved.target_nodes(bc.target()).iter().filter_map(|l| mesh.label_dofs(*l)) { - for dof in 1..=count as u8 { - let Some(v) = bc.value(dof) else { continue }; - let index = first + usize::from(dof - 1); + for slot in resolved.target_nodes(bc.target()).iter().filter_map(|l| mesh.layout.slot(*l)) { + for kind in mesh.layout.kinds(slot).iter() { + let Some(v) = bc.value(kind.number()) else { continue }; + let index = mesh.dof(slot, kind)?; if !active.get(index).copied().unwrap_or(false) { continue; } @@ -207,14 +208,27 @@ pub fn prescribed( Ok(out) } -fn add(f: &mut [f64], first: usize, v: [f64; 3]) { - for (k, c) in v.iter().enumerate() { - if let Some(x) = f.get_mut(first + k) { +/// Add a vector to the three degrees of freedom `dofs`. +fn add(f: &mut [f64], dofs: [usize; 3], v: [f64; 3]) { + for (d, c) in dofs.iter().zip(v) { + if let Some(x) = f.get_mut(*d) { *x += c; } } } +/// The translations (or rotations) of layout slot `slot`, if it has them. +fn vector_dofs(mesh: &Mesh, slot: usize, rotation: bool) -> Option<[usize; 3]> { + let kind = |c| if rotation { DofKind::rotation(c) } else { DofKind::translation(c) }; + let [a, b, c] = [0, 1, 2].map(|c| kind(c).and_then(|k| mesh.layout.dof(slot, k))); + Some([a?, b?, c?]) +} + +/// The translations of mesh node `node`. +fn translations(mesh: &Mesh, node: usize) -> Result<[usize; 3]> { + vector_dofs(mesh, node, false).ok_or_else(|| Error::Internal(format!("mesh node {node} has no translations"))) +} + /// The load vector of a static step: nodal forces and tractions as the model distributes them, /// forces and moments on coupled reference points, consistent pressure loads, gravity and body /// forces. @@ -230,15 +244,14 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ return Err(Error::Invalid(format!("load `{name}` acts on reference point `{rp}`, which no coupling ties to the mesh"))); } // Moments act on the rotations, which only reference points have. - let offset = if matches!(load, Load::Moment { .. }) { 3 } else { 0 }; + let rotation = matches!(load, Load::Moment { .. }); for (label, v) in resolved.nodal_loads.get(name).map(Vec::as_slice).unwrap_or_default() { - let (first, count) = - mesh.label_dofs(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: node {label} belongs to no solid element")))?; - if offset + 3 > count { - return Err(Error::Invalid(format!("moment `{name}`: solid nodes have no rotations; apply it to a coupled reference point"))); - } - for (c, x) in v.iter().enumerate() { - let dof = first + offset + c; + 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")) + })?; + for (&dof, x) in dofs.iter().zip(v) { if *x != 0.0 && !active.get(dof).copied().unwrap_or(false) { return Err(Error::Invalid(format!( "load `{name}`: {} is not coupled to the mesh (its coupling leaves that direction free)", @@ -246,7 +259,7 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ ))); } } - add(&mut f, first + offset, *v); + add(&mut f, dofs, *v); } } Load::Pressure { set, magnitude } => { @@ -262,7 +275,7 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ let xyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); let forces = face::pressure(&xyz, *magnitude).map_err(|e| Error::element(fc.element, e))?; for (i, v) in idx.iter().zip(forces) { - add(&mut f, 3 * i, v); + add(&mut f, translations(mesh, *i)?, v); } } } @@ -273,13 +286,13 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ continue; } let rho = e.density.ok_or_else(|| Error::Invalid(format!("load `{name}`: element {} has no density", e.label)))?; - body(&mut f, e, acceleration.map(|a| rho * a))?; + body(&mut f, mesh, e, acceleration.map(|a| rho * a))?; } } Load::BodyForce { set, vector } => { for label in resolved.sets.get(set).map(|r| r.elements.as_slice()).unwrap_or_default() { let e = mesh.element(*label).ok_or_else(|| Error::Invalid(format!("load `{name}`: element {label} is not a solid element")))?; - body(&mut f, e, *vector)?; + body(&mut f, mesh, e, *vector)?; } } } @@ -287,10 +300,10 @@ pub fn loads(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step, activ Ok(f) } -fn body(f: &mut [f64], e: &Elem, b: [f64; 3]) -> Result<()> { +fn body(f: &mut [f64], mesh: &Mesh, e: &Elem, b: [f64; 3]) -> Result<()> { let forces = e.solid()?.body_force(b).map_err(|err| Error::element(e.label, err))?; for (node, v) in e.nodes.iter().zip(forces) { - add(f, 3 * node, v); + add(f, translations(mesh, *node)?, v); } Ok(()) } diff --git a/docs/elements.md b/docs/elements.md index ef0b727..0ede686 100644 --- a/docs/elements.md +++ b/docs/elements.md @@ -9,6 +9,12 @@ then shells, special-purpose and multiphysics elements. Status: **done** (on main, verified against analytic solutions), **in progress** (an agent task is open), **planned**. +Every element builds on the shared foundation of docs/fem.md: it declares its nodal degrees of +freedom (`DofKind`: translations, rotations, temperature, pressure) and condensed internal ones +(`ElementDofs`), the solvers number them with `DofLayout`, and its interpolation comes from the +shape-function library (`Basis`: line2/3, tri3/6, quad4/8/9, tet4/10, wedge6/15, pyramid5, +hex8/20/27 with first and second derivatives). + ## Point and discrete elements | Element | Use | Status | @@ -66,8 +72,8 @@ open), **planned**. | hex8 (full, incompatible modes), hex20 (full, reduced) | Structured and swept meshes | done | | wedge6 | Swept meshes from triangles, transitions | done | | hex8 one-point with hourglass control, nodal-pressure tet4 (explicit) | Explicit dynamics | done | -| wedge15, pyramid5/13 | Quadratic transitions, tet–hex transition | planned (M5) | -| Hybrid (mixed u–p) solids | Nearly incompressible materials (rubber) | planned (M8) | +| wedge15, pyramid5/13 | Quadratic transitions, tet–hex transition | planned (M5); wedge15 and pyramid5 shape functions done | +| Hybrid (mixed u–p) solids | Nearly incompressible materials (rubber) | in progress (FEM foundation, phase C); pressure DOF kind done | | Modified quadratic tet (tet10 for contact/plasticity) | Robust quadratic tets | planned | ## Special-purpose and multiphysics @@ -77,6 +83,6 @@ open), **planned**. | Cohesive (2D, 3D) | Adhesive layers, delamination | planned | | Gasket | Through-thickness gasket behaviour | planned | | Infinite elements | Unbounded domains | planned | -| Heat transfer (all shapes), coupled temperature–displacement | Thermal and thermal-stress analyses | planned (M8) | +| Heat transfer (all shapes), coupled temperature–displacement | Thermal and thermal-stress analyses | in progress (FEM foundation, phase B); temperature DOF kind done | | Acoustic | Fluid acoustics, structural–acoustic coupling | planned | | Rigid elements | Rigid bodies from meshes | planned (M7) | diff --git a/docs/fem.md b/docs/fem.md new file mode 100644 index 0000000..38d0a86 --- /dev/null +++ b/docs/fem.md @@ -0,0 +1,97 @@ +# FEM foundation: degrees of freedom and shape functions + +What every element in `crates/fem` and every solver builds on. Code: `crates/fem/src/dof.rs`, +`crates/fem/src/basis.rs`, `crates/fem/src/jet.rs`. + +## Degrees of freedom + +Every unknown is a pair (node, `DofKind`): + +| `DofKind` | Meaning | Model / deck number | +|---|---|---| +| `Ux`, `Uy`, `Uz` | translations along X, Y, Z | 1, 2, 3 | +| `Rx`, `Ry`, `Rz` | rotations about X, Y, Z | 4, 5, 6 | +| `P` | pressure (mixed u–p elements) | 8 | +| `T` | temperature | 11 | + +The enum order (Ux, Uy, Uz, Rx, Ry, Rz, T, P) is the numbering order within a node. +`DofKind::from_number` / `number` convert to the 1-based numbers boundary conditions use; +`from_direction` / `direction` to the structural directions 0–5. + +**Elements declare their unknowns** with `ElementDofs { node: DofSet, internal: Internal }`: +the kinds every node of the element carries (`DofSet::TRANSLATIONS` for solids, `PLANE` for +2D elements, `STRUCTURAL` for shells and beams, `TEMPERATURE` for heat transfer, translations +plus `PRESSURE` for continuous-pressure mixed elements) and the internal unknowns the element +condenses out itself (`incompatible_modes`, discontinuous `pressure`), which never reach the +global system. An element's matrices are ordered node by node, its `node` kinds in `DofKind` +order within a node (a shell's 24 × 24 stiffness: node 1 Ux Uy Uz Rx Ry Rz, node 2 …). Solid +elements: `Formulation::dofs()` / `Solid::dofs()`. + +**Numbering** (`DofLayout`): collect what every node needs with `DofLayout::builder()` +(`element(labels, dofs)`, `node(label, kinds)` for constraints and boundary conditions that +need a kind no element gives, `point(label, kinds)` for reference points), then `build()`. A +node carries the union of what was asked for. Unknowns are numbered mesh nodes in label order, +each node's kinds in `DofKind` order, then the points in the order added; a slot's unknowns are +contiguous, so a node pair is one dense block in the assembled matrix. The same model always +gives the same numbering. + +Lookups: `node(label)` / `slot(label)` (a slot is a position: mesh nodes `0..nodes()`, points +`point_slot(k)`), `dof(slot, kind)`, `range(slot)`, `kinds(slot)`, `locate(dof) -> (slot, +kind)`, `element_dofs(slots, kinds)`, `describe(dof)`, `stride()` (uniform layouts). + +Solvers: `modelcraft_solver::system::Mesh::layout` numbers the implicit system (mesh nodes, +then six per coupled reference point); assembly (`sparse::add_block(k, layout, slots, kinds, +data)`, `add_scalar`, `add_diagonal`) works for any mix of kinds per node, so a solid sharing +nodes with a shell adds its translation block into the shell nodes' six-unknown blocks. +Boundary conditions and loads address (node, kind): a boundary condition holds the kinds a node +carries, and a moment needs rotations. The explicit solver has solid nodes only and checks that +its layout is the uniform node-major one (`NODE_DOFS = 3`). + +## Shape functions + +`Basis` is every family's interpolation: + +| Family | Reference element | Nodes | +|---|---|---| +| `Line2`, `Line3` | ξ ∈ [−1, 1] | ends −1, 1; line3 adds the middle | +| `Tri3`, `Tri6` | ξ, η ≥ 0, ξ + η ≤ 1 | corners (0,0), (1,0), (0,1); mid-sides of 1-2, 2-3, 3-1 | +| `Quad4`, `Quad8`, `Quad9` | [−1, 1]² | corners counter-clockwise from (−1,−1); mid-sides of 1-2, 2-3, 3-4, 4-1; centre | +| `Tet4`, `Tet10` | unit tetrahedron | corners (0,0,0), (1,0,0), (0,1,0), (0,0,1); mid-sides of 1-2, 2-3, 3-1, 1-4, 2-4, 3-4 | +| `Wedge6`, `Wedge15` | triangle × ζ ∈ [−1, 1] | bottom 1-3, top 4-6; mid-sides of the bottom edges, top edges, vertical edges | +| `Pyramid5` | base [−1, 1]² at ζ = 0, apex (0, 0, 1) | base corners as quad4, then the apex | +| `Hex8`, `Hex20`, `Hex27` | [−1, 1]³ | corners 1-4 at ζ = −1, 5-8 above; hex20 mid-sides (bottom, top, vertical edges); hex27 face centres in `ElementKind::Hex8.faces()` order, then the centre | + +Node orders are the model's (`ElementKind::edges()` / `faces()` are the single source): +corners, mid-side nodes in edge order, face nodes, body nodes. `Basis::of(ElementKind)` maps a +model element; `nodes()` gives natural node coordinates; `dim()`, `count()`, `corners()`, +`degree()` (the complete polynomial degree reproduced: 1 for linear families including the +rational pyramid, 2 for the quadratic ones). + +Evaluation at a natural point `xi: [f64; 3]` (unused coordinates are zero): + +- `eval(xi) -> Shape`: values `n` and first derivatives `dn[i] = ∂Nᵢ/∂(ξ, η, ζ)`. tet4, tet10, + wedge6, hex8, hex20 use the hand-written functions of `shape.rs` (the solid kernels' path). +- `eval2(xi) -> Shape2`: adds `d2n[i][a][b] = ∂²Nᵢ/∂ξₐ∂ξ_b`, from second-order forward-mode + automatic differentiation of the same formulas (`jet::Jet`): exact to round-off. +- `gradients(&shape, xyz) -> Gradients`: physical `dx[i] = ∇ₓNᵢ = J⁻¹ ∇_ξNᵢ` and `det J`, with + J_ab = ∂x_b/∂ξₐ; an inverted or collapsed element is `Error::Degenerate`. +- `hessians(&shape2, xyz) -> Hessians`: physical `d2x[i]` = + J⁻¹ (H_ξ − Σ_b (∂N/∂x_b) X_b″) J⁻ᵀ with X_b″ = Σᵢ x_{i,b} H_ξ(Nᵢ), the map's second + derivative (T. J. R. Hughes, *The Finite Element Method*, 2000, §3.9, §4.2), and + `laplacian(i)`. Without the X″ term a field linear in x would have a nonzero Hessian on a + curved element. + +A family of dimension d maps onto the first d physical coordinates (plane elements in x–y, +lines along x). Elements embedded in a higher dimension (shells, beams, trusses in 3D) pass +node coordinates in their own local frame. + +## Verification (`basis::tests`, `dof::tests`, `jet::tests`) + +For every family: Kronecker delta at the nodes; partition of unity (Σ N = 1, Σ ∇N = 0, +Σ H = 0); exact reproduction of every monomial up to `degree()` with its gradient and Hessian +in natural coordinates; first and second derivatives against central differences; +on an affine map, physical gradients, Hessians and Laplacians exact for a complete polynomial +of the family's degree; on a curved map, linear fields exact with zero Hessian, and the +physical Hessian of a quadratic interpolant consistent with finite differences of its physical +gradient. The solver and explicit migration to `DofLayout` was checked bit for bit: every +result of the solver, explicit, engine and io test suites hashes identically before and after.