diff --git a/crates/fem/src/constitutive.rs b/crates/fem/src/constitutive.rs new file mode 100644 index 0000000..d801387 --- /dev/null +++ b/crates/fem/src/constitutive.rs @@ -0,0 +1,356 @@ +//! The material point of finite-strain analysis: from the deformation gradient F (and the +//! material's history) to the stress and its consistent tangent. +//! +//! **Stress measure.** A material returns the **Kirchhoff stress** τ = J σ (J = det F, σ the +//! Cauchy stress) and the **spatial tangent** c of the Lie derivative of τ, L_v τ = c : d, with +//! d the rate of deformation. Both are 6-vectors and 6 × 6 matrices in the order xx, yy, zz, +//! xy, xz, yz, c taking engineering shears like the small-strain elasticity matrices (so the +//! linear elasticity D is c at F = I). This is the total Lagrangian formulation written with +//! spatial quantities: τ = F S Fᵀ and c_ijkl = F_iI F_jJ F_kK F_lL C_IJKL push the second +//! Piola–Kirchhoff stress S and its material tangent C = ∂S/∂E forward, so the internal forces +//! ∫ Bₓᵀ τ dV₀ and the stiffness ∫ (Bₓᵀ c Bₓ + ∇ₓN·τ∇ₓN) dV₀, with the shape function gradients +//! of the current configuration and integrals over the reference volume, are exactly the total +//! Lagrangian ∫ B₀ᵀ S dV₀ and its linearization (J. Bonet, R. D. Wood, *Nonlinear Continuum +//! Mechanics for Finite Element Analysis*, 2nd ed., Cambridge 2008, §7.4 and §9.2; G. A. +//! Holzapfel, *Nonlinear Solid Mechanics*, Wiley 2000, §6.6). Finite-strain plasticity (the +//! exponential return map of J. C. Simo, CMAME 99:61–112, 1992) gives τ and c directly, which is +//! why the interface is spatial. +//! +//! **History.** A material keeps [`MaterialModel::history_len`] numbers per integration point. +//! An evaluation reads the history committed at the end of the last converged increment and +//! writes the updated one to a separate slice; the solver commits it when the increment +//! converges, so a rejected iteration or a cut-back increment leaves no trace. The elastic +//! materials here keep none. +//! +//! **Materials:** +//! - **St Venant–Kirchhoff**: S = D E with E = ½(FᵀF − I) the Green–Lagrange strain and D the +//! small-strain elasticity (isotropic or orthotropic in the global axes): an elastic material +//! in a nonlinear step. Exact for large rotations with small strains; it softens without bound +//! in large compression, so large strains want a hyperelastic model. +//! - **Neo-Hookean** (compressible, with the isochoric–volumetric split of P. J. Flory, +//! *Thermodynamic relations for high elastic materials*, Trans. Faraday Soc. 57:829–838, 1961): +//! W = μ/2 (Ī₁ − 3) + U(J), Ī₁ = J^(−2/3) tr b, U = κ/2 (J − 1)². Then +//! τ = μ dev b̄ + J U′ I with b̄ = J^(−2/3) b, and +//! c = ⅔ μ tr b̄ ℙ − ⅔ (τ_iso ⊗ I + I ⊗ τ_iso) + J (p̃ I ⊗ I − 2p 𝕀), with p = U′(J), +//! p̃ = p + J U″(J), 𝕀 the symmetric fourth-order identity and ℙ = 𝕀 − ⅓ I ⊗ I (Holzapfel 2000, +//! §6.6, eqs. 6.165–6.169 and 6.193, in Kirchhoff form). The split is what the mixed (hybrid) +//! elements need: they take the isochoric part from [`MaterialModel::isochoric`] and the +//! volumetric one from [`MaterialModel::volumetric`] with the element's own dilatation. + +use fieldcraft_model::{Elastic, Hyperelastic}; + +use crate::dense::det3; +use crate::material::{Elasticity, elasticity, stress}; +use crate::rotation::{IDENTITY, Mat3, mul, transpose, voigt}; +use crate::{Error, Result}; + +/// A 6 × 6 tangent in Voigt order, engineering-shear columns. +pub type Tangent = [[f64; 6]; 6]; + +/// The (row, column) of each Voigt component. +pub const PAIRS: [(usize, usize); 6] = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]; + +/// The Kirchhoff stress and the spatial tangent of a material point. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Response { + pub stress: [f64; 6], + pub tangent: Tangent, +} + +/// The constitutive model of a solid in a nonlinear step. +#[derive(Debug, Clone, PartialEq)] +pub enum MaterialModel { + StVenantKirchhoff(Box), + NeoHookean { shear: f64, bulk: f64 }, +} + +/// The 6 × 6 matrix of the symmetric fourth-order identity 𝕀 (engineering-shear columns). +fn identity4() -> Tangent { + let mut i4 = [[0.0; 6]; 6]; + for k in 0..6 { + i4[k][k] = if k < 3 { 1.0 } else { 0.5 }; + } + i4 +} + +/// The matrix T with (F A Fᵀ) = T a for a symmetric A of Voigt vector a (tensor components): +/// T[(ij)][(IJ)] = F_iI F_jJ + F_iJ F_jI for I ≠ J and F_iI F_jI for I = J. A tangent pushes +/// forward as T C Tᵀ. +pub fn push_forward(f: &Mat3) -> Tangent { + let mut t = [[0.0; 6]; 6]; + for (a, &(i, j)) in PAIRS.iter().enumerate() { + for (b, &(ii, jj)) in PAIRS.iter().enumerate() { + t[a][b] = if ii == jj { f[i][ii] * f[j][ii] } else { f[i][ii] * f[j][jj] + f[i][jj] * f[j][ii] }; + } + } + t +} + +/// T C Tᵀ. +fn sandwich(t: &Tangent, c: &Tangent) -> Tangent { + let mut tc = [[0.0; 6]; 6]; + for a in 0..6 { + for b in 0..6 { + tc[a][b] = (0..6).map(|k| t[a][k] * c[k][b]).sum(); + } + } + let mut out = [[0.0; 6]; 6]; + for a in 0..6 { + for b in 0..6 { + out[a][b] = (0..6).map(|k| tc[a][k] * t[b][k]).sum(); + } + } + out +} + +/// det F, refusing an inverted or collapsed point. +fn jacobian(f: &Mat3) -> Result { + let j = det3(f); + if j > 0.0 && j.is_finite() { Ok(j) } else { Err(Error::Degenerate(format!("the material is inverted (J = det F = {j:.3e} is not positive)"))) } +} + +impl MaterialModel { + /// The model of a material's elastic behaviour in a nonlinear step: neo-Hookean when the + /// material asks for it (its isotropic constants give μ and κ), St Venant–Kirchhoff + /// otherwise. + pub fn from_model(elastic: &Elastic, hyperelastic: Option<&Hyperelastic>) -> Result { + match hyperelastic { + None => Ok(MaterialModel::StVenantKirchhoff(Box::new(elasticity(elastic)?))), + Some(Hyperelastic::NeoHookean {}) => { + let (young, poisson) = crate::material::isotropic_constants(elastic, "neo-Hookean")?; + let (shear, bulk) = Hyperelastic::moduli(young, poisson); + Ok(MaterialModel::NeoHookean { shear, bulk }) + } + } + } + + /// The model of an element whose elasticity is `d` (the small-strain one) and whose material + /// has the large-strain option `hyperelastic`: neo-Hookean takes μ = D₄₄ and κ = D₁₂ + ⅔ μ of + /// an isotropic `d`, St Venant–Kirchhoff `d` itself. + pub fn from_elasticity(d: &Elasticity, hyperelastic: Option<&Hyperelastic>) -> MaterialModel { + match hyperelastic { + None => MaterialModel::StVenantKirchhoff(Box::new(*d)), + Some(Hyperelastic::NeoHookean {}) => MaterialModel::NeoHookean { shear: d[3][3], bulk: d[0][1] + 2.0 / 3.0 * d[3][3] }, + } + } + + pub fn name(&self) -> &'static str { + match self { + MaterialModel::StVenantKirchhoff(_) => "St Venant–Kirchhoff", + MaterialModel::NeoHookean { .. } => "neo-Hookean", + } + } + + /// History variables per integration point. + pub fn history_len(&self) -> usize { + 0 + } + + /// Whether the strain energy splits into isochoric and volumetric parts (what mixed + /// elements need). + pub fn is_decoupled(&self) -> bool { + matches!(self, MaterialModel::NeoHookean { .. }) + } + + /// The Kirchhoff stress and spatial tangent at deformation gradient `f`, from the history + /// `committed` (written, updated, to `updated`). + pub fn response(&self, f: &Mat3, committed: &[f64], updated: &mut [f64]) -> Result { + updated.iter_mut().zip(committed).for_each(|(u, c)| *u = *c); + let j = jacobian(f)?; + match self { + MaterialModel::StVenantKirchhoff(d) => { + let c = mul(&transpose(f), f); + let e = [0.5 * (c[0][0] - 1.0), 0.5 * (c[1][1] - 1.0), 0.5 * (c[2][2] - 1.0), c[0][1], c[0][2], c[1][2]]; + let s = stress(d, &e); + let t = push_forward(f); + Ok(Response { stress: std::array::from_fn(|a| (0..6).map(|b| t[a][b] * s[b]).sum()), tangent: sandwich(&t, d) }) + } + MaterialModel::NeoHookean { .. } => { + let iso = self.isochoric(f)?; + let (p, dp) = self.volumetric(j).ok_or_else(|| Error::Invalid("volumetric response".into()))?; + Ok(add_pressure(&iso, j, p, j * dp)) + } + } + } + + /// The isochoric part of a decoupled material (τ_iso = dev τ̄ and its tangent). + pub fn isochoric(&self, f: &Mat3) -> Result { + let &MaterialModel::NeoHookean { shear: mu, .. } = self else { + return Err(Error::NotSupported(format!("the isochoric–volumetric split of a {} material", self.name()))); + }; + let j = jacobian(f)?; + let scale = j.powf(-2.0 / 3.0); + let b = mul(f, &transpose(f)); + let bbar = voigt(&b).map(|v| scale * v); + let tr = bbar[0] + bbar[1] + bbar[2]; + let mut tau = bbar.map(|v| mu * v); + for v in tau.iter_mut().take(3) { + *v -= mu * tr / 3.0; + } + let i4 = identity4(); + let mut c = [[0.0; 6]; 6]; + for a in 0..6 { + for b in 0..6 { + let (ia, ib) = (if a < 3 { 1.0 } else { 0.0 }, if b < 3 { 1.0 } else { 0.0 }); + let projector = i4[a][b] - ia * ib / 3.0; + c[a][b] = 2.0 / 3.0 * mu * tr * projector - 2.0 / 3.0 * (tau[a] * ib + ia * tau[b]); + } + } + Ok(Response { stress: tau, tangent: c }) + } + + /// U′(J) and U″(J) of a decoupled material's volumetric energy. + pub fn volumetric(&self, j: f64) -> Option<(f64, f64)> { + match *self { + MaterialModel::NeoHookean { bulk, .. } => Some((bulk * (j - 1.0), bulk)), + MaterialModel::StVenantKirchhoff(_) => None, + } + } +} + +/// An isochoric response with a pressure p acting on dilatation J: τ + J p I, and the tangent +/// c + J (p̃ I ⊗ I − 2p 𝕀) with p̃ = p + `j_dp` (J dp/dJ, zero when p is held). +pub fn add_pressure(iso: &Response, j: f64, p: f64, j_dp: f64) -> Response { + let mut out = *iso; + let i4 = identity4(); + for a in 0..3 { + out.stress[a] += j * p; + } + for a in 0..6 { + for b in 0..6 { + let ii = if a < 3 && b < 3 { 1.0 } else { 0.0 }; + out.tangent[a][b] += j * ((p + j_dp) * ii - 2.0 * p * i4[a][b]); + } + } + out +} + +/// The logarithmic (Hencky) strain ln V = ½ ln b of a deformation gradient, as a strain 6-vector +/// with engineering shears; `None` when F is not finite or inverted. +pub fn logarithmic_strain(f: &Mat3) -> Option<[f64; 6]> { + let (l, v) = crate::dense::symmetric_eigen(&mul(f, &transpose(f)))?; + if !l.iter().all(|x| *x > 0.0) { + return None; + } + let ln = l.map(|x| 0.5 * x.ln()); + let e = |i: usize, j: usize| (0..3).map(|k| ln[k] * v[i][k] * v[j][k]).sum::(); + Some([e(0, 0), e(1, 1), e(2, 2), 2.0 * e(0, 1), 2.0 * e(0, 2), 2.0 * e(1, 2)]) +} + +/// The deformation gradient F = I + Σ u_a ⊗ ∇₀N_a of nodal displacements and reference +/// gradients. +pub fn deformation_gradient(u: &[[f64; 3]], gradients: &[[f64; 3]]) -> Mat3 { + let mut f = IDENTITY; + for (ua, g) in u.iter().zip(gradients) { + for i in 0..3 { + for jj in 0..3 { + f[i][jj] += ua[i] * g[jj]; + } + } + } + f +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::dense::inverse3; + use crate::material::isotropic; + + fn deformation() -> Mat3 { + [[1.2, 0.15, -0.05], [-0.1, 0.9, 0.2], [0.07, -0.12, 1.05]] + } + + /// S = F⁻¹ τ F⁻ᵀ. + fn pull_back_stress(f: &Mat3, tau: &[f64; 6]) -> [f64; 6] { + crate::rotation::rotate(&inverse3(f).unwrap().0, tau) + } + + /// The tangent's consistency: dS = C : dE with C the pull-back of c, checked by central + /// differences of S along several directions dF (dE = sym(Fᵀ dF), engineering shears). + fn check_tangent(m: &MaterialModel, f: &Mat3) { + let r = m.response(f, &[], &mut []).unwrap(); + let (inv, _) = inverse3(f).unwrap(); + let back = sandwich(&push_forward(&inv), &r.tangent); + let scale = back.iter().flatten().fold(0.0f64, |a, v| a.max(v.abs())); + let h = 1e-6; + for k in 0..9 { + let mut df = [[0.0; 3]; 3]; + df[k / 3][k % 3] = 1.0; + df[(k + 1) % 3][(k + 2) % 3] += 0.3; + let at = |s: f64| { + let fp: Mat3 = std::array::from_fn(|i| std::array::from_fn(|j| f[i][j] + s * df[i][j])); + pull_back_stress(&fp, &m.response(&fp, &[], &mut []).unwrap().stress) + }; + let (sp, sm) = (at(h), at(-h)); + let ftdf = mul(&transpose(f), &df); + let de = [ftdf[0][0], ftdf[1][1], ftdf[2][2], ftdf[0][1] + ftdf[1][0], ftdf[0][2] + ftdf[2][0], ftdf[1][2] + ftdf[2][1]]; + for a in 0..6 { + let numeric = (sp[a] - sm[a]) / (2.0 * h); + let analytic: f64 = (0..6).map(|b| back[a][b] * de[b]).sum(); + assert!((numeric - analytic).abs() < 1e-6 * scale, "{} {k} {a}: {numeric} vs {analytic}", m.name()); + } + } + } + + #[test] + fn tangents_are_consistent_with_the_stress() { + let d = isotropic(1000.0, 0.3).unwrap(); + for m in [MaterialModel::StVenantKirchhoff(Box::new(d)), MaterialModel::NeoHookean { shear: 300.0, bulk: 5000.0 }] { + check_tangent(&m, &deformation()); + check_tangent(&m, &IDENTITY); + // Symmetric (hyperelastic): major symmetry of c. + let c = m.response(&deformation(), &[], &mut []).unwrap().tangent; + for a in 0..6 { + for b in 0..6 { + assert!((c[a][b] - c[b][a]).abs() < 1e-9 * 5000.0, "{a} {b}"); + } + } + } + let ortho = crate::material::orthotropic([1000.0, 300.0, 500.0], [0.3, 0.2, 0.25], [100.0, 120.0, 90.0]).unwrap(); + check_tangent(&MaterialModel::StVenantKirchhoff(Box::new(ortho)), &deformation()); + } + + #[test] + fn small_strain_limit_is_linear_elasticity_and_rotations_do_not_stress() { + let (e, nu) = (1000.0, 0.3); + let d = isotropic(e, nu).unwrap(); + let neo = MaterialModel::from_model(&Elastic::Isotropic { young: e, poisson: nu }, Some(&Hyperelastic::NeoHookean {})).unwrap(); + for m in [MaterialModel::StVenantKirchhoff(Box::new(d)), neo] { + let c = m.response(&IDENTITY, &[], &mut []).unwrap().tangent; + for a in 0..6 { + for b in 0..6 { + assert!((c[a][b] - d[a][b]).abs() < 1e-9 * e, "{}: {a} {b}", m.name()); + } + } + // A rigid rotation (90° about an inclined axis) leaves the stress at zero. + let (s, cth) = (1.0f64, 0.0f64); + let axis = [1.0 / 3f64.sqrt(); 3]; + let k = crate::rotation::skew(&axis); + let k2 = mul(&k, &k); + let rot: Mat3 = std::array::from_fn(|i| std::array::from_fn(|j| IDENTITY[i][j] + s * k[i][j] + (1.0 - cth) * k2[i][j])); + let tau = m.response(&rot, &[], &mut []).unwrap().stress; + assert!(tau.iter().all(|v| v.abs() < 1e-12 * e), "{}: {tau:?}", m.name()); + } + } + + #[test] + fn uniaxial_stretch_of_neo_hookean_and_log_strain() { + // Stretch λ along x with the lateral stretch μ chosen so that σ_yy = 0: check σ_xx + // against W directly, σ = (μ/J^(5/3)) (λ² − (λ² + 2 l²)/3) + κ (J − 1). + let (mu, kappa) = (300.0, 5000.0); + let m = MaterialModel::NeoHookean { shear: mu, bulk: kappa }; + let (lam, l) = (1.5f64, 0.85f64); + let f = [[lam, 0.0, 0.0], [0.0, l, 0.0], [0.0, 0.0, l]]; + let j = lam * l * l; + let tau = m.response(&f, &[], &mut []).unwrap().stress; + let i1 = lam * lam + 2.0 * l * l; + let sxx = mu / j.powf(5.0 / 3.0) * (lam * lam - i1 / 3.0) + kappa * (j - 1.0); + let syy = mu / j.powf(5.0 / 3.0) * (l * l - i1 / 3.0) + kappa * (j - 1.0); + assert!((tau[0] / j - sxx).abs() < 1e-10 * mu && (tau[1] / j - syy).abs() < 1e-10 * mu, "{tau:?}"); + let e = logarithmic_strain(&f).unwrap(); + assert!((e[0] - lam.ln()).abs() < 1e-14 && (e[1] - l.ln()).abs() < 1e-14 && e[3].abs() < 1e-14, "{e:?}"); + assert!(m.response(&[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, -1.0]], &[], &mut []).is_err(), "inverted"); + assert!(MaterialModel::StVenantKirchhoff(Box::new(isotropic(1.0, 0.3).unwrap())).isochoric(&IDENTITY).is_err()); + } +} diff --git a/crates/fem/src/dense.rs b/crates/fem/src/dense.rs index 294aa6b..4d4137e 100644 --- a/crates/fem/src/dense.rs +++ b/crates/fem/src/dense.rs @@ -113,6 +113,57 @@ pub fn positive_inverse(j: &[[f64; 3]; 3]) -> Option<([[f64; 3]; 3], f64)> { (det > 1e-12 * norms).then_some((inv, det)) } +/// The eigenvalues and unit eigenvectors (the columns of the matrix returned) of a symmetric +/// 3 × 3 matrix, by cyclic Jacobi rotations (G. H. Golub, C. F. Van Loan, *Matrix +/// Computations*, 4th ed., Johns Hopkins 2013, §8.5): each rotation zeroes one off-diagonal +/// entry with the stable choice of angle, and the sweeps stop when the off-diagonal part is +/// round-off relative to the diagonal (at most 32 sweeps; three or four suffice). Accurate to a +/// few ulps of the largest eigenvalue; `None` for non-finite input. +pub fn symmetric_eigen(a: &[[f64; 3]; 3]) -> Option<([f64; 3], [[f64; 3]; 3])> { + if a.iter().flatten().any(|v| !v.is_finite()) { + return None; + } + let mut m = [[a[0][0], 0.5 * (a[0][1] + a[1][0]), 0.5 * (a[0][2] + a[2][0])], [0.0; 3], [0.0; 3]]; + m[1] = [m[0][1], a[1][1], 0.5 * (a[1][2] + a[2][1])]; + m[2] = [m[0][2], m[1][2], a[2][2]]; + let mut v = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]; + for _ in 0..32 { + let off = m[0][1] * m[0][1] + m[0][2] * m[0][2] + m[1][2] * m[1][2]; + let diag = m[0][0] * m[0][0] + m[1][1] * m[1][1] + m[2][2] * m[2][2]; + if off <= 1e-36 * diag || off == 0.0 { + break; + } + for (p, q) in [(0, 1), (0, 2), (1, 2)] { + let apq = m[p][q]; + if apq == 0.0 { + continue; + } + // tan of the rotation angle, the smaller root of t² + 2θt − 1 = 0. + let theta = (m[q][q] - m[p][p]) / (2.0 * apq); + let sign = if theta >= 0.0 { 1.0 } else { -1.0 }; + let t = sign / (theta.abs() + (theta * theta + 1.0).sqrt()); + let c = 1.0 / (t * t + 1.0).sqrt(); + let s = t * c; + for k in 0..3 { + let (mkp, mkq) = (m[k][p], m[k][q]); + m[k][p] = c * mkp - s * mkq; + m[k][q] = s * mkp + c * mkq; + } + for k in 0..3 { + let (mpk, mqk) = (m[p][k], m[q][k]); + m[p][k] = c * mpk - s * mqk; + m[q][k] = s * mpk + c * mqk; + } + for row in &mut v { + let (vp, vq) = (row[p], row[q]); + row[p] = c * vp - s * vq; + row[q] = s * vp + c * vq; + } + } + } + Some(([m[0][0], m[1][1], m[2][2]], v)) +} + 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]) @@ -143,4 +194,29 @@ mod tests { } assert!(inverse3(&[[1.0, 2.0, 3.0], [2.0, 4.0, 6.0], [0.0, 0.0, 1.0]]).is_none()); } + + #[test] + fn symmetric_eigen_decomposes() { + for a in [ + [[4.0, 1.0, -2.0], [1.0, 3.0, 0.5], [-2.0, 0.5, 1.0]], + [[2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 5.0]], + [[1.0, 1e-9, 0.0], [1e-9, 1.0, 0.0], [0.0, 0.0, 1.0]], + [[1e6, 3.0, 2.0], [3.0, 1e-3, 0.0], [2.0, 0.0, -7.0]], + ] { + let (l, v) = symmetric_eigen(&a).unwrap(); + let scale = l.iter().fold(0.0f64, |m, x| m.max(x.abs())); + for k in 0..3 { + // A v_k = λ_k v_k and the v_k are orthonormal. + for i in 0..3 { + let av: f64 = (0..3).map(|j| a[i][j] * v[j][k]).sum(); + assert!((av - l[k] * v[i][k]).abs() < 1e-12 * scale, "{a:?}: {l:?}"); + } + for m in 0..3 { + let d: f64 = (0..3).map(|i| v[i][k] * v[i][m]).sum(); + assert!((d - if k == m { 1.0 } else { 0.0 }).abs() < 1e-13); + } + } + } + assert!(symmetric_eigen(&[[f64::NAN, 0.0, 0.0], [0.0; 3], [0.0; 3]]).is_none()); + } } diff --git a/crates/fem/src/face.rs b/crates/fem/src/face.rs index 7bf6684..4403f15 100644 --- a/crates/fem/src/face.rs +++ b/crates/fem/src/face.rs @@ -21,6 +21,113 @@ pub fn pressure(xyz: &[[f64; 3]], magnitude: f64) -> Result> { Ok(normal_weights(xyz)?.into_iter().map(|w| w.map(|c| -magnitude * c)).collect()) } +/// A pressure that follows the deformed face (positive pushes in): the nodal forces +/// fᵢ = −p ∫ Nᵢ (x,ξ × x,η) dξ dη on the current positions x = X + u (reference coordinates `xyz`, +/// displacements `u`, in the face's node order), and the load stiffness K = −∂f/∂u, 3n × 3n +/// row-major, from Δ(x,ξ × x,η) = [x,ξ]× Δx,η − [x,η]× Δx,ξ: +/// K_ab = p ∫ N_a (N_b,η [x,ξ]× − N_b,ξ [x,η]×) dξ dη (K. Schweizerhof, E. Ramm, *Displacement +/// dependent pressure loads in nonlinear finite element analyses*, Comput. Struct. +/// 18(6):1099–1114, 1984). K is not symmetric on its own; over a closed surface, or one whose +/// boundary is held, the assembled skew parts cancel. Same rule as [`pressure`], so at u = 0 the +/// forces are its forces. +pub fn follower_pressure(xyz: &[[f64; 3]], u: &[[f64; 3]], magnitude: f64) -> Result<(Vec<[f64; 3]>, Vec)> { + let n = xyz.len(); + if u.len() != n { + return Err(Error::Invalid(format!("{} displacements for a {n}-node face", u.len()))); + } + let x: Vec<[f64; 3]> = xyz.iter().zip(u).map(|(a, b)| [a[0] + b[0], a[1] + b[1], a[2] + b[2]]).collect(); + if x.iter().flatten().any(|c| !c.is_finite()) || !magnitude.is_finite() { + return Err(Error::Invalid("face positions and pressure must be finite".into())); + } + let points = fieldcraft_model::nodal::face_points(n).map_err(|e| Error::Invalid(e.to_string()))?; + let mut f = vec![[0.0; 3]; n]; + let mut k = vec![0.0; 9 * n * n]; + for q in &points { + let (mut tx, mut te) = ([0.0; 3], [0.0; 3]); + for ((p, dx), de) in x.iter().zip(&q.dxi).zip(&q.deta) { + for c in 0..3 { + tx[c] += dx * p[c]; + te[c] += de * p[c]; + } + } + let normal = [tx[1] * te[2] - tx[2] * te[1], tx[2] * te[0] - tx[0] * te[2], tx[0] * te[1] - tx[1] * te[0]]; + let (sx, se) = (crate::rotation::skew(&tx), crate::rotation::skew(&te)); + for (a, na) in q.n.iter().enumerate() { + if let Some(fa) = f.get_mut(a) { + for c in 0..3 { + fa[c] -= magnitude * q.weight * na * normal[c]; + } + } + for (b, (db_xi, db_eta)) in q.dxi.iter().zip(&q.deta).enumerate() { + for r in 0..3 { + for c in 0..3 { + let v = magnitude * q.weight * na * (db_eta * sx[r][c] - db_xi * se[r][c]); + if let Some(e) = k.get_mut((3 * a + r) * 3 * n + 3 * b + c) { + *e += v; + } + } + } + } + } + } + Ok((f, k)) +} + +/// A pressure that follows a deformed edge of a 2D element (positive pushes into a +/// counter-clockwise element): the nodal forces (x, y) fᵢ = −p ∫ Nᵢ e (t_y, −t_x) ds on the +/// current positions, with t the tangent and e the out-of-plane extent (a thickness, or 2π r at +/// the current radius when revolved), and the load stiffness K = −∂f/∂u, 2n × 2n row-major (the +/// revolved extent depends on the radial motion too). +pub fn follower_edge_pressure( + xyz: &[[f64; 3]], + u: &[[f64; 2]], + extent: fieldcraft_model::Extent, + magnitude: f64, +) -> Result<(Vec<[f64; 2]>, Vec)> { + let n = xyz.len(); + if u.len() != n { + return Err(Error::Invalid(format!("{} displacements for a {n}-node edge", u.len()))); + } + let x: Vec<[f64; 2]> = xyz.iter().zip(u).map(|(a, b)| [a[0] + b[0], a[1] + b[1]]).collect(); + if x.iter().flatten().any(|c| !c.is_finite()) || !magnitude.is_finite() { + return Err(Error::Invalid("edge positions and pressure must be finite".into())); + } + let revolved = matches!(extent, fieldcraft_model::Extent::Revolved); + let points = fieldcraft_model::nodal::edge_points(n).map_err(|e| Error::Invalid(e.to_string()))?; + let mut f = vec![[0.0; 2]; n]; + let mut k = vec![0.0; 4 * n * n]; + let two_pi = 2.0 * std::f64::consts::PI; + for q in &points { + let (mut t, mut pos) = ([0.0; 2], [0.0; 2]); + for ((p, dn), ni) in x.iter().zip(&q.dxi).zip(&q.n) { + for c in 0..2 { + t[c] += dn * p[c]; + pos[c] += ni * p[c]; + } + } + let e = extent.at([pos[0], pos[1], 0.0]); + let w = magnitude * q.weight; + for (a, na) in q.n.iter().enumerate() { + if let Some(fa) = f.get_mut(a) { + fa[0] -= w * na * e * t[1]; + fa[1] += w * na * e * t[0]; + } + for (b, (nb, db)) in q.n.iter().zip(&q.dxi).enumerate() { + // ∂e/∂u_b (radial) = 2π N_b when revolved. + let de = if revolved { two_pi * nb } else { 0.0 }; + // K = −∂f/∂u. + let entries = [(0, 0, w * na * de * t[1]), (0, 1, w * na * e * db), (1, 0, -w * na * (e * db + de * t[0])), (1, 1, 0.0)]; + for (r, c, v) in entries { + if let Some(x) = k.get_mut((2 * a + r) * 2 * n + 2 * b + c) { + *x += v; + } + } + } + } + } + Ok((f, k)) +} + #[cfg(test)] mod tests { use super::*; @@ -57,6 +164,51 @@ mod tests { assert!(pressure(&tri[..2], 1.0).is_err()); } + /// Central differences of the follower forces against the load stiffness, on a distorted + /// quadratic face and a revolved quadratic edge; at u = 0 the forces are the dead ones. + #[test] + fn follower_pressure_stiffness_is_the_derivative_of_its_forces() { + let face = [[0., 0., 0.], [2., 0., 0.1], [2.2, 3., 0.], [0., 3., -0.2], [1., -0.1, 0.], [2.1, 1.5, 0.], [1., 3.1, 0.1], [0.1, 1.5, 0.]]; + let u: Vec<[f64; 3]> = (0..8).map(|i| [0.1 * (i as f64).sin(), 0.05 * (i as f64).cos(), 0.2 * ((i * i) as f64).sin()]).collect(); + let (f0, _) = follower_pressure(&face, &[[0.0; 3]; 8], 3.0).unwrap(); + let dead = pressure(&face, 3.0).unwrap(); + assert!(f0.iter().flatten().zip(dead.iter().flatten()).all(|(a, b)| (a - b).abs() < 1e-13), "{f0:?}"); + let (_, k) = follower_pressure(&face, &u, 3.0).unwrap(); + let h = 1e-6; + for j in 0..24 { + let shifted = |s: f64| { + let mut v = u.clone(); + v[j / 3][j % 3] += s; + follower_pressure(&face, &v, 3.0).unwrap().0 + }; + let (p, m) = (shifted(h), shifted(-h)); + for i in 0..24 { + let numeric = -(p[i / 3][i % 3] - m[i / 3][i % 3]) / (2.0 * h); + assert!((numeric - k[i * 24 + j]).abs() < 1e-7, "{i} {j}: {numeric} vs {}", k[i * 24 + j]); + } + } + let edge = [[1.0, 0.0, 0.0], [1.5, 2.0, 0.0], [1.3, 1.0, 0.0]]; + let ue = [[0.1, 0.05], [-0.2, 0.1], [0.05, -0.03]]; + for extent in [fieldcraft_model::Extent::Revolved, fieldcraft_model::Extent::Thickness(2.0)] { + let (f0, _) = follower_edge_pressure(&edge, &[[0.0; 2]; 3], extent, 2.0).unwrap(); + let dead = fieldcraft_model::nodal::edge_pressure(&edge, extent, 2.0).unwrap(); + assert!(f0.iter().zip(&dead).all(|(a, b)| (a[0] - b[0]).abs() < 1e-12 && (a[1] - b[1]).abs() < 1e-12)); + let (_, k) = follower_edge_pressure(&edge, &ue, extent, 2.0).unwrap(); + for j in 0..6 { + let shifted = |s: f64| { + let mut v = ue; + v[j / 2][j % 2] += s; + follower_edge_pressure(&edge, &v, extent, 2.0).unwrap().0 + }; + let (p, m) = (shifted(h), shifted(-h)); + for i in 0..6 { + let numeric = -(p[i / 2][i % 2] - m[i / 2][i % 2]) / (2.0 * h); + assert!((numeric - k[i * 6 + j]).abs() < 1e-7, "{extent:?} {i} {j}: {numeric} vs {}", k[i * 6 + j]); + } + } + } + } + #[test] fn pressure_on_a_curved_face_sums_to_the_projected_area() { // A quadratic face on a quarter cylinder r = 2, z in 0..1, normal outward: the total diff --git a/crates/fem/src/finite.rs b/crates/fem/src/finite.rs new file mode 100644 index 0000000..fad60c1 --- /dev/null +++ b/crates/fem/src/finite.rs @@ -0,0 +1,561 @@ +//! Finite-strain (total Lagrangian) solid and 2D elements: internal forces and the consistent +//! tangent stiffness for nodal displacements u, and the Cauchy stress and logarithmic strain at +//! the integration points. +//! +//! **Formulation.** Total Lagrangian: the integrals run over the reference (undeformed) element +//! with its own quadrature, F = I + Σ u_a ⊗ ∇₀N_a at each point, and the material returns the +//! Kirchhoff stress τ and its spatial tangent c ([`crate::constitutive`]). With the shape +//! function gradients of the current configuration ∇ₓN_a = F⁻ᵀ ∇₀N_a, +//! - internal forces f_a = ∫ Bₓ,aᵀ τ dV₀, +//! - tangent K_ab = ∫ (Bₓ,aᵀ c Bₓ,b + (∇ₓN_a · τ ∇ₓN_b) I) dV₀ (material plus geometric, or +//! initial-stress, stiffness), +//! +//! with Bₓ the small-strain operator of the current gradients. This is the push-forward of the +//! textbook ∫ B₀ᵀ S dV₀ and ∫ (B₀ᵀ C B₀ + ∇₀N·S∇₀N I) dV₀, identical term by term (Bonet & Wood +//! 2008, §9.2–9.4; K.-J. Bathe, *Finite Element Procedures*, 2nd ed. 2014, §6.3). A rigid motion +//! leaves τ at zero exactly, whatever its size. +//! +//! **Hybrid (mixed u–p) solids** with a decoupled material W = W̄(C̄) + U(J): the three-field +//! formulation of J. C. Simo, R. L. Taylor, K. S. Pister, *Variational and projection methods +//! for the volume constraint in finite deformation elasto-plasticity*, CMAME 51:177–208, 1985, +//! with the element's discontinuous pressure functions N_p (the linear hybrid elements' fields, +//! [`crate::hybrid`]): the dilatation θ = N_p·θ̂ is the projection of J, M θ̂ = ∫ N_p J dV₀ with +//! M = ∫ N_p N_pᵀ dV₀; the pressure p = N_p·p̂ the projection of U′(θ), M p̂ = ∫ N_p U′(θ) dV₀. +//! Then τ = τ_iso(F) + J p I, and with G = ∫ N_p J (∇ₓN)ᵀ dV₀ and H = ∫ N_p U″(θ) N_pᵀ dV₀, +//! K = ∫ (Bₓᵀ c* Bₓ + ∇ₓN·τ∇ₓN I) dV₀ + Gᵀ M⁻¹ H M⁻¹ G with c* = c_iso + J p (I ⊗ I − 2𝕀), +//! symmetric; at u = 0 it is the linear hybrid stiffness K_dev + κ G M⁻¹ Gᵀ. The constant +//! pressure case is the mean dilatation method of J. C. Nagtegaal, D. M. Parks, J. R. Rice +//! (CMAME 4:153–177, 1974) at finite strain. St Venant–Kirchhoff has no such split, so hybrid +//! elements of an elastic material are refused in nonlinear steps. +//! +//! **2D elements**: plane strain (F₃₃ = 1) and axisymmetric (F₃₃ = r/R, the hoop stretch, with +//! the hoop column N/r of the current radius in Bₓ and a geometric term τ_θθ N_a N_b / r² on the +//! radial motion; Bonet & Wood 2008, §9.5). Plane stress (whose thickness change needs a local +//! iteration for σ_zz = 0) is refused. +//! +//! **Incompatible modes** (hex8, quad4) are refused in nonlinear steps: their finite-strain form +//! (enhanced assumed strain, J. C. Simo, F. Armero, IJNME 33:1413–1449, 1992) needs the internal +//! modes as element unknowns carried between iterations, and it develops spurious hourglass +//! modes in large compression (P. Wriggers, S. Reese, CMAME 135:201–209, 1996), which want a +//! stabilised variant of their own. + +use crate::constitutive::{MaterialModel, Response, add_pressure, logarithmic_strain}; +use crate::dense::{Cholesky, Matrix, inverse3}; +use crate::element::{Formulation, PointStress, Solid, accumulate_symmetric, b_block}; +use crate::hybrid::pressure_basis; +use crate::plane::{Plane, Theory}; +use crate::rotation::{IDENTITY, Mat3}; +use crate::{Error, Result}; + +/// Internal forces (element degrees of freedom, node by node) and, when asked for, the tangent. +#[derive(Debug, Clone, PartialEq)] +pub struct Finite { + pub force: Vec, + pub tangent: Option, +} + +/// What one integration point gives the element. +struct PointKinematics { + weight: f64, + f: Mat3, + j: f64, + /// ∇ₓN of every node. + gradients: Vec<[f64; 3]>, + /// Reference position. + x0: [f64; 3], +} + +/// The checked history slice of one integration point. +fn point_history(history: &[f64], len: usize, point: usize) -> &[f64] { + history.get(point * len..(point + 1) * len).unwrap_or_default() +} + +/// Add f_a += w Bₓ,aᵀ τ and, with `k`, the material and geometric tangent. +fn add_point(force: &mut [f64], k: Option<&mut Matrix>, gradients: &[[f64; 3]], r: &Response, w: f64) { + for (a, g) in gradients.iter().enumerate() { + let b = b_block(g); + for c in 0..3 { + if let Some(x) = force.get_mut(3 * a + c) { + *x += w * (0..6).map(|s| b[s][c] * r.stress[s]).sum::(); + } + } + } + let Some(k) = k else { return }; + accumulate_symmetric(&mut k.data, gradients, &r.tangent, w); + let tau = crate::rotation::tensor(&r.stress); + let n = k.n; + for (a, ga) in gradients.iter().enumerate() { + let tg: [f64; 3] = std::array::from_fn(|i| (0..3).map(|jj| tau[i][jj] * ga[jj]).sum()); + for (b, gb) in gradients.iter().enumerate() { + let s = w * (tg[0] * gb[0] + tg[1] * gb[1] + tg[2] * gb[2]); + for c in 0..3 { + if let Some(x) = k.data.get_mut((3 * a + c) * n + 3 * b + c) { + *x += s; + } + } + } + } +} + +fn inverted(kind: &str, j: f64) -> Error { + Error::Degenerate(format!("a {kind} element is inverted at an integration point (J = det F = {j:.3e})")) +} + +/// One displacement vector per node, all finite. +fn check_displacements(u: &[[f64; N]], nodes: usize) -> Result<()> { + if u.len() != nodes { + return Err(Error::Invalid(format!("{} displacements for a {nodes}-node element", u.len()))); + } + if u.iter().flatten().any(|v| !v.is_finite()) { + return Err(Error::Invalid("displacements must be finite".into())); + } + Ok(()) +} + +/// Check the element's history slices against the material and its points. +fn check_history(m: &MaterialModel, points: usize, committed: &[f64], updated: &[f64]) -> Result { + let len = m.history_len(); + if committed.len() != points * len || updated.len() != points * len { + return Err(Error::Invalid(format!("{} history values for {points} points of {len}", committed.len()))); + } + Ok(len) +} + +impl Solid<'_> { + /// Refuse what has no finite-strain form here. + fn check_finite(&self, m: &MaterialModel, u: &[[f64; 3]]) -> Result<()> { + check_displacements(u, self.ty.nodes)?; + match self.formulation { + Formulation::IncompatibleModes => Err(Error::NotSupported(format!( + "{} incompatible modes in nonlinear steps (their finite-strain form, enhanced assumed strain, is not implemented; use full integration, hybrid or hex20)", + self.ty.kind.name() + ))), + Formulation::Hybrid { .. } if !m.is_decoupled() => Err(Error::NotSupported(format!( + "hybrid elements of a {} material in nonlinear steps (the mixed formulation needs a volumetric–isochoric split of the strain energy; give the material the neo-Hookean option)", + m.name() + ))), + _ => Ok(()), + } + } + + /// F, J and the current gradients at every stiffness point. + fn kinematics(&self, u: &[[f64; 3]]) -> Result> { + let n = self.ty.nodes; + let mut out = Vec::new(); + for p in self.rule().points() { + let g = self.gradients(p.xi)?; + let reference = g.dx.get(..n).unwrap_or_default(); + let f = crate::constitutive::deformation_gradient(u, reference); + let (inv, j) = inverse3(&f).ok_or_else(|| inverted(self.ty.kind.name(), 0.0))?; + if !(j > 0.0 && j.is_finite()) { + return Err(inverted(self.ty.kind.name(), j)); + } + let gradients = reference.iter().map(|d| std::array::from_fn(|i| (0..3).map(|jj| d[jj] * inv[jj][i]).sum())).collect(); + out.push(PointKinematics { weight: g.det * p.weight, f, j, gradients, x0: self.position_of(g.shape.values()) }); + } + Ok(out) + } + + /// The effective response (Kirchhoff stress and tangent) at every point, the forces, and the + /// tangent when `tangent`. + fn finite_points( + &self, + m: &MaterialModel, + u: &[[f64; 3]], + committed: &[f64], + updated: &mut [f64], + tangent: bool, + ) -> Result<(Finite, Vec<(PointKinematics, Response)>)> { + self.check_finite(m, u)?; + let points = self.kinematics(u)?; + let len = check_history(m, points.len(), committed, updated)?; + let ndof = 3 * self.ty.nodes; + let mut force = vec![0.0; ndof]; + let mut k = tangent.then(|| Matrix::zeros(ndof)); + let responses: Vec = if let Formulation::Hybrid { .. } = self.formulation { + self.mixed_responses(m, &points, k.as_mut())? + } else { + let mut out = Vec::with_capacity(points.len()); + for (i, p) in points.iter().enumerate() { + let updated_point = updated.get_mut(i * len..(i + 1) * len).unwrap_or_default(); + out.push(m.response(&p.f, point_history(committed, len, i), updated_point)?); + } + out + }; + for (p, r) in points.iter().zip(&responses) { + add_point(&mut force, k.as_mut(), &p.gradients, r, p.weight); + } + Ok((Finite { force, tangent: k }, points.into_iter().zip(responses).collect())) + } + + /// The three-field mixed responses τ_iso + J p I with their tangents, and the dilatation + /// stiffness Gᵀ M⁻¹ H M⁻¹ G added to `k`. + fn mixed_responses(&self, m: &MaterialModel, points: &[PointKinematics], k: Option<&mut Matrix>) -> Result> { + let (np, frame) = (self.np(), self.frame()); + let basis: Vec<[f64; 4]> = points.iter().map(|p| pressure_basis(np, p.x0, &frame)).collect(); + let mut mm = Matrix::zeros(np); + let mut theta = vec![0.0; np]; + for (p, nb) in points.iter().zip(&basis) { + for a in 0..np { + theta[a] += p.weight * nb[a] * p.j; + for b in 0..np { + mm.add(a, b, p.weight * nb[a] * nb[b]); + } + } + } + let chol = Cholesky::new(&mm)?; + chol.solve(&mut theta); + let mut pressure = vec![0.0; np]; + let mut h = Matrix::zeros(np); + for (p, nb) in points.iter().zip(&basis) { + let th: f64 = (0..np).map(|a| nb[a] * theta[a]).sum(); + if !(th > 0.0 && th.is_finite()) { + return Err(inverted(self.ty.kind.name(), th)); + } + let (dp, ddp) = m.volumetric(th).ok_or_else(|| Error::Invalid("volumetric response".into()))?; + for a in 0..np { + pressure[a] += p.weight * nb[a] * dp; + for b in 0..np { + h.add(a, b, p.weight * ddp * nb[a] * nb[b]); + } + } + } + chol.solve(&mut pressure); + let mut out = Vec::with_capacity(points.len()); + for (p, nb) in points.iter().zip(&basis) { + let pr: f64 = (0..np).map(|a| nb[a] * pressure[a]).sum(); + out.push(add_pressure(&m.isochoric(&p.f)?, p.j, pr, 0.0)); + } + if let Some(k) = k { + // A = M⁻¹ G, np × ndof; K += Aᵀ H A. + let ndof = k.n; + let mut a_cols = vec![vec![0.0; np]; ndof]; + for (p, nb) in points.iter().zip(&basis) { + for (node, g) in p.gradients.iter().enumerate() { + for c in 0..3 { + if let Some(col) = a_cols.get_mut(3 * node + c) { + for (x, n) in col.iter_mut().zip(nb) { + *x += p.weight * n * p.j * g[c]; + } + } + } + } + } + for col in &mut a_cols { + chol.solve(col); + } + let ha: Vec> = a_cols.iter().map(|col| h.mul_vec(col)).collect(); + for i in 0..ndof { + for jj in i..ndof { + let v: f64 = a_cols.get(i).zip(ha.get(jj)).map_or(0.0, |(x, y)| x.iter().zip(y).map(|(p, q)| p * q).sum()); + k.add(i, jj, v); + if i != jj { + k.add(jj, i, v); + } + } + } + } + Ok(out) + } + + /// Internal forces and (with `tangent`) the tangent stiffness at nodal displacements `u`, + /// with the material history `committed` at the start of the increment (points × the + /// material's history length) written, updated, to `updated`. + pub fn finite_strain(&self, m: &MaterialModel, u: &[[f64; 3]], committed: &[f64], updated: &mut [f64], tangent: bool) -> Result { + Ok(self.finite_points(m, u, committed, updated, tangent)?.0) + } + + /// The Cauchy stress σ = τ/J and the logarithmic strain at every integration point (at its + /// reference position). + pub fn finite_strain_stresses(&self, m: &MaterialModel, u: &[[f64; 3]], committed: &[f64]) -> Result> { + let mut scratch = committed.to_vec(); + let (_, points) = self.finite_points(m, u, committed, &mut scratch, false)?; + points + .iter() + .map(|(p, r)| { + let strain = logarithmic_strain(&p.f).ok_or_else(|| inverted(self.ty.kind.name(), p.j))?; + Ok(PointStress { xyz: p.x0, strain, stress: r.stress.map(|v| v / p.j) }) + }) + .collect() + } +} + +impl Plane<'_> { + fn check_finite(&self, u: &[[f64; 2]]) -> Result<()> { + check_displacements(u, self.ty.nodes)?; + if let Theory::PlaneStress { .. } = self.theory { + return Err(Error::NotSupported( + "plane stress elements in nonlinear steps (the thickness change needs a local iteration; use plane strain, axisymmetric or solid elements)" + .into(), + )); + } + if self.formulation == Formulation::IncompatibleModes { + return Err(Error::NotSupported(format!( + "2D {} incompatible modes in nonlinear steps (their finite-strain form, enhanced assumed strain, is not implemented; use full integration or quad8)", + self.ty.kind.name() + ))); + } + Ok(()) + } + + /// Per point: (dv, F, J, current columns (∂N/∂x, ∂N/∂y, N/r), reference position, response). + #[allow(clippy::type_complexity)] + fn finite_points( + &self, + m: &MaterialModel, + u: &[[f64; 2]], + committed: &[f64], + updated: &mut [f64], + tangent: bool, + ) -> Result<(Finite, Vec<(Mat3, f64, [f64; 3], Response)>)> { + self.check_finite(u)?; + let n = self.ty.nodes; + let rule = self.rule().points(); + let len = check_history(m, rule.len(), committed, updated)?; + let ndof = 2 * n; + let mut force = vec![0.0; ndof]; + let mut k = tangent.then(|| Matrix::zeros(ndof)); + let axisymmetric = self.theory == Theory::Axisymmetric; + let mut out = Vec::with_capacity(rule.len()); + for (i, p) in rule.iter().enumerate() { + let at = self.at(p)?; + let cols = at.g.get(..n).unwrap_or_default(); + let mut f = IDENTITY; + for (ua, g) in u.iter().zip(cols) { + for r in 0..2 { + for c in 0..2 { + f[r][c] += ua[r] * g[c]; + } + } + if axisymmetric { + f[2][2] += ua[0] * g[2]; + } + } + let det2 = f[0][0] * f[1][1] - f[0][1] * f[1][0]; + let j = det2 * f[2][2]; + if !(j > 0.0 && j.is_finite() && f[2][2] > 0.0) { + return Err(inverted(&format!("2D {}", self.ty.kind.name()), j)); + } + let inv = [[f[1][1] / det2, -f[0][1] / det2], [-f[1][0] / det2, f[0][0] / det2]]; + let current: Vec<[f64; 3]> = + cols.iter().map(|g| [g[0] * inv[0][0] + g[1] * inv[1][0], g[0] * inv[0][1] + g[1] * inv[1][1], g[2] / f[2][2]]).collect(); + let updated_point = updated.get_mut(i * len..(i + 1) * len).unwrap_or_default(); + let r = m.response(&f, point_history(committed, len, i), updated_point)?; + let tau = [r.stress[0], r.stress[1], r.stress[2], r.stress[3]]; + for (a, g) in current.iter().enumerate() { + if let Some(x) = force.get_mut(2 * a) { + *x += at.dv * (g[0] * tau[0] + g[2] * tau[2] + g[1] * tau[3]); + } + if let Some(x) = force.get_mut(2 * a + 1) { + *x += at.dv * (g[1] * tau[1] + g[0] * tau[3]); + } + } + if let Some(k) = k.as_mut() { + let c4: [[f64; 4]; 4] = std::array::from_fn(|a| std::array::from_fn(|b| r.tangent[a][b])); + crate::plane::accumulate(&mut k.data, ¤t, &c4, at.dv); + for (a, ga) in current.iter().enumerate() { + for (b, gb) in current.iter().enumerate() { + let s = ga[0] * gb[0] * tau[0] + ga[1] * gb[1] * tau[1] + (ga[0] * gb[1] + ga[1] * gb[0]) * tau[3]; + k.add(2 * a, 2 * b, at.dv * (s + tau[2] * ga[2] * gb[2])); + k.add(2 * a + 1, 2 * b + 1, at.dv * s); + } + } + } + out.push((f, j, at.x, r)); + } + Ok((Finite { force, tangent: k }, out)) + } + + /// Internal forces (u, v per node) and (with `tangent`) the tangent stiffness at nodal + /// displacements `u`, as for [`Solid::finite_strain`]. + pub fn finite_strain(&self, m: &MaterialModel, u: &[[f64; 2]], committed: &[f64], updated: &mut [f64], tangent: bool) -> Result { + Ok(self.finite_points(m, u, committed, updated, tangent)?.0) + } + + /// The Cauchy stress and logarithmic strain at every integration point, as the solid + /// elements' 6-vectors (axisymmetric: rr, axial, hoop, r-axial). + pub fn finite_strain_stresses(&self, m: &MaterialModel, u: &[[f64; 2]], committed: &[f64]) -> Result> { + let mut scratch = committed.to_vec(); + let (_, points) = self.finite_points(m, u, committed, &mut scratch, false)?; + points + .iter() + .map(|(f, j, x, r)| { + let strain = logarithmic_strain(f).ok_or_else(|| inverted(&format!("2D {}", self.ty.kind.name()), *j))?; + Ok(PointStress { xyz: *x, strain, stress: r.stress.map(|v| v / j) }) + }) + .collect() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::material::isotropic; + use fieldcraft_model::ElementKind; + + fn hex20() -> Vec<[f64; 3]> { + let mut x: Vec<[f64; 3]> = crate::shape::HEX_CORNERS.iter().map(|c| [c[0] + 1.0, 0.5 * c[1] + 1.0, 0.8 * c[2]]).collect(); + let edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), (0, 4), (1, 5), (2, 6), (3, 7)]; + for (a, b) in edges { + x.push(std::array::from_fn(|k| 0.5 * (x[a][k] + x[b][k]))); + } + // Distort a corner and a mid-side node. + x[6][0] += 0.2; + x[9][2] += 0.05; + x + } + + fn tet10() -> Vec<[f64; 3]> { + let c = [[0.0, 0.0, 0.0], [1.2, 0.1, 0.0], [0.2, 1.0, 0.1], [0.1, 0.2, 0.9]]; + let mut x = c.to_vec(); + for (a, b) in [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)] { + x.push(std::array::from_fn(|k| 0.5 * (c[a][k] + c[b][k]))); + } + x + } + + fn displacement(n: usize, scale: f64) -> Vec<[f64; 3]> { + (0..n).map(|i| [scale * (i as f64 * 1.3).sin(), scale * (i as f64 * 0.7).cos(), scale * (i as f64 * 2.1).sin()]).collect() + } + + /// Central differences of the internal forces against the tangent. + fn check_tangent(solid: &Solid, m: &MaterialModel, u: &[[f64; 3]]) { + let k = solid.finite_strain(m, u, &[], &mut [], true).unwrap().tangent.unwrap(); + let scale = k.data.iter().fold(0.0f64, |a, v| a.max(v.abs())); + let n = 3 * u.len(); + let h = 1e-6; + for j in 0..n { + let shifted = |s: f64| { + let mut v = u.to_vec(); + v[j / 3][j % 3] += s; + solid.finite_strain(m, &v, &[], &mut [], false).unwrap().force + }; + let (p, q) = (shifted(h), shifted(-h)); + for i in 0..n { + let numeric = (p[i] - q[i]) / (2.0 * h); + assert!((numeric - k.get(i, j)).abs() < 1e-6 * scale, "{:?} {i} {j}: {numeric} vs {}", solid.ty.kind, k.get(i, j)); + } + } + for i in 0..n { + for j in 0..n { + assert!((k.get(i, j) - k.get(j, i)).abs() <= 1e-12 * scale, "symmetric"); + } + } + } + + #[test] + fn solid_tangents_are_the_derivative_of_the_forces() { + let d = isotropic(1000.0, 0.3).unwrap(); + let neo = MaterialModel::NeoHookean { shear: 300.0, bulk: 3000.0 }; + let svk = MaterialModel::StVenantKirchhoff(Box::new(d)); + let (x20, x10) = (hex20(), tet10()); + let x8: Vec<[f64; 3]> = x20[..8].to_vec(); + let x6 = vec![[0., 0., 0.], [1., 0., 0.], [0., 1., 0.1], [0., 0., 1.], [1.1, 0., 1.], [0., 1., 1.2]]; + let x4 = vec![x10[0], x10[1], x10[2], x10[3]]; + for (kind, x) in + [(ElementKind::Tet4, &x4), (ElementKind::Tet10, &x10), (ElementKind::Wedge6, &x6), (ElementKind::Hex8, &x8), (ElementKind::Hex20, &x20)] + { + for m in [&svk, &neo] { + let solid = Solid::new(kind, x, Formulation::Full).unwrap(); + check_tangent(&solid, m, &displacement(x.len(), 0.08)); + } + } + let reduced = Solid::new(ElementKind::Hex20, &x20, Formulation::Reduced).unwrap(); + check_tangent(&reduced, &neo, &displacement(20, 0.08)); + for (kind, x, np) in [(ElementKind::Hex8, &x8, 1), (ElementKind::Hex20, &x20, 4), (ElementKind::Tet10, &x10, 1)] { + let hybrid = Solid::new(kind, x, Formulation::Hybrid { pressures: np }).unwrap(); + check_tangent(&hybrid, &neo, &displacement(x.len(), 0.08)); + assert!(hybrid.finite_strain(&svk, &displacement(x.len(), 0.0), &[], &mut [], false).is_err(), "SVK has no split"); + } + let im = Solid::new(ElementKind::Hex8, &x8, Formulation::IncompatibleModes).unwrap(); + assert!(matches!(im.finite_strain(&svk, &displacement(8, 0.0), &[], &mut [], true), Err(Error::NotSupported(_)))); + } + + #[test] + fn at_zero_displacement_the_tangent_is_the_linear_stiffness() { + let (e, nu) = (1000.0, 0.45); + let d = isotropic(e, nu).unwrap(); + let neo = MaterialModel::NeoHookean { shear: e / (2.0 * (1.0 + nu)), bulk: e / (3.0 * (1.0 - 2.0 * nu)) }; + let x20 = hex20(); + for (f, m) in [ + (Formulation::Full, MaterialModel::StVenantKirchhoff(Box::new(d))), + (Formulation::Full, neo.clone()), + (Formulation::Hybrid { pressures: 4 }, neo), + ] { + let solid = Solid::new(ElementKind::Hex20, &x20, f).unwrap(); + let r = solid.finite_strain(&m, &vec![[0.0; 3]; 20], &[], &mut [], true).unwrap(); + let linear = solid.stiffness(&d).unwrap(); + let k = r.tangent.unwrap(); + for (a, b) in k.data.iter().zip(&linear.data) { + assert!((a - b).abs() < 1e-9 * e, "{f:?}: {a} vs {b}"); + } + assert!(r.force.iter().all(|v| v.abs() < 1e-12 * e)); + } + } + + #[test] + fn rigid_rotation_strains_and_stresses_nothing() { + let x20 = hex20(); + let rot = [[0.0, -0.6, 0.8], [0.6, 0.64, 0.48], [-0.8, 0.48, 0.36]]; + let u: Vec<[f64; 3]> = x20.iter().map(|p| std::array::from_fn(|i| (0..3).map(|jj| rot[i][jj] * p[jj]).sum::() + 5.0 - p[i])).collect(); + let neo = MaterialModel::NeoHookean { shear: 300.0, bulk: 3000.0 }; + for f in [Formulation::Full, Formulation::Hybrid { pressures: 4 }] { + let solid = Solid::new(ElementKind::Hex20, &x20, f).unwrap(); + let s = solid.finite_strain_stresses(&neo, &u, &[]).unwrap(); + assert!(s.iter().all(|p| p.stress.iter().chain(&p.strain).all(|v| v.abs() < 1e-12 * 3000.0)), "{f:?}: {:?}", s[0]); + let r = solid.finite_strain(&neo, &u, &[], &mut [], false).unwrap(); + assert!(r.force.iter().all(|v| v.abs() < 1e-11 * 3000.0)); + } + // A rotation that inverts the element is refused. + let flip: Vec<[f64; 3]> = x20.iter().map(|p| [-2.0 * p[0], 0.0, 0.0]).collect(); + let solid = Solid::new(ElementKind::Hex20, &x20, Formulation::Full).unwrap(); + assert!(matches!(solid.finite_strain(&neo, &flip, &[], &mut [], false), Err(Error::Degenerate(_)))); + } + + #[test] + fn plane_strain_and_axisymmetric_tangents_are_the_derivative_of_the_forces() { + let quad8 = [ + [1.0, 0.0, 0.0], + [2.0, 0.1, 0.0], + [2.2, 1.0, 0.0], + [0.9, 1.1, 0.0], + [1.5, 0.0, 0.0], + [2.1, 0.55, 0.0], + [1.5, 1.1, 0.0], + [0.95, 0.5, 0.0], + ]; + let neo = MaterialModel::NeoHookean { shear: 300.0, bulk: 3000.0 }; + let svk = MaterialModel::StVenantKirchhoff(Box::new(isotropic(1000.0, 0.3).unwrap())); + for theory in [Theory::PlaneStrain { thickness: 2.0 }, Theory::Axisymmetric] { + for (kind, x) in [(ElementKind::Quad8, &quad8[..]), (ElementKind::Quad4, &quad8[..4]), (ElementKind::Tri3, &quad8[..3])] { + let plane = Plane::new(kind, x, Formulation::Full, theory).unwrap(); + for m in [&neo, &svk] { + let u: Vec<[f64; 2]> = displacement(x.len(), 0.06).iter().map(|v| [v[0], v[1]]).collect(); + let k = plane.finite_strain(m, &u, &[], &mut [], true).unwrap().tangent.unwrap(); + let scale = k.data.iter().fold(0.0f64, |a, v| a.max(v.abs())); + let h = 1e-6; + for j in 0..2 * x.len() { + let shifted = |s: f64| { + let mut v = u.clone(); + v[j / 2][j % 2] += s; + plane.finite_strain(m, &v, &[], &mut [], false).unwrap().force + }; + let (p, q) = (shifted(h), shifted(-h)); + for i in 0..2 * x.len() { + let numeric = (p[i] - q[i]) / (2.0 * h); + assert!((numeric - k.get(i, j)).abs() < 1e-6 * scale, "{theory:?} {kind:?} {i} {j}: {numeric} vs {}", k.get(i, j)); + } + } + } + // At rest the tangent is the linear stiffness. + let d = isotropic(1000.0, 0.3).unwrap(); + let k = plane.finite_strain(&svk, &vec![[0.0; 2]; x.len()], &[], &mut [], true).unwrap().tangent.unwrap(); + for (a, b) in k.data.iter().zip(&plane.stiffness(&d).unwrap().data) { + assert!((a - b).abs() < 1e-9 * 1000.0); + } + } + } + let stress = Plane::new(ElementKind::Quad4, &quad8[..4], Formulation::Full, Theory::PlaneStress { thickness: 1.0 }).unwrap(); + assert!(matches!(stress.finite_strain(&svk, &[[0.0; 2]; 4], &[], &mut [], true), Err(Error::NotSupported(_)))); + } +} diff --git a/crates/fem/src/hybrid.rs b/crates/fem/src/hybrid.rs index 1814d2c..b140b71 100644 --- a/crates/fem/src/hybrid.rs +++ b/crates/fem/src/hybrid.rs @@ -111,7 +111,7 @@ impl PressureField { /// 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] { +pub(crate) 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 { @@ -124,7 +124,7 @@ fn pressure_basis(np: usize, x: [f64; 3], frame: &([f64; 3], f64)) -> [f64; 4] { impl Solid<'_> { /// The pressure unknowns of this (hybrid) element. - fn np(&self) -> usize { + pub(crate) fn np(&self) -> usize { match self.formulation { Formulation::Hybrid { pressures } => pressures.min(4), _ => 0, @@ -133,7 +133,7 @@ impl Solid<'_> { /// 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) { + pub(crate) 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]); diff --git a/crates/fem/src/lib.rs b/crates/fem/src/lib.rs index aaf88d8..c2d15da 100644 --- a/crates/fem/src/lib.rs +++ b/crates/fem/src/lib.rs @@ -13,8 +13,13 @@ //! - [`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; +//! - [`constitutive`]: the material point of finite-strain analysis (Kirchhoff stress and +//! spatial tangent from the deformation gradient, with history): St Venant–Kirchhoff and +//! neo-Hookean; //! - [`element`]: the element table, stiffness (full, reduced, incompatible modes), consistent //! and lumped mass, body forces, integration-point stresses and their extrapolation to nodes; +//! - [`finite`]: total Lagrangian solid and 2D elements (internal forces, the consistent tangent +//! with its geometric part, mixed u–p at finite strain, Cauchy stress and logarithmic strain); //! - [`plane`]: 2D elements (plane stress, plane strain, axisymmetric) on tri3, tri6, quad4 and //! quad8, with two displacement components per node; //! - [`shell`]: flat MITC4 quadrilateral and MITC3 triangle shells with a drilling rotation, six @@ -44,10 +49,12 @@ pub mod basis; pub mod beam; +pub mod constitutive; pub mod dense; pub mod dof; pub mod element; pub mod face; +pub mod finite; pub mod frame; pub mod heat; pub mod hybrid; diff --git a/crates/fem/src/plane.rs b/crates/fem/src/plane.rs index 4a5cfe2..7c1a2c5 100644 --- a/crates/fem/src/plane.rs +++ b/crates/fem/src/plane.rs @@ -189,18 +189,18 @@ pub fn formulation(kind: ElementKind, f: &fieldcraft_model::Formulation) -> Resu } /// One column of the strain operator: (∂N/∂x, ∂N/∂y, hoop N/r), for a node or an internal mode. -type Column = [f64; 3]; +pub(crate) type Column = [f64; 3]; /// What the kernels need at one integration point. -struct At { - shape: Shape, +pub(crate) struct At { + pub(crate) shape: Shape, /// Strain columns of the nodes. - g: [Column; MAX_NODES], - det: f64, - x: [f64; 3], + pub(crate) g: [Column; MAX_NODES], + pub(crate) det: f64, + pub(crate) x: [f64; 3], /// det J times the point's weight times the out-of-plane extent: the point's share of the /// volume. - dv: f64, + pub(crate) dv: f64, } /// The incompatible modes' count of internal degrees of freedom: two modes, two directions. @@ -210,14 +210,14 @@ const INCOMPATIBLE_DOFS: usize = 4; #[derive(Debug, Clone, Copy)] pub struct Plane<'a> { pub ty: &'static ElementType, - xyz: &'a [[f64; 3]], - formulation: Formulation, - theory: Theory, + pub(crate) xyz: &'a [[f64; 3]], + pub(crate) formulation: Formulation, + pub(crate) theory: Theory, } /// out[2i + r][2j + c] += w (Bᵢᵀ D Bⱼ)[r][c] over every pair of columns (row-major, `cols` wide), /// computed for j ≥ i and mirrored, so the result is exactly symmetric. -fn accumulate(out: &mut [f64], columns: &[Column], d: &Elasticity2, w: f64) { +pub(crate) fn accumulate(out: &mut [f64], columns: &[Column], d: &Elasticity2, w: f64) { let cols = 2 * columns.len(); let block = |g: &Column| -> [[f64; 2]; 4] { [[g[0], 0.0], [0.0, g[1]], [g[2], 0.0], [g[1], g[0]]] }; for (j, gj) in columns.iter().enumerate() { @@ -293,7 +293,7 @@ impl<'a> Plane<'a> { } } - fn degenerate(&self) -> Error { + pub(crate) fn degenerate(&self) -> Error { Error::Degenerate(format!("a 2D {} element is inverted, clockwise or degenerate (its Jacobian is not positive)", self.ty.kind.name())) } @@ -310,7 +310,7 @@ impl<'a> Plane<'a> { (j, j[0][0] * j[1][1] - j[0][1] * j[1][0]) } - fn at(&self, p: &Point) -> Result { + pub(crate) fn at(&self, p: &Point) -> Result { let basis = Basis::of(self.ty.kind); let shape = basis.eval(p.xi); let grads = basis.gradients(&shape, self.xyz).map_err(|_| self.degenerate())?; diff --git a/crates/model/src/nodal.rs b/crates/model/src/nodal.rs index e3dfa1c..12aed36 100644 --- a/crates/model/src/nodal.rs +++ b/crates/model/src/nodal.rs @@ -124,6 +124,55 @@ pub fn integrate_face(xyz: &[[f64; 3]], mut visit: impl FnMut(&[f64], [f64; 3])) Ok(()) } +/// A quadrature point of a reference face or edge: the shape function values, their +/// derivatives along the face's natural coordinates ξ and η (an edge has ξ only; its `deta` is +/// empty) and the weight. +#[derive(Debug, Clone, PartialEq)] +pub struct QuadraturePoint { + pub n: Vec, + pub dxi: Vec, + pub deta: Vec, + pub weight: f64, +} + +/// The quadrature points of a face of `nodes` nodes (3, 4, 6 or 8), with the rule and shape +/// functions of [`integrate_face`]: what a load on the deformed face (follower pressure) +/// integrates with. +pub fn face_points(nodes: usize) -> Result> { + let (corners, shape): (usize, fn(usize, f64, f64) -> Shape) = match nodes { + 3 | 6 => (3, triangle), + 4 | 8 => (4, quadrilateral), + n => return Err(Error::Invalid(format!("a face has 3, 4, 6 or 8 nodes, not {n}"))), + }; + Ok(rule(corners) + .into_iter() + .map(|(xi, eta, weight)| { + let s = shape(nodes, xi, eta); + QuadraturePoint { n: s.n, dxi: s.dxi, deta: s.deta, weight } + }) + .collect()) +} + +/// The quadrature points of an element edge of `nodes` nodes (2 or 3), with the rule and shape +/// functions of [`edge_normal_weights`]. +pub fn edge_points(nodes: usize) -> Result> { + if !matches!(nodes, 2 | 3) { + return Err(Error::Invalid(format!("an edge has 2 or 3 nodes, not {nodes}"))); + } + let g = (0.6f64).sqrt(); + Ok([(-g, 5.0 / 9.0), (0.0, 8.0 / 9.0), (g, 5.0 / 9.0)] + .into_iter() + .map(|(s, weight)| { + let (n, dn) = if nodes == 2 { + (vec![(1.0 - s) / 2.0, (1.0 + s) / 2.0], vec![-0.5, 0.5]) + } else { + (vec![s * (s - 1.0) / 2.0, s * (s + 1.0) / 2.0, 1.0 - s * s], vec![s - 0.5, s + 0.5, -2.0 * s]) + }; + QuadraturePoint { n, dxi: dn, deta: Vec::new(), weight } + }) + .collect()) +} + /// ∫ Nᵢ n dA for every node of a face: the consistent nodal forces of a unit pressure pulling /// along the face's right-hand (outward, for solids) normal. pub fn face_normal_weights(xyz: &[[f64; 3]]) -> Result> {