From 5ebb9f16dcd3a69f8eb895ac97cdd4b6ad455fc0 Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Tue, 13 Oct 2026 22:20:21 -0400 Subject: [PATCH] =?UTF-8?q?Hybrid=20(mixed=20u=E2=80=93p)=20solids:=20hex8?= =?UTF-8?q?=20and=20tet10=20with=20a=20constant=20pressure,=20hex20=20with?= =?UTF-8?q?=20a=20linear=20one=20(#49)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The isotropic elasticity splits into D_dev + K m mᵀ; an element-wise discontinuous pressure is condensed out (p̂ = K M⁻¹ Gᵀ u), leaving K_dev + K G M⁻¹ Gᵀ (Hughes §4.3–4.5). Q2–P1disc and P2–P0 are inf-sup stable; Q1–P0 does not lock; tet4 and wedge6 still lock and are refused, as are orthotropic materials (named at the section). ElementDofs counts the condensed pressures. Verified: element patch, rigid modes and totals for every hybrid kind, a block in uniaxial strain at ν = 0.4999 (exact), Lamé's cylinder in plane strain at ν = 0.4999 (hex20 1e-5, tet10 1.3e-3, hex8 2.5e-3 vs the locking displacement hex8 at −80 %). Every non-hybrid result of the test suites hashes identically before and after. Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 --- ROADMAP.md | 2 +- crates/fem/src/element.rs | 75 ++++++--- crates/fem/src/hybrid.rs | 277 ++++++++++++++++++++++++++++++++++ crates/fem/src/lib.rs | 1 + crates/fem/src/plane.rs | 1 + crates/fem/src/tests.rs | 21 ++- crates/solver/src/system.rs | 8 +- crates/solver/tests/hybrid.rs | 140 +++++++++++++++++ docs/elements.md | 2 +- docs/fem.md | 16 ++ docs/solvers/literature.md | 14 ++ 11 files changed, 525 insertions(+), 32 deletions(-) create mode 100644 crates/fem/src/hybrid.rs create mode 100644 crates/solver/tests/hybrid.rs diff --git a/ROADMAP.md b/ROADMAP.md index ff7477f..869756b 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -26,7 +26,7 @@ hex meshes today and it gives the best elements where stresses are read. | M5 Structured hex | 2D quad meshing of faces (`modelcraft-mesh::quad`: transfinite mapped meshes, free all-quad meshes of polygons with holes on prescribed boundary nodes, exact planar quad validity, shared `SizeField` and even edge divisions; done for planar faces); mapped hex of six-sided regions; sweep meshing; conforming interfaces and tet–hex transition; automatic decomposition into sweepable and mappable regions, proposed as SolidCraft splits | started | | M6 Automatic hex | Octree/grid-dual all-hex for any part; hex-dominant by tet recombination; polycube for blocky parts; smoothing, untangling, boundary padding, layers through thin walls; mesh editing and orphan-mesh import | planned | | M7 Assemblies and contact | Part instances from SolidCraft occurrences; tie between non-matching meshes; general and pairwise contact with automatic pair detection, friction; bolt pretension, connectors; amplitudes, multiple steps; nonlinear static steps in decks | planned | -| M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; heat transfer (solids done: steady and θ-method transient conduction with flux, convection and sources, thermal strain in static steps (solids and 2D) and explicit dynamics (with plasticity), coupled temperature–displacement, `solve.heat`; deck output, 2D/shell/beam conduction planned) | started | +| M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; hybrid (mixed u–p) hex8/hex20/tet10 for nearly incompressible materials (done); heat transfer (solids done: steady and θ-method transient conduction with flux, convection and sources, thermal strain in static steps (solids and 2D) and explicit dynamics (with plasticity), coupled temperature–displacement, `solve.heat`; deck output, 2D/shell/beam conduction planned) | started | | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned | | M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `modelcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (X2: done — updated Lagrangian elements, Jaumann/Hughes–Winget and Green–Naghdi/Flanagan–Taylor rates, J2 bilinear/table plasticity, Johnson–Cook strength with heating and damage, erosion; materials gain `plastic`, `damage`, `erosion` with `model.material.plasticity`/`damage`; Dienes shear, rigid rotation, J2 uniaxial and Taylor-bar tests); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | started | | M11 Equations of state and shells | Polynomial, ideal-gas, Mie–Grüneisen, JWL with programmed burn, tabulated; energy coupling; shock verification (Rankine–Hugoniot, Sod, Noh, Sedov); Belytschko–Lin–Tsay shells with self-contact on quad meshes | planned | diff --git a/crates/fem/src/element.rs b/crates/fem/src/element.rs index 790c192..6c1f1a7 100644 --- a/crates/fem/src/element.rs +++ b/crates/fem/src/element.rs @@ -146,6 +146,11 @@ pub enum Formulation { Full, Reduced, IncompatibleModes, + /// Mixed displacement–pressure with `pressures` discontinuous pressure unknowns condensed + /// out of the element ([`crate::hybrid`]). + Hybrid { + pressures: usize, + }, } impl Formulation { @@ -153,14 +158,23 @@ impl Formulation { /// 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 } } + let pressure = if let Formulation::Hybrid { pressures } = self { pressures } else { 0 }; + ElementDofs { node: DofSet::TRANSLATIONS, internal: Internal { incompatible_modes, pressure } } } /// The formulation a model section asks for, if this element type has it. pub fn from_section(kind: ElementKind, f: &modelcraft_model::Formulation) -> Result { let ty = element_type(kind)?; if f.hybrid { - return Err(Error::NotSupported("hybrid (mixed displacement–pressure) elements".into())); + if f.incompatible_modes || f.integration == modelcraft_model::Integration::Reduced { + return Err(Error::NotSupported("hybrid elements with incompatible modes or reduced integration".into())); + } + return crate::hybrid::pressures(kind).map(|pressures| Formulation::Hybrid { pressures }).ok_or_else(|| { + Error::NotSupported(format!( + "hybrid {} elements (with one pressure per element the linear tetrahedron and wedge still lock; use hex8, hex20 or tet10)", + kind.name() + )) + }); } if f.incompatible_modes { return if ty.incompatible_modes { @@ -189,10 +203,10 @@ pub struct PointStress { } /// Shape function gradients in physical coordinates at one point. -struct Gradients { - shape: Shape, - dx: [[f64; 3]; MAX_NODES], - det: f64, +pub(crate) struct Gradients { + pub(crate) shape: Shape, + pub(crate) dx: [[f64; 3]; MAX_NODES], + pub(crate) det: f64, } /// A thermal strain field over an element: the expansion coefficients along X, Y, Z and the @@ -219,8 +233,8 @@ impl ThermalStrain { #[derive(Debug, Clone, Copy)] pub struct Solid<'a> { pub ty: &'static ElementType, - xyz: &'a [[f64; 3]], - formulation: Formulation, + pub(crate) xyz: &'a [[f64; 3]], + pub(crate) formulation: Formulation, } /// The stiffness blocks of the incompatible-modes element. @@ -245,6 +259,7 @@ impl<'a> Solid<'a> { Formulation::Full => true, Formulation::Reduced => ty.reduced.is_some(), Formulation::IncompatibleModes => ty.incompatible_modes, + Formulation::Hybrid { pressures } => crate::hybrid::pressures(kind) == Some(pressures), }; if !ok { return Err(Error::NotSupported(format!("{formulation:?} formulation of {} elements", kind.name()))); @@ -283,7 +298,7 @@ impl<'a> Solid<'a> { j } - fn gradients(&self, xi: [f64; 3]) -> Result { + pub(crate) fn gradients(&self, xi: [f64; 3]) -> Result { let shape = (self.ty.shape)(xi); let j = self.jacobian(&shape); let (inv, det) = positive_inverse(&j).ok_or_else(|| self.degenerate())?; @@ -302,11 +317,15 @@ impl<'a> Solid<'a> { /// Position of a natural point. pub fn position(&self, xi: [f64; 3]) -> [f64; 3] { - let s = (self.ty.shape)(xi); + self.position_of((self.ty.shape)(xi).values()) + } + + /// The position Σ Nᵢ xᵢ for shape values `n`. + pub(crate) fn position_of(&self, n: &[f64]) -> [f64; 3] { let mut p = [0.0; 3]; - for (n, x) in s.values().iter().zip(self.xyz) { + for (ni, x) in n.iter().zip(self.xyz) { for k in 0..3 { - p[k] += n * x[k]; + p[k] += ni * x[k]; } } p @@ -364,6 +383,9 @@ impl<'a> Solid<'a> { /// The stiffness matrix, 3n × 3n, degrees of freedom ordered node by node (x, y, z). pub fn stiffness(&self, d: &Elasticity) -> Result { + if let Formulation::Hybrid { .. } = self.formulation { + return self.hybrid_stiffness(d); + } let n = 3 * self.ty.nodes; if self.formulation == Formulation::IncompatibleModes { let IncompatibleParts { kuu: mut k, kua, kaa, .. } = self.incompatible_parts(d, true)?; @@ -444,7 +466,7 @@ impl<'a> Solid<'a> { } /// The thermal strain α ΔT at a point with shape values `n` (xx, yy, zz; no shear). - fn thermal_strain(thermal: &ThermalStrain, n: &[f64]) -> [f64; 6] { + pub(crate) fn thermal_strain(thermal: &ThermalStrain, n: &[f64]) -> [f64; 6] { let dt: f64 = n.iter().zip(&thermal.delta).map(|(a, b)| a * b).sum(); [thermal.alpha[0] * dt, thermal.alpha[1] * dt, thermal.alpha[2] * dt, 0.0, 0.0, 0.0] } @@ -454,6 +476,9 @@ impl<'a> Solid<'a> { /// are condensed out as in the stiffness, f = f_u − K_ua K_aa⁻¹ f_a. pub fn thermal_load(&self, d: &Elasticity, thermal: &ThermalStrain) -> Result> { thermal.check(self.ty.nodes)?; + if let Formulation::Hybrid { .. } = self.formulation { + return self.hybrid_thermal_load(d, thermal); + } if self.formulation == Formulation::IncompatibleModes { let IncompatibleParts { kua, kaa, j0, .. } = self.incompatible_parts(d, false)?; let (mut f, mut fa) = self.thermal_forces(d, thermal, Some(&j0))?; @@ -511,6 +536,8 @@ impl<'a> Solid<'a> { if let Some(t) = thermal { t.check(self.ty.nodes)?; } + // Hybrid: the element's pressure field, and σ = D_dev (ε − ε_th) + p m. + let hybrid = if let Formulation::Hybrid { .. } = self.formulation { Some(self.pressure_field(d, u, thermal)?) } else { None }; // Incompatible modes: recover the internal amplitudes from K_aa α + K_au u = f_a (the // thermal load on the modes; zero without one). let alpha = if self.formulation == Formulation::IncompatibleModes { @@ -552,18 +579,18 @@ impl<'a> Solid<'a> { } } let strain = [grad[0][0], grad[1][1], grad[2][2], grad[0][1] + grad[1][0], grad[0][2] + grad[2][0], grad[1][2] + grad[2][1]]; - let mut xyz = [0.0; 3]; - for (n, x) in g.shape.values().iter().zip(self.xyz) { - for k in 0..3 { - xyz[k] += n * x[k]; - } - } - let sigma = match thermal { + let xyz = self.position_of(g.shape.values()); + // The mechanical strain, less the thermal strain when there is one. + let mechanical = match thermal { Some(t) => { let eth = Self::thermal_strain(t, g.shape.values()); - stress(d, &[0, 1, 2, 3, 4, 5].map(|k| strain[k] - eth[k])) + [0, 1, 2, 3, 4, 5].map(|k| strain[k] - eth[k]) } - None => stress(d, &strain), + None => strain, + }; + let sigma = match &hybrid { + Some(h) => h.stress(&mechanical, xyz), + None => stress(d, &mechanical), }; out.push(PointStress { xyz, strain, stress: sigma }); } @@ -606,7 +633,7 @@ pub(crate) fn symmetrize(k: &mut Matrix) { } /// The strain-displacement block of one gradient g: ε = B u for a single node's (x, y, z). -fn b_block(g: &[f64; 3]) -> [[f64; 3]; 6] { +pub(crate) fn b_block(g: &[f64; 3]) -> [[f64; 3]; 6] { [[g[0], 0.0, 0.0], [0.0, g[1], 0.0], [0.0, 0.0, g[2]], [g[1], g[0], 0.0], [g[2], 0.0, g[0]], [0.0, g[2], g[1]]] } @@ -638,7 +665,7 @@ fn accumulate(out: &mut [f64], cols: usize, ga: &[[f64; 3]], gb: &[[f64; 3]], d: /// The symmetric case of [`accumulate`] (rows and columns from the same gradients, `out` square): /// computes the blocks with j ≥ i and mirrors them, so the result is exactly symmetric. -fn accumulate_symmetric(out: &mut [f64], g: &[[f64; 3]], d: &Elasticity, w: f64) { +pub(crate) fn accumulate_symmetric(out: &mut [f64], g: &[[f64; 3]], d: &Elasticity, w: f64) { let cols = 3 * g.len(); for (j, gj) in g.iter().enumerate() { let bj = b_block(gj); diff --git a/crates/fem/src/hybrid.rs b/crates/fem/src/hybrid.rs new file mode 100644 index 0000000..4dc3a98 --- /dev/null +++ b/crates/fem/src/hybrid.rs @@ -0,0 +1,277 @@ +//! Hybrid (mixed displacement–pressure) solids for nearly incompressible materials. +//! +//! The isotropic elasticity splits into a deviatoric part and the bulk modulus, +//! D = D_dev + K m mᵀ with m = (1, 1, 1, 0, 0, 0) and K = mᵀ D m / 9. The pressure (here the +//! mean stress, positive in tension) is an independent field p = N_p · p̂, discontinuous between +//! elements, and the volumetric strain is enforced weakly (the perturbed Lagrangian of +//! T. J. R. Hughes, *The Finite Element Method*, Dover 2000, §4.3–4.4): +//! +//! K_dev = ∫ Bᵀ D_dev B dV, G = ∫ Bᵀ m N_pᵀ dV, M = ∫ N_p N_pᵀ dV, +//! [K_dev G; Gᵀ −M/K] [u; p̂] = [f; 0], +//! +//! so p̂ = K M⁻¹ Gᵀ u is condensed out element by element and the element stiffness is +//! K_dev + K G M⁻¹ Gᵀ: symmetric positive definite for any ν < ½, and bounded volumetric +//! stiffness per element however close ν is to ½ — the constraint count no longer grows with +//! the nodes, which is what locks the displacement elements. With a constant pressure this is +//! the B-bar / mean dilatation element (Hughes 2000, §4.5; J. C. Nagtegaal, D. M. Parks, +//! J. R. Rice, *On numerically accurate finite element solutions in the fully plastic range*, +//! CMAME 4(2):153–177, 1974). +//! +//! Pressure fields and stability (the inf-sup or LBB condition, F. Brezzi, M. Fortin, *Mixed +//! and Hybrid Finite Element Methods*, Springer 1991, ch. VI; K.-J. Bathe, *Finite Element +//! Procedures*, 2nd ed., 2014, §4.5): +//! - **hex20**: linear pressure 1, x, y, z (Q2–P1disc, inf-sup stable); +//! - **tet10**: constant pressure (P2–P0, inf-sup stable; P2–P1disc is not on general meshes); +//! - **tet10**: a single element has three more zero-energy modes (dilatations with a linear +//! volumetric strain, which its deviatoric stiffness does not see and one pressure controls +//! only in the mean); neighbouring elements restrain them, as they do the reduced hex20's; +//! - **hex8**: constant pressure (Q1–P0). It passes the patch test and does not lock, but is +//! not inf-sup stable: on some meshes and constraints the element pressures show a +//! checkerboard mode while displacements stay good; nodal stresses are averaged anyway; +//! - **tet4, wedge6**: refused — their cross-section is the linear triangle, and one pressure +//! per element still leaves too few displacement unknowns per constraint: they lock (verified +//! on the thick cylinder; a nodal-pressure formulation would be needed). +//! +//! A continuous nodal pressure (Taylor–Hood, P on the corner nodes) makes the global system a +//! symmetric indefinite saddle point that the Cholesky factorization cannot solve; it waits for +//! an indefinite (LDLᵀ) factorization. + +use modelcraft_model::ElementKind; + +use crate::dense::{Cholesky, Matrix}; +use crate::element::{Formulation, Solid, ThermalStrain, accumulate_symmetric, b_block}; +use crate::material::Elasticity; +use crate::{Error, Result}; + +/// The pressure unknowns of the hybrid formulation of an element kind, if it has one. +pub fn pressures(kind: ElementKind) -> Option { + match kind { + ElementKind::Hex20 => Some(4), + ElementKind::Hex8 | ElementKind::Tet10 => Some(1), + _ => None, + } +} + +/// The deviatoric elasticity and the bulk modulus of an isotropic elasticity; orthotropic +/// materials have no such split and are refused. +pub fn split(d: &Elasticity) -> Result<(Elasticity, f64)> { + let scale = d.iter().flatten().fold(0.0f64, |m, v| m.max(v.abs())); + let same = |a: f64, b: f64| (a - b).abs() <= 1e-12 * scale; + let (c11, c12, c44) = (d[0][0], d[0][1], d[3][3]); + let isotropic = (0..3).all(|i| same(d[i][i], c11) && same(d[i + 3][i + 3], c44)) + && (0..3).all(|i| (0..3).all(|j| i == j || same(d[i][j], c12))) + && same(c44, 0.5 * (c11 - c12)) + && (0..6).all(|i| (0..6).all(|j| (i < 3 && j < 3) || i == j || d[i][j] == 0.0)); + if !isotropic { + return Err(Error::NotSupported("hybrid elements of orthotropic materials".into())); + } + let bulk = (c11 + 2.0 * c12) / 3.0; + let mut dev = *d; + for row in dev.iter_mut().take(3) { + for v in row.iter_mut().take(3) { + *v -= bulk; + } + } + Ok((dev, bulk)) +} + +/// The parts of a hybrid element: G, 3n × np row-major; the factor of M; and K_dev when +/// asked for. +struct Mixed { + kdev: Option, + g: Vec, + m: Cholesky, + /// With a thermal strain: ∫ N_p tr ε_th dV (np) and ∫ Bᵀ D_dev ε_th dV (3n). + h: Vec, + f_dev: Vec, +} + +/// The pressure field of a hybrid element under nodal displacements, for its stresses. +pub(crate) struct PressureField { + dev: Elasticity, + bulk: f64, + /// M⁻¹ Gᵀ u, the coefficients of the volumetric strain's projection. + coefficients: [f64; 4], + np: usize, + frame: ([f64; 3], f64), +} + +impl PressureField { + /// D_dev ε + p m at `x`. + pub(crate) fn stress(&self, strain: &[f64; 6], x: [f64; 3]) -> [f64; 6] { + let basis = pressure_basis(self.np, x, &self.frame); + let p = self.bulk * self.coefficients.iter().zip(basis).map(|(a, b)| a * b).sum::(); + let mut s = crate::material::stress(&self.dev, strain); + for v in s.iter_mut().take(3) { + *v += p; + } + s + } +} + +/// The pressure functions at a point: 1, then (x − c)/h, (y − c)/h, (z − c)/h for the linear +/// field (`np` = 4). +fn pressure_basis(np: usize, x: [f64; 3], frame: &([f64; 3], f64)) -> [f64; 4] { + let (c, h) = frame; + let mut out = [1.0, 0.0, 0.0, 0.0]; + if np == 4 && *h > 0.0 { + for k in 0..3 { + out[k + 1] = (x[k] - c[k]) / h; + } + } + out +} + +impl Solid<'_> { + /// The pressure unknowns of this (hybrid) element. + fn np(&self) -> usize { + match self.formulation { + Formulation::Hybrid { pressures } => pressures.min(4), + _ => 0, + } + } + + /// The element's centre (mean of its corners) and size (largest corner distance from it), + /// which scale the linear pressure functions. + fn frame(&self) -> ([f64; 3], f64) { + let corners = self.xyz.get(..self.ty.corners).unwrap_or_default(); + let n = corners.len().max(1) as f64; + let c = corners.iter().fold([0.0; 3], |a, x| [a[0] + x[0] / n, a[1] + x[1] / n, a[2] + x[2] / n]); + let size = corners.iter().map(|x| ((x[0] - c[0]).powi(2) + (x[1] - c[1]).powi(2) + (x[2] - c[2]).powi(2)).sqrt()).fold(0.0, f64::max); + (c, size) + } + + /// G and M (and K_dev with `kdev`, and the thermal terms with `thermal`) in one pass over + /// the stiffness rule. + fn mixed(&self, dev: &Elasticity, kdev: bool, thermal: Option<&ThermalStrain>) -> Result { + let (np, nodes, frame) = (self.np(), self.ty.nodes, self.frame()); + let mut kdev = kdev.then(|| Matrix::zeros(3 * nodes)); + let mut g = vec![0.0; 3 * nodes * np]; + let mut m = Matrix::zeros(np); + let mut h = vec![0.0; np]; + let mut f_dev = vec![0.0; 3 * nodes]; + for p in self.rule().points() { + let gr = self.gradients(p.xi)?; + let w = gr.det * p.weight; + let gu = gr.dx.get(..nodes).unwrap_or_default(); + if let Some(k) = kdev.as_mut() { + accumulate_symmetric(&mut k.data, gu, dev, w); + } + let x = if np == 4 { self.position_of(gr.shape.values()) } else { [0.0; 3] }; + let np_at = pressure_basis(np, x, &frame); + if let Some(t) = thermal { + let eth = Self::thermal_strain(t, gr.shape.values()); + let volumetric = eth[0] + eth[1] + eth[2]; + for (hk, nk) in h.iter_mut().zip(np_at) { + *hk += w * nk * volumetric; + } + let sig = crate::material::stress(dev, ð); + for (i, gi) in gu.iter().enumerate() { + let b = b_block(gi); + for c in 0..3 { + if let Some(f) = f_dev.get_mut(3 * i + c) { + *f += w * (0..6).map(|s| b[s][c] * sig[s]).sum::(); + } + } + } + } + // Bᵢᵀ m is node i's gradient. + for (i, gi) in gu.iter().enumerate() { + for c in 0..3 { + for k in 0..np { + if let Some(v) = g.get_mut((3 * i + c) * np + k) { + *v += w * gi[c] * np_at[k]; + } + } + } + } + for a in 0..np { + for b in 0..np { + m.add(a, b, w * np_at[a] * np_at[b]); + } + } + } + Ok(Mixed { kdev, g, m: Cholesky::new(&m)?, h, f_dev }) + } + + /// K_dev + K G M⁻¹ Gᵀ, symmetric by construction (H = M⁻¹ Gᵀ once, upper triangle mirrored). + pub(crate) fn hybrid_stiffness(&self, d: &Elasticity) -> Result { + let (dev, bulk) = split(d)?; + let Mixed { kdev, g, m, .. } = self.mixed(&dev, true, None)?; + let mut k = kdev.ok_or_else(|| Error::Invalid("hybrid stiffness".into()))?; + let (n, np) = (k.n, self.np()); + let row = |i: usize| g.get(i * np..(i + 1) * np).unwrap_or_default(); + let h: Vec> = (0..n) + .map(|j| { + let mut y = row(j).to_vec(); + m.solve(&mut y); + y + }) + .collect(); + for i in 0..n { + for (j, hj) in h.iter().enumerate().skip(i) { + let v = bulk * row(i).iter().zip(hj).map(|(a, b)| a * b).sum::(); + k.add(i, j, v); + if i != j { + k.add(j, i, v); + } + } + } + Ok(k) + } + + /// The equivalent nodal forces of a thermal strain: ∫ Bᵀ D_dev ε_th dV for the deviatoric + /// part, and K G M⁻¹ ∫ N_p tr ε_th dV for the volumetric one through the pressure. + pub(crate) fn hybrid_thermal_load(&self, d: &Elasticity, thermal: &ThermalStrain) -> Result> { + let (dev, bulk) = split(d)?; + let np = self.np(); + let Mixed { g, m, mut h, f_dev, .. } = self.mixed(&dev, false, Some(thermal))?; + m.solve(&mut h); + let f: Vec = f_dev + .iter() + .enumerate() + .map(|(i, fd)| fd + bulk * (0..np).map(|k| g.get(i * np + k).copied().unwrap_or(0.0) * h.get(k).copied().unwrap_or(0.0)).sum::()) + .collect(); + Ok(f.as_chunks::<3>().0.to_vec()) + } + + /// The element's pressure field under the nodal displacements `u` and a thermal strain: + /// p̂ = K M⁻¹ (Gᵀ u − ∫ N_p tr ε_th dV). + pub(crate) fn pressure_field(&self, d: &Elasticity, u: &[[f64; 3]], thermal: Option<&ThermalStrain>) -> Result { + let (dev, bulk) = split(d)?; + let np = self.np(); + let Mixed { g, m, h, .. } = self.mixed(&dev, false, thermal)?; + let mut pc: Vec = h.iter().map(|x| -x).collect(); + for (i, ui) in u.iter().flatten().enumerate() { + for (k, x) in pc.iter_mut().enumerate() { + *x += g.get(i * np + k).copied().unwrap_or(0.0) * ui; + } + } + m.solve(&mut pc); + let mut coefficients = [0.0; 4]; + for (c, v) in coefficients.iter_mut().zip(pc) { + *c = v; + } + Ok(PressureField { dev, bulk, coefficients, np, frame: self.frame() }) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::material::{isotropic, orthotropic}; + + #[test] + fn the_split_is_deviatoric_plus_bulk() { + let d = isotropic(1000.0, 0.3).unwrap(); + let (dev, k) = split(&d).unwrap(); + assert!((k - 1000.0 / (3.0 * (1.0 - 0.6))).abs() < 1e-9); + // D_dev m = 0. + for r in 0..6 { + assert!((0..3).map(|c| dev[r][c]).sum::().abs() < 1e-9); + } + assert!(split(&orthotropic([1000.0, 800.0, 600.0], [0.3, 0.25, 0.2], [400.0, 300.0, 350.0]).unwrap()).is_err()); + assert_eq!(pressures(ElementKind::Tet4), None); + assert_eq!(pressures(ElementKind::Hex20), Some(4)); + } +} diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index f0ddcdc..ad41bf0 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -50,6 +50,7 @@ pub mod element; pub mod face; pub mod frame; pub mod heat; +pub mod hybrid; pub mod jet; pub mod material; pub mod plane; diff --git a/crates/fem/src/plane.rs b/crates/fem/src/plane.rs index df23620..bd19975 100644 --- a/crates/fem/src/plane.rs +++ b/crates/fem/src/plane.rs @@ -270,6 +270,7 @@ impl<'a> Plane<'a> { Formulation::Full => true, Formulation::Reduced => ty.reduced.is_some(), Formulation::IncompatibleModes => ty.incompatible_modes, + Formulation::Hybrid { .. } => false, }; if !ok { return Err(Error::NotSupported(format!("{formulation:?} formulation of 2D {} elements", kind.name()))); diff --git a/crates/fem/src/tests.rs b/crates/fem/src/tests.rs index ac7cf9a..75d15f7 100644 --- a/crates/fem/src/tests.rs +++ b/crates/fem/src/tests.rs @@ -34,7 +34,7 @@ fn distorted(kind: ElementKind, curved: bool) -> Vec<[f64; 3]> { xyz } -const CASES: [(ElementKind, Formulation); 7] = [ +const CASES: [(ElementKind, Formulation); 10] = [ (ElementKind::Tet4, Formulation::Full), (ElementKind::Tet10, Formulation::Full), (ElementKind::Wedge6, Formulation::Full), @@ -42,6 +42,9 @@ const CASES: [(ElementKind, Formulation); 7] = [ (ElementKind::Hex8, Formulation::IncompatibleModes), (ElementKind::Hex20, Formulation::Full), (ElementKind::Hex20, Formulation::Reduced), + (ElementKind::Tet10, Formulation::Hybrid { pressures: 1 }), + (ElementKind::Hex8, Formulation::Hybrid { pressures: 1 }), + (ElementKind::Hex20, Formulation::Hybrid { pressures: 4 }), ]; fn rigid_modes(xyz: &[[f64; 3]]) -> Vec> { @@ -75,8 +78,11 @@ fn stiffness_is_symmetric_with_exactly_six_rigid_body_modes() { } // Holding six degrees of freedom (node 1 fully, node 2 in y and z, node 3 in z) leaves a // positive definite matrix: no other zero-energy mode. Reduced hex20 has spurious modes - // on a single element, so it is skipped. - if f != Formulation::Reduced { + // on a single element, and so has hybrid tet10: its deviatoric stiffness is blind to the + // four dilatations with a linear volumetric strain, and its constant pressure controls + // one; neighbouring elements restrain the other three (P2–P0, crate::hybrid). Both are + // skipped. + if f != Formulation::Reduced && !(kind == ElementKind::Tet10 && matches!(f, Formulation::Hybrid { .. })) { let mut held = k.clone(); for dof in [0, 1, 2, 4, 5, 8] { held.add(dof, dof, 1e3 * scale); @@ -95,7 +101,9 @@ fn constant_strain_is_reproduced_exactly() { let strain = [a[0][0], a[1][1], a[2][2], a[0][1] + a[1][0], a[0][2] + a[2][0], a[1][2] + a[2][1]]; for d in [isotropic(210e3, 0.3).unwrap(), orthotropic([140e3, 10e3, 12e3], [0.3, 0.25, 0.4], [5e3, 4e3, 3.8e3]).unwrap()] { let want = stress(&d, &strain); - for (kind, f) in CASES { + // Hybrid elements need an isotropic material. + let hybrid_ok = crate::hybrid::split(&d).is_ok(); + for (kind, f) in CASES.into_iter().filter(|(_, f)| hybrid_ok || !matches!(f, Formulation::Hybrid { .. })) { // Curved mid-sides keep the field in the isoparametric span only for straight edges. let xyz = distorted(kind, false); let u: Vec<[f64; 3]> = xyz.iter().map(|x| [0, 1, 2].map(|i| (0..3).map(|j| a[i][j] * x[j]).sum())).collect(); @@ -213,9 +221,12 @@ fn section_formulations() { fn uniform_thermal_strain_is_free_expansion_and_its_load_balances_it() { // A uniform ΔT with orthotropic α: the free expansion u = diag(α ΔT) x is stress-free, and // the thermal load equals K u for it (the element is in equilibrium with no other load). - let d = orthotropic([1000.0, 800.0, 600.0], [0.3, 0.25, 0.2], [400.0, 300.0, 350.0]).unwrap(); + let ortho = orthotropic([1000.0, 800.0, 600.0], [0.3, 0.25, 0.2], [400.0, 300.0, 350.0]).unwrap(); + let iso = isotropic(1000.0, 0.4999).unwrap(); let alpha = [1.2e-5, 2.0e-5, 0.7e-5]; for (kind, f) in CASES { + // Hybrid elements take isotropic (here nearly incompressible) materials. + let d = if matches!(f, Formulation::Hybrid { .. }) { iso } else { ortho }; let xyz = distorted(kind, true); let solid = Solid::new(kind, &xyz, f).unwrap(); let th = crate::ThermalStrain { alpha, delta: vec![80.0; xyz.len()] }; diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 8b5a9b1..0e53ef7 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -17,7 +17,7 @@ use modelcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, ElementDofs, PointS use modelcraft_model::nodal::edge_pressure; use std::collections::BTreeSet; -use modelcraft_model::{ElementKind, Integration, Load, Model, Resolved, SectionKind, Step, Target, centrifugal_acceleration}; +use modelcraft_model::{Elastic, ElementKind, Integration, Load, Model, Resolved, SectionKind, Step, Target, centrifugal_acceleration}; use crate::analysis::MassKind; use crate::{Error, Result}; @@ -398,6 +398,12 @@ impl Mesh { .ok_or_else(|| Error::Invalid(format!("section `{name}`: material `{}` is missing", section.material)))?; let elastic = material.elastic.as_ref().ok_or_else(|| Error::Invalid(format!("material `{}` has no elastic constants", section.material)))?; + if section.formulation.hybrid && !matches!(elastic, Elastic::Isotropic { .. }) { + return Err(Error::NotSupported(format!( + "hybrid elements of the orthotropic material `{}` (section `{name}`): the pressure split needs isotropic elasticity", + section.material + ))); + } let set = resolved.sets.get(§ion.set).ok_or_else(|| Error::Invalid(format!("set `{}` is not resolved", section.set)))?; let bad = |e: modelcraft_fem::Error| match e { modelcraft_fem::Error::NotSupported(m) => Error::NotSupported(format!("section `{name}`: {m}")), diff --git a/crates/solver/tests/hybrid.rs b/crates/solver/tests/hybrid.rs new file mode 100644 index 0000000..b91600b --- /dev/null +++ b/crates/solver/tests/hybrid.rs @@ -0,0 +1,140 @@ +//! Hybrid (mixed u–p) solids for nearly incompressible materials, against closed forms: a block +//! in uniaxial strain (exact for every element) and Lamé's thick cylinder in plane strain at +//! ν = 0.4999, where the displacement elements lock and the hybrid ones do not. Run with +//! `--nocapture` to see the errors. + +mod common; + +use common::*; +use modelcraft_model::*; +use modelcraft_solver::{Limits, solve_static}; + +const NU: f64 = 0.4999; + +fn hybrid() -> Formulation { + Formulation { hybrid: true, ..Formulation::default() } +} + +fn hybrids() -> Vec<(ElementKind, &'static str)> { + vec![(ElementKind::Hex8, "hex8 hybrid"), (ElementKind::Hex20, "hex20 hybrid"), (ElementKind::Tet10, "tet10 hybrid")] +} + +#[test] +fn block_in_uniaxial_strain() { + // Held normal to the bottom and the four sides, pressed on top: ε_zz = −p/M with the + // constrained modulus M = E(1 − ν)/((1 + ν)(1 − 2ν)); σ_xx = σ_yy = −p ν/(1 − ν). + let (e, p) = (10.0, 1.0); + let m = e * (1.0 - NU) / ((1.0 + NU) * (1.0 - 2.0 * NU)); + for (kind, name) in hybrids() { + let mut g = grid([2, 2, 2], kind, |x| x); + g.define(iso(e, NU), None, hybrid(), Procedure::Static {}); + g.symmetry("x0", Symmetry::X, |x| near(x[0], 0.0)); + g.symmetry("x1", Symmetry::X, |x| near(x[0], 1.0)); + g.symmetry("y0", Symmetry::Y, |x| near(x[1], 0.0)); + g.symmetry("y1", Symmetry::Y, |x| near(x[1], 1.0)); + g.symmetry("z0", Symmetry::Z, |x| near(x[2], 0.0)); + g.face_set("top", |x| near(x[2], 1.0)); + g.load("press", Load::Pressure { set: "top".into(), magnitude: p }); + let r = solve_static(&g.model, None, &Limits::default()).unwrap(); + for (i, l) in r.nodes.iter().enumerate() { + let z = g.model.node(*l).unwrap().xyz[2]; + assert!((r.displacement[i][2] + p * z / m).abs() < 1e-9 * p / m, "{name} node {l}: {:?}", r.displacement[i]); + } + for s in &r.stress { + assert!((s[2] + p).abs() < 1e-6 && (s[0] + p * NU / (1.0 - NU)).abs() < 1e-6 && (s[1] - s[0]).abs() < 1e-6, "{name}: {s:?}"); + } + } +} + +/// The inner radial displacement of a quarter of a thick cylinder a = 1, b = 2 in plane strain +/// under internal pressure 1, and Lamé's value. +fn cylinder(kind: ElementKind, formulation: Formulation, n: [usize; 3]) -> (f64, f64) { + let (e, p, a, b) = (1000.0, 1.0, 1.0, 2.0); + let mut g = grid(n, kind, |x| { + let (r, t) = (a + (b - a) * x[0], std::f64::consts::FRAC_PI_2 * x[1]); + [r * t.cos(), r * t.sin(), 0.2 * x[2]] + }); + g.define(iso(e, NU), None, formulation, Procedure::Static {}); + g.symmetry("y0", Symmetry::Y, |x| near(x[1], 0.0)); + g.symmetry("x0", Symmetry::X, |x| near(x[1], 1.0)); + g.symmetry("z0", Symmetry::Z, |x| near(x[2], 0.0)); + g.symmetry("z1", Symmetry::Z, |x| near(x[2], 1.0)); + g.face_set("bore", |x| near(x[0], 0.0)); + g.load("pressure", Load::Pressure { set: "bore".into(), magnitude: p }); + let r = solve_static(&g.model, None, &Limits::default()).unwrap(); + let inner = g.nodes_where(|x| near(x[0], 0.0) && near(x[1], 0.0) && near(x[2], 0.0))[0]; + let u = r.displacement[r.nodes.binary_search(&inner).unwrap()][0]; + let exact = (1.0 + NU) * p * a * a / (e * (b * b - a * a)) * ((1.0 - 2.0 * NU) * a + b * b / a); + (u, exact) +} + +#[test] +fn nearly_incompressible_thick_cylinder_does_not_lock() { + for (kind, name, n, tol) in [ + (ElementKind::Hex8, "hex8 hybrid", [8, 12, 1], 0.01), + (ElementKind::Hex20, "hex20 hybrid", [4, 6, 1], 1e-4), + (ElementKind::Tet10, "tet10 hybrid", [4, 6, 1], 5e-3), + ] { + let (u, exact) = cylinder(kind, hybrid(), n); + eprintln!("{name}: u(a) = {u:.6e}, Lamé {exact:.6e}, error {:.2e}", (u - exact) / exact); + assert!(((u - exact) / exact).abs() < tol, "{name}: {u} vs {exact}"); + } + // The displacement hex8 locks: far too stiff on the same mesh. + let (u, exact) = cylinder(ElementKind::Hex8, Formulation::default(), [8, 12, 1]); + eprintln!("hex8 (displacement): u(a) = {u:.6e}, Lamé {exact:.6e}, error {:.2e}", (u - exact) / exact); + assert!(u < 0.7 * exact, "full hex8 should lock at ν = {NU}: {u} vs {exact}"); +} + +/// Why a one-cell hybrid model of `kind` and `elastic` cannot be solved. +fn refusal(kind: ElementKind, elastic: Elastic) -> String { + let mut g = grid([1, 1, 1], kind, |x| x); + g.define(elastic, None, hybrid(), Procedure::Static {}); + g.fixed("base", |x| near(x[2], 0.0)); + g.face_set("top", |x| near(x[2], 1.0)); + g.load("press", Load::Pressure { set: "top".into(), magnitude: 1.0 }); + solve_static(&g.model, None, &Limits::default()).unwrap_err().to_string() +} + +#[test] +fn hybrid_refusals() { + for (kind, needle) in [(ElementKind::Tet4, "hybrid tet4"), (ElementKind::Wedge6, "hybrid wedge6")] { + let e = refusal(kind, iso(10.0, NU)); + assert!(e.contains(needle), "{e}"); + } + let ortho = Elastic::EngineeringConstants { young: [10.0, 20.0, 30.0], poisson: [0.3, 0.2, 0.1], shear: [5.0, 6.0, 7.0] }; + let e = refusal(ElementKind::Hex8, ortho); + assert!(e.contains("orthotropic material `Mat`"), "{e}"); +} + +#[test] +fn heated_confined_block_carries_the_bulk_thermal_stress() { + // Held normal to all six faces and heated by ΔT: σ = −3 K α ΔT in every direction, the + // pressure the hybrid elements carry as their own unknown. + let (e, alpha, dt) = (10.0, 1e-4, 30.0); + let bulk = e / (3.0 * (1.0 - 2.0 * NU)); + for (kind, name) in hybrids() { + let mut g = grid([2, 2, 2], kind, |x| x); + g.define(iso(e, NU), None, hybrid(), Procedure::Static {}); + let thermal = Thermal { expansion: Some(Directional::Isotropic(alpha)), ..Thermal::default() }; + g.model.replace_material("Mat", Material { elastic: Some(iso(e, NU)), thermal: Some(thermal), ..Material::default() }).unwrap(); + for (set, plane, k, v) in [ + ("x0", Symmetry::X, 0, 0.0), + ("x1", Symmetry::X, 0, 1.0), + ("y0", Symmetry::Y, 1, 0.0), + ("y1", Symmetry::Y, 1, 1.0), + ("z0", Symmetry::Z, 2, 0.0), + ("z1", Symmetry::Z, 2, 1.0), + ] { + g.symmetry(set, plane, |x| near(x[k], v)); + } + g.node_set("everything", |_| true); + g.load("warm", Load::Temperature { target: Target::Set("everything".into()), value: dt }); + let r = solve_static(&g.model, None, &Limits::default()).unwrap(); + let want = -3.0 * bulk * alpha * dt; + for s in &r.stress { + for k in 0..3 { + assert!((s[k] - want).abs() < 1e-6 * want.abs(), "{name}: {s:?} vs {want}"); + } + } + } +} diff --git a/docs/elements.md b/docs/elements.md index 4eee6b6..b5646aa 100644 --- a/docs/elements.md +++ b/docs/elements.md @@ -73,7 +73,7 @@ hex8/20/27 with first and second derivatives). | 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); 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 | +| Hybrid (mixed u–p) solids | Nearly incompressible materials (rubber) | done for hex8 (constant pressure), hex20 (linear pressure) and tet10 (constant pressure), discontinuous pressure condensed, isotropic materials; tet4/wedge6 refused (they lock); continuous-pressure (Taylor–Hood) elements planned with an indefinite factorization | | Modified quadratic tet (tet10 for contact/plasticity) | Robust quadratic tets | planned | ## Special-purpose and multiphysics diff --git a/docs/fem.md b/docs/fem.md index 175925f..02aa909 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -133,3 +133,19 @@ or per axis, and `reference_temperature`), `BoundaryCondition::Temperature`, loa in a static step), procedures `HeatTransfer { transient }` and `CoupledTemperatureDisplacement`. Commands: `solve.heat`, and `solve.static` for coupled steps. The deck writers refuse heat transfer by name for now. +## Hybrid (mixed u–p) solids + +A section with `"formulation": {"hybrid": true}` gives `Formulation::Hybrid { pressures }` +(`crate::hybrid`): the isotropic elasticity splits as D = D_dev + K m mᵀ, the mean stress is an +independent field discontinuous between elements (hex20: 1, x, y, z; hex8 and tet10: constant), +and it is condensed out element by element, p̂ = K M⁻¹ Gᵀ u, leaving the stiffness +K_dev + K G M⁻¹ Gᵀ. `ElementDofs.internal.pressure` counts the condensed unknowns; nothing +global changes, and Cholesky still applies. Stresses are D_dev ε + p m with the element's +pressure field. The choice of pressure space follows the inf-sup analysis (docs in +`crates/fem/src/hybrid.rs`): Q2–P1disc and P2–P0 are stable; Q1–P0 does not lock but may show +checkerboard pressures; linear tetrahedra and wedges lock even so and are refused. A thermal +strain splits the same way: ∫ Bᵀ D_dev ε_th on the nodes and, through the pressure, +K G M⁻¹ ∫ N_p tr ε_th, so a heated nearly incompressible solid takes its bulk stress −3KαΔT +when confined. A continuous +nodal pressure (`DofKind::P`, Taylor–Hood) needs a symmetric indefinite factorization and is +planned. diff --git a/docs/solvers/literature.md b/docs/solvers/literature.md index d75ddd9..39865bf 100644 --- a/docs/solvers/literature.md +++ b/docs/solvers/literature.md @@ -1141,3 +1141,17 @@ update are written from (cited again in the module docs): constrained thermal expansion (σ = −EαΔT), the bimetal strip (S. Timoshenko, *Analysis of bi-metal thermostats*, J. Opt. Soc. Am. 11(3):233–255, 1925), and the coupled step against the heat step followed by the static step. + +## Part L: Hybrid (mixed u–p) solids, as implemented + +- Perturbed-Lagrangian u–p formulation with element-wise condensed discontinuous pressure: + T. J. R. Hughes, *The Finite Element Method*, Dover 2000, §4.3–4.5 (B-bar and mean + dilatation as the constant-pressure case); J. C. Nagtegaal, D. M. Parks, J. R. Rice, CMAME + 4(2):153–177, 1974 (locking of displacement elements near incompressibility). +- Pressure spaces and the inf-sup condition: F. Brezzi, M. Fortin, *Mixed and Hybrid Finite + Element Methods*, Springer 1991, ch. VI; K.-J. Bathe, *Finite Element Procedures*, 2nd ed., + 2014, §4.5 (Q2–P1disc and P2–P0 stable, Q1–P0 checkerboard-prone, P1–P0 locks). +- Verification (`crates/solver/tests/hybrid.rs`, analytic): a block in uniaxial strain at + ν = 0.4999 (exact), Lamé's thick cylinder in plane strain at ν = 0.4999 (hex20 1e-5, + tet10 1e-3, hex8 3e-3 relative error in the bore displacement, where the displacement hex8 + is 80 % too stiff). Cook's membrane is a published benchmark, not part of the analytic gate.