1b-agent/OneCAD-Simulation
1
0 Fork 0
forked from 1b-agent/fieldcraft-upstream
Code Pull Requests Activity

Explicit X2: large deformation, J2 and Johnson-Cook plasticity, damage and erosion (M10)

Browse Source 
- model: materials gain solver-neutral plastic (bilinear isotropic/kinematic,
  hardening table, Johnson-Cook with adiabatic heating), damage (Johnson-Cook
  D1-D5) and erosion (plastic strain); explicit steps gain nonlinear_geometry
  (default on) and stress_rate (jaumann | green_naghdi).
- fem: rotation (Hughes-Winget, Flanagan-Taylor), plasticity (radial return,
  capped bracketed Newton for Johnson-Cook, damage; the update takes the
  mechanical strain increment so thermal strains can be subtracted), closed-form
  hex8 gradients for deformed configurations.
- explicit: updated Lagrangian hex8/tet4, per-element material state, erosion
  with reported energy, plastic work in the energy account.
- solver: warns that the linear solver uses the elastic constants of plastic
  materials (results hash-identical).
- decks: Abaqus *PLASTIC (+ Johnson-Cook keywords), Nastran MATS1/TABLES1,
  LS-DYNA *MAT_PLASTIC_KINEMATIC / *MAT_PIECEWISE_LINEAR_PLASTICITY; what a
  format cannot express is reported as Unsupported.
- engine: model.material.plasticity, model.material.damage (catalogue ready),
  explicit element probe, eroded elements in the solve summary.

Step commits before the rebase onto the DOF framework: branch explicit-x2-steps.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm
This commit is contained in:
jyvikoandClaude Opus 5.5 authored and Irakli Kourtis committed 2026-10-13 20:52:57 -04:00
1 parent 4c0fcc80ad
commit 50764def18
38 files changed
+2839 -254
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

crates/fem/src/lib.rs
+6
View File
@@ -20,6 +20,10 @@
//! - [`face`]: consistent nodal forces of pressure on element faces;
//! - [`uniform`]: uniform-strain (one-point) hexahedron and tetrahedron kinematics, hourglass
//! shape vectors and element frequency bounds for explicit dynamics;
//! - [`rotation`]: incrementally objective stress updates (Hughes–Winget Jaumann, Flanagan–Taylor
//! Green–Naghdi);
//! - [`plasticity`]: J2 radial return with isotropic, kinematic and Johnson–Cook hardening,
//! Johnson–Cook damage;
//! - [`dense`]: the small dense linear algebra inside an element.
//!
//! Conventions: element matrices are ordered node by node, each node's degrees of freedom in
@@ -38,7 +42,9 @@ pub mod face;
pub mod jet;
pub mod material;
pub mod plane;
pub mod plasticity;
pub mod quadrature;
pub mod rotation;
pub mod shape;
pub mod uniform;
crates/fem/src/plasticity.rs
+484
View File
@@ -0,0 +1,484 @@
//! Hypoelastic–plastic stress update of an isotropic material: von Mises (J2) plasticity with
//! isotropic, linear kinematic or Johnson–Cook hardening, integrated by radial return, and
//! Johnson–Cook damage.
//!
//! **Radial return** (J. C. Simo, T. J. R. Hughes, *Computational Inelasticity*, Springer 1998,
//! §3.3, Box 3.1 and 3.2; R. D. Krieg, D. B. Krieg, J. Press. Vessel Technol. 99(4):510–515,
//! 1977, for its accuracy): with shear modulus G and bulk modulus K, the pressure takes the
//! volumetric increment elastically, p ← p + K Δθ, and the deviator the elastic trial
//! s* = s + 2G Δe. With the back stress α (kinematic hardening), ξ* = s* − α and the trial
//! von Mises stress q* = √(3/2)‖ξ*‖. If q* exceeds the current yield stress, the equivalent
//! plastic strain increment Δε_p solves the consistency condition
//! q* − 3G Δε_p − H_k Δε_p = σ_y(ε_p + Δε_p, Δε_p / Δt, T), and with the flow direction
//! n = ξ*/‖ξ*‖: s = s* − 2G √(3/2) Δε_p n, α ← α + √(2/3) H_k Δε_p n (Prager's linear kinematic
//! rule; H_k the kinematic modulus, so the uniaxial stress rises by H_k Δε_p either way).
//! - linear hardening: Δε_p = (q* − σ_y)/(3G + H) in closed form;
//! - a hardening table: piecewise linear, solved exactly segment by segment;
//! - Johnson–Cook (rate and temperature dependent): a bracketed scalar Newton iteration with
//! bisection safeguard on [0, (q* − σ_y)/(3G)] (the root is in it because σ_y only grows with
//! Δε_p), capped at [`MAX_ITERATIONS`].
//!
//! **Johnson–Cook** flow stress σ_y = (A + B ε_pⁿ)(1 + C ln ε̇*)(1 − T*ᵐ) (G. R. Johnson,
//! W. H. Cook, Proc. 7th Int. Symp. Ballistics, 1983) with ε̇* = ε̇_p / ε̇₀ taken as at least 1
//! (the logarithm stays ≥ 0, which also removes its singularity at rest) and
//! T* = (T − T_room)/(T_melt − T_room) clamped to [0, 1]. The temperature is that at the start
//! of the increment; plastic work heats the material adiabatically, ΔT = χ q Δε_p / (ρ c_p).
//!
//! **Johnson–Cook damage** (G. R. Johnson, W. H. Cook, Eng. Fract. Mech. 21(1):31–48, 1985):
//! D ← D + Δε_p / ε_f with ε_f = (D1 + D2 e^{D3 σ*})(1 + D4 ln ε̇*)(1 + D5 T*), σ* the mean
//! stress over the von Mises stress after the increment. A fracture strain below
//! [`MIN_FRACTURE_STRAIN`] counts as that (some published sets give negative values at high
//! triaxiality, where the point fails at once). The point has failed at D ≥ 1.
//!
//! The strain increment the update takes is the **mechanical** one: a caller with initial or
//! thermal strains subtracts their increment first.
//!
//! Stress 6-vectors hold tensor components (xx, yy, zz, xy, xz, yz), strain increments
//! engineering shears.
use modelcraft_model::{Damage, Elastic, Hardening, Plastic};
use crate::material::von_mises;
use crate::{Error, Result};
/// Most iterations of the rate-dependent return mapping.
pub const MAX_ITERATIONS: u32 = 100;
/// The smallest Johnson–Cook fracture strain.
pub const MIN_FRACTURE_STRAIN: f64 = 1e-6;
/// The flow stress of a plastic material.
#[derive(Debug, Clone, PartialEq)]
pub enum Flow {
/// σ_y = y0 + h ε_p (isotropic) or a surface of radius y0 moving with modulus h (kinematic).
Linear {
y0: f64,
h: f64,
kinematic: bool,
},
/// Isotropic hardening by (plastic strain, yield stress) points, constant after the last.
Table(Vec<(f64, f64)>),
JohnsonCook {
a: f64,
b: f64,
n: f64,
c: f64,
m: f64,
rate0: f64,
room: f64,
melt: f64,
},
}
/// A hypoelastic, possibly plastic, isotropic material.
#[derive(Debug, Clone, PartialEq)]
pub struct J2 {
pub shear: f64,
pub bulk: f64,
pub flow: Option<Flow>,
/// Johnson–Cook damage constants D1–D5.
pub damage: Option<[f64; 5]>,
/// χ / c_p: the temperature rise per plastic work per mass (0 without heating).
pub heating: f64,
}
/// The state of a material point.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PointState {
pub stress: [f64; 6],
pub back: [f64; 6],
pub plastic_strain: f64,
pub temperature: f64,
pub damage: f64,
}
impl PointState {
pub fn new(temperature: f64) -> PointState {
PointState { stress: [0.0; 6], back: [0.0; 6], plastic_strain: 0.0, temperature, damage: 0.0 }
}
}
/// What one update did.
#[derive(Debug, Clone, Copy, Default, PartialEq)]
pub struct Increment {
/// The equivalent plastic strain increment.
pub plastic_strain: f64,
/// Plastic work per volume, q Δε_p.
pub plastic_work: f64,
}
/// ‖s‖² of a symmetric tensor's 6-vector (tensor components).
fn norm2(s: &[f64; 6]) -> f64 {
s[0] * s[0] + s[1] * s[1] + s[2] * s[2] + 2.0 * (s[3] * s[3] + s[4] * s[4] + s[5] * s[5])
}
impl Flow {
/// The flow stress at plastic strain `e`, plastic strain rate `rate` and temperature `t`, and
/// its derivative with respect to the plastic strain increment over an increment `dt` (for
/// the Newton iteration).
fn stress(&self, e: f64, rate: f64, t: f64, dt: f64) -> (f64, f64) {
match *self {
Flow::Linear { y0, h, kinematic } => {
if kinematic {
(y0, 0.0)
} else {
(y0 + h * e, h)
}
}
Flow::Table(ref points) => table(points, e),
Flow::JohnsonCook { a, b, n, c, m, rate0, room, melt } => {
let hardening = a + b * e.max(0.0).powf(n);
let d_hardening = if e > 0.0 {
b * n * e.powf(n - 1.0)
} else if n < 1.0 && b > 0.0 {
f64::INFINITY
} else {
b * n
};
let ratio = rate / rate0;
let (rate_term, d_rate) = if ratio > 1.0 && dt > 0.0 { (1.0 + c * ratio.ln(), c / (rate * dt)) } else { (1.0, 0.0) };
let homologous = ((t - room) / (melt - room)).clamp(0.0, 1.0);
let thermal = if homologous > 0.0 { 1.0 - homologous.powf(m) } else { 1.0 };
(hardening * rate_term * thermal, (d_hardening * rate_term + hardening * d_rate) * thermal)
}
}
}
/// The kinematic hardening modulus.
fn kinematic(&self) -> f64 {
match *self {
Flow::Linear { h, kinematic: true, .. } => h,
_ => 0.0,
}
}
/// The rate sensitivity's reference rate and the homologous temperature (for damage).
fn rate_and_temperature(&self, t: f64) -> (f64, f64) {
match *self {
Flow::JohnsonCook { rate0, room, melt, .. } => (rate0, ((t - room) / (melt - room)).clamp(0.0, 1.0)),
_ => (1.0, 0.0),
}
}
}
/// The yield stress of a hardening table at plastic strain `e`, and its slope.
fn table(points: &[(f64, f64)], e: f64) -> (f64, f64) {
for w in points.windows(2) {
if let [(e0, s0), (e1, s1)] = *w
&& e < e1
{
let slope = (s1 - s0) / (e1 - e0);
return (s0 + slope * (e - e0).max(0.0), slope);
}
}
(points.last().map_or(0.0, |p| p.1), 0.0)
}
impl J2 {
/// The material of a model's isotropic elastic constants with its plasticity and damage.
pub fn from_model(elastic: &Elastic, plastic: Option<&Plastic>, damage: Option<&Damage>) -> Result<J2> {
let Elastic::Isotropic { young, poisson } = *elastic else {
return Err(Error::NotSupported("plasticity of an orthotropic material".into()));
};
let shear = young / (2.0 * (1.0 + poisson));
let bulk = young / (3.0 * (1.0 - 2.0 * poisson));
let mut heating = 0.0;
let flow = plastic.map(|p| match p {
Plastic::Bilinear { yield_stress, hardening_modulus, hardening } => {
Flow::Linear { y0: *yield_stress, h: *hardening_modulus, kinematic: *hardening == Hardening::Kinematic }
}
Plastic::Table { points } => Flow::Table(points.iter().map(|p| (p.strain, p.stress)).collect()),
Plastic::JohnsonCook { a, b, n, c, m, reference_strain_rate, room_temperature, melt_temperature, specific_heat, taylor_quinney } => {
heating = specific_heat.map_or(0.0, |cp| taylor_quinney / cp);
Flow::JohnsonCook {
a: *a,
b: *b,
n: *n,
c: *c,
m: *m,
rate0: *reference_strain_rate,
room: *room_temperature,
melt: *melt_temperature,
}
}
});
let damage = damage.map(|Damage::JohnsonCook { d1, d2, d3, d4, d5 }| [*d1, *d2, *d3, *d4, *d5]);
if !(shear > 0.0 && bulk > 0.0 && shear.is_finite() && bulk.is_finite()) {
return Err(Error::Invalid(format!("E = {young}, ν = {poisson} are not a stable material")));
}
Ok(J2 { shear, bulk, flow, damage, heating })
}
/// The temperature a point starts at.
pub fn initial_temperature(&self) -> f64 {
match self.flow {
Some(Flow::JohnsonCook { room, .. }) => room,
_ => 0.0,
}
}
/// Update `state` over an increment of length `dt` with the deviatoric strain increment
/// `dev` (engineering shears; its trace is ignored) and the volumetric strain increment
/// `vol`, at the current density `density` (for heating).
pub fn update(&self, state: &mut PointState, dev: &[f64; 6], vol: f64, dt: f64, density: f64) -> Increment {
let s = &mut state.stress;
let mean = (s[0] + s[1] + s[2]) / 3.0;
let tr = (dev[0] + dev[1] + dev[2]) / 3.0;
let g2 = 2.0 * self.shear;
// Trial deviator.
let mut trial = [s[0] - mean + g2 * (dev[0] - tr), s[1] - mean + g2 * (dev[1] - tr), s[2] - mean + g2 * (dev[2] - tr), 0.0, 0.0, 0.0];
for k in 3..6 {
trial[k] = s[k] + self.shear * dev[k];
}
let pressure = mean + self.bulk * vol;
let mut out = Increment::default();
if let Some(flow) = &self.flow {
let xi: [f64; 6] = std::array::from_fn(|k| trial[k] - state.back[k]);
let norm = norm2(&xi).sqrt();
let q = (1.5f64).sqrt() * norm;
let (yield_now, _) = flow.stress(state.plastic_strain, 0.0, state.temperature, dt);
if q > yield_now && norm > 0.0 {
let hk = flow.kinematic();
let de = self.plastic_increment(flow, state, q, yield_now, hk, dt);
let n: [f64; 6] = xi.map(|x| x / norm);
let radial = g2 * (1.5f64).sqrt() * de;
let shift = (2.0f64 / 3.0).sqrt() * hk * de;
for k in 0..6 {
trial[k] -= radial * n[k];
state.back[k] += shift * n[k];
}
state.plastic_strain += de;
let q_new = q - (3.0 * self.shear + hk) * de;
out = Increment { plastic_strain: de, plastic_work: q_new * de };
if self.heating > 0.0 && density > 0.0 {
state.temperature += self.heating * out.plastic_work / density;
}
}
}
for k in 0..3 {
trial[k] += pressure;
}
state.stress = trial;
if let (Some(d), Some(flow)) = (self.damage, &self.flow)
&& out.plastic_strain > 0.0
{
let q = von_mises(&state.stress);
let triaxiality = if q > 0.0 { pressure / q } else { 0.0 };
let (rate0, homologous) = flow.rate_and_temperature(state.temperature);
let rate = if dt > 0.0 { out.plastic_strain / dt / rate0 } else { 0.0 };
let log_rate = if rate > 1.0 { rate.ln() } else { 0.0 };
let fracture = (d[0] + d[1] * (d[2] * triaxiality).exp()) * (1.0 + d[3] * log_rate) * (1.0 + d[4] * homologous);
state.damage += out.plastic_strain / fracture.max(MIN_FRACTURE_STRAIN);
}
out
}
/// The plastic strain increment of a trial state that yields.
fn plastic_increment(&self, flow: &Flow, state: &PointState, q: f64, yield_now: f64, hk: f64, dt: f64) -> f64 {
let g3 = 3.0 * self.shear;
let (e0, t) = (state.plastic_strain, state.temperature);
match flow {
Flow::Linear { h, kinematic, .. } => (q - yield_now) / (g3 + if *kinematic { hk } else { *h }),
Flow::Table(points) => {
// Walk the segments from the current one: on a segment of slope h starting at
// (e_a, s_a), q − 3G Δ = s_a + h (e0 + Δ − e_a).
for w in points.windows(2) {
let [(ea, sa), (eb, sb)] = *w else { continue };
if eb <= e0 {
continue;
}
let h = (sb - sa) / (eb - ea);
let at = sa + h * (e0 - ea);
let de = (q - at) / (g3 + h);
if e0 + de <= eb {
return de.max(0.0);
}
}
// Beyond the table: perfectly plastic at the last stress.
let last = points.last().map_or(yield_now, |p| p.1);
((q - last) / g3).max(0.0)
}
Flow::JohnsonCook { .. } => {
let residual = |de: f64| {
let rate = if dt > 0.0 { de / dt } else { 0.0 };
let (y, dy) = flow.stress(e0 + de, rate, t, dt);
(q - g3 * de - y, -g3 - dy)
};
let (mut lo, mut hi) = (0.0, ((q - yield_now) / g3).max(0.0));
let mut de = 0.5 * hi;
for _ in 0..MAX_ITERATIONS {
let (r, dr) = residual(de);
if r > 0.0 {
lo = de;
} else {
hi = de;
}
if r.abs() <= 1e-13 * q || hi - lo <= 1e-15 * hi {
break;
}
let newton = de - r / dr;
de = if newton.is_finite() && newton > lo && newton < hi { newton } else { 0.5 * (lo + hi) };
}
de
}
}
}
}
#[cfg(test)]
mod tests {
use super::*;
use modelcraft_model::HardeningPoint;
const E: f64 = 200_000.0;
const NU: f64 = 0.3;
fn material(plastic: Plastic) -> J2 {
J2::from_model(&Elastic::Isotropic { young: E, poisson: NU }, Some(&plastic), None).unwrap()
}
/// Uniaxial stress at a constant strain rate: the axial strain increments by `de` while the
/// lateral strain increment is found by a secant iteration on the point update itself so
/// the lateral stresses vanish. Returns (axial strain, axial stress, plastic strain) per step.
fn uniaxial(m: &J2, steps: &[f64], dt: f64) -> Vec<(f64, f64, f64)> {
let mut s = PointState::new(m.initial_temperature());
let mut out = Vec::new();
let mut strain = 0.0;
for de in steps {
let run = |lateral: f64| {
let mut t = s;
let d = [*de, lateral * de, lateral * de, 0.0, 0.0, 0.0];
m.update(&mut t, &d, d[0] + d[1] + d[2], dt, 7.8e-9);
t
};
let (mut a, mut b) = (-NU, -0.5);
let (mut fa, mut fb) = (run(a).stress[1], run(b).stress[1]);
for _ in 0..100 {
if fb == fa || fb.abs() < 1e-12 {
break;
}
let c = b - fb * (b - a) / (fb - fa);
(a, fa) = (b, fb);
b = c;
fb = run(b).stress[1];
}
s = run(b);
strain += de;
out.push((strain, s.stress[0], s.plastic_strain));
}
out
}
#[test]
fn bilinear_uniaxial_curve_and_unloading() {
let (y0, h) = (250.0, 2000.0);
for kinematic in [false, true] {
let hardening = if kinematic { Hardening::Kinematic } else { Hardening::Isotropic };
let m = material(Plastic::Bilinear { yield_stress: y0, hardening_modulus: h, hardening });
// Load to 1% strain, unload to −1%.
let mut steps = vec![1e-4; 100];
steps.extend(vec![-1e-4; 200]);
let curve = uniaxial(&m, &steps, 1e-3);
// Loading: σ = E ε below yield, then σ = y0 + H ε_p with ε = σ/E + ε_p, that is
// σ = y0 + E_t (ε − y0/E), E_t = E H/(E + H).
let et = E * h / (E + h);
let (e1, s1, p1) = curve[99];
let exact = y0 + et * (e1 - y0 / E);
assert!((s1 - exact).abs() < 1e-9 * exact, "{kinematic}: {s1} vs {exact}");
assert!((p1 - (s1 - y0) / h).abs() < 1e-12, "{p1}");
// Unloading is elastic until the reversed yield: isotropic at −σ1, kinematic at
// σ1 − 2 y0 (the surface keeps its size and moves).
let reverse = if kinematic { s1 - 2.0 * y0 } else { -s1 };
for (e, s, _) in &curve[100..] {
let elastic = s1 - E * (e1 - e);
if elastic > reverse + 1.0 {
assert!((s - elastic).abs() < 1e-8 * s1, "{kinematic}: {s} vs elastic {elastic}");
}
}
let (_, s_end, _) = curve[299];
assert!(s_end < reverse, "{kinematic}: {s_end} vs {reverse}");
eprintln!("bilinear kinematic = {kinematic}: σ(1%) = {s1:.6} (exact {exact:.6}); reversed yield at {reverse:.3}; σ(−1%) = {s_end:.4}");
}
}
#[test]
fn table_and_johnson_cook_follow_their_flow_curves() {
let points = vec![
HardeningPoint { strain: 0.0, stress: 200.0 },
HardeningPoint { strain: 0.01, stress: 300.0 },
HardeningPoint { strain: 0.05, stress: 340.0 },
];
let m = material(Plastic::Table { points });
let curve = uniaxial(&m, &vec![2e-4; 400], 1e-3);
for (_, s, p) in &curve {
if *p > 0.0 {
let want = table(&[(0.0, 200.0), (0.01, 300.0), (0.05, 340.0)], *p).0;
assert!((s - want).abs() < 1e-8 * want, "{s} vs {want} at {p}");
}
}
assert!(curve.last().unwrap().2 > 0.05, "past the table's end");
// Johnson–Cook at a constant plastic rate: the stress is the flow stress at that rate.
let (a, b, n, c) = (350.0, 275.0, 0.36, 0.022);
let jc = Plastic::JohnsonCook {
a,
b,
n,
c,
m: 1.0,
reference_strain_rate: 1.0,
room_temperature: 293.0,
melt_temperature: 1356.0,
specific_heat: None,
taylor_quinney: 0.9,
};
let m = material(jc);
let dt = 1e-6;
let curve = uniaxial(&m, &vec![1e-3; 200], dt);
let mut worst: f64 = 0.0;
for w in curve.windows(2).skip(20) {
let (_, s, p) = w[1];
let rate = (p - w[0].2) / dt;
let want = (a + b * p.powf(n)) * (1.0 + c * rate.ln());
worst = worst.max((s - want).abs() / want);
}
eprintln!("Johnson–Cook uniaxial: largest relative error of the flow stress {worst:.2e}");
assert!(worst < 1e-9, "{worst}");
}
#[test]
fn heating_softens_and_damage_accumulates() {
let jc = |cp: Option<f64>| Plastic::JohnsonCook {
a: 350.0,
b: 275.0,
n: 0.36,
c: 0.0,
m: 1.0,
reference_strain_rate: 1.0,
room_temperature: 293.0,
melt_temperature: 1356.0,
specific_heat: cp,
taylor_quinney: 0.9,
};
let el = Elastic::Isotropic { young: E, poisson: NU };
let damage = Damage::JohnsonCook { d1: 0.5, d2: 0.0, d3: 0.0, d4: 0.0, d5: 0.0 };
let cold = J2::from_model(&el, Some(&jc(None)), Some(&damage)).unwrap();
let hot = J2::from_model(&el, Some(&jc(Some(3.83e8))), Some(&damage)).unwrap();
let a = uniaxial(&cold, &vec![1e-3; 300], 1e-3);
let b = uniaxial(&hot, &vec![1e-3; 300], 1e-3);
assert!(b.last().unwrap().1 < a.last().unwrap().1);
// Constant fracture strain 0.5: D = ε_p / 0.5.
let mut s = PointState::new(293.0);
for _ in 0..300 {
cold.update(&mut s, &[1e-3, -5e-4, -5e-4, 0.0, 0.0, 0.0], 0.0, 1e-3, 7.8e-9);
}
assert!((s.damage - s.plastic_strain / 0.5).abs() < 1e-12, "{} vs {}", s.damage, s.plastic_strain);
}
#[test]
fn hostile_increments_stay_finite() {
let m = material(Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 0.0, hardening: Hardening::Isotropic });
let mut s = PointState::new(0.0);
m.update(&mut s, &[1e3, -1e3, 0.0, 5e2, 0.0, 0.0], 0.0, 0.0, 0.0);
assert!(s.stress.iter().all(|x| x.is_finite()) && (von_mises(&s.stress) - 250.0).abs() < 1e-9 * 250.0);
}
}
crates/fem/src/rotation.rs
+263
View File
@@ -0,0 +1,263 @@
//! Incrementally objective stress updates for large deformation: the rotation that carries a
//! stress (and every other tensor of the material state) from one configuration to the next.
//!
//! Increments are given as the gradient of the displacement increment with respect to the
//! **midpoint configuration**, ΔG = ∂Δu/∂x_{n+½} (= Δt L at the midpoint). Its symmetric part is
//! the strain increment ΔD and its skew part the spin increment ΔW.
//!
//! **Jaumann rate, Hughes–Winget integration** (T. J. R. Hughes, J. Winget, *Finite rotation
//! effects in numerical integration of rate constitutive equations arising in large-deformation
//! analysis*, IJNME 15(12):1862–1867, 1980): the stress is rotated with
//! Q = (I − ½ΔW)⁻¹(I + ½ΔW) and the constitutive increment added: σ_{n+1} = Q σ_n Qᵀ + Δσ(ΔD).
//! Because ΔG is taken on the midpoint configuration, a rigid rotation x_{n+1} = R x_n gives
//! ΔG = 2(R − I)(R + I)⁻¹, which is skew (ΔD = 0), and its Cayley transform Q is R itself: a
//! rigid rotation of any size rotates the stress exactly and strains nothing.
//!
//! **Green–Naghdi rate, Flanagan–Taylor integration** (D. P. Flanagan, L. M. Taylor, *An
//! accurate numerical algorithm for stress integration with finite rotations*, CMAME
//! 62(3):305–320, 1987): the rotation R of the polar decomposition F = V R is integrated
//! incrementally instead of decomposing F. From V̇ = L V − V Ω with Ω = Ṙ Rᵀ and V symmetric,
//! Ω = W + X where the skew X solves V X + X V = D V − V D, that is, with axial vectors,
//! [(tr V) I − V] x = axial(D V − V D). R is advanced with the Cayley transform of Ω Δt, V with
//! V̇, and the material is updated in the unrotated frame: σ̂ = Rᵀ σ R, d̂ = Rᵀ D R. A material
//! whose axes are those of the undeformed body (orthotropic elasticity) is updated correctly in
//! that frame.
//!
//! Tensors are symmetric 6-vectors in the order xx, yy, zz, xy, xz, yz; strain-like vectors carry
//! engineering shears (γ = 2ε), stress-like ones tensor components.
use crate::dense::inverse3;
/// A 3 × 3 matrix, row-major.
pub type Mat3 = [[f64; 3]; 3];
pub const IDENTITY: Mat3 = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
pub fn mul(a: &Mat3, b: &Mat3) -> Mat3 {
let mut c = [[0.0; 3]; 3];
for i in 0..3 {
for j in 0..3 {
c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
}
}
c
}
pub fn transpose(a: &Mat3) -> Mat3 {
[[a[0][0], a[1][0], a[2][0]], [a[0][1], a[1][1], a[2][1]], [a[0][2], a[1][2], a[2][2]]]
}
/// A symmetric tensor from its 6-vector of tensor components.
pub fn tensor(s: &[f64; 6]) -> Mat3 {
[[s[0], s[3], s[4]], [s[3], s[1], s[5]], [s[4], s[5], s[2]]]
}
/// The 6-vector of tensor components of (the symmetric part of) a matrix.
pub fn voigt(a: &Mat3) -> [f64; 6] {
[a[0][0], a[1][1], a[2][2], 0.5 * (a[0][1] + a[1][0]), 0.5 * (a[0][2] + a[2][0]), 0.5 * (a[1][2] + a[2][1])]
}
/// Q s Qᵀ for a stress-like 6-vector.
pub fn rotate(q: &Mat3, s: &[f64; 6]) -> [f64; 6] {
voigt(&mul(&mul(q, &tensor(s)), &transpose(q)))
}
/// Q e Qᵀ for a strain-like 6-vector (engineering shears).
pub fn rotate_strain(q: &Mat3, e: &[f64; 6]) -> [f64; 6] {
let half = [e[0], e[1], e[2], 0.5 * e[3], 0.5 * e[4], 0.5 * e[5]];
let r = rotate(q, &half);
[r[0], r[1], r[2], 2.0 * r[3], 2.0 * r[4], 2.0 * r[5]]
}
/// Q v for a vector.
pub fn apply(q: &Mat3, v: &[f64; 3]) -> [f64; 3] {
[0, 1, 2].map(|i| q[i][0] * v[0] + q[i][1] * v[1] + q[i][2] * v[2])
}
/// The strain increment ΔD (engineering shears) and the spin increment ΔW of a displacement
/// increment gradient.
pub fn split(g: &Mat3) -> ([f64; 6], Mat3) {
let d = [g[0][0], g[1][1], g[2][2], g[0][1] + g[1][0], g[0][2] + g[2][0], g[1][2] + g[2][1]];
let mut w = [[0.0; 3]; 3];
for i in 0..3 {
for j in 0..3 {
w[i][j] = 0.5 * (g[i][j] - g[j][i]);
}
}
(d, w)
}
/// The Cayley transform (I − ½A)⁻¹(I + ½A) of a skew matrix: an exact rotation. `None` only
/// for non-finite input (I − ½A is invertible for every real skew A).
pub fn cayley(a: &Mat3) -> Option<Mat3> {
let mut minus = IDENTITY;
let mut plus = IDENTITY;
for i in 0..3 {
for j in 0..3 {
minus[i][j] -= 0.5 * a[i][j];
plus[i][j] += 0.5 * a[i][j];
}
}
let (inv, _) = inverse3(&minus)?;
let q = mul(&inv, &plus);
q.iter().flatten().all(|v| v.is_finite()).then_some(q)
}
/// The axial vector a of a skew matrix A (A v = a × v).
pub fn axial(a: &Mat3) -> [f64; 3] {
[0.5 * (a[2][1] - a[1][2]), 0.5 * (a[0][2] - a[2][0]), 0.5 * (a[1][0] - a[0][1])]
}
/// The skew matrix of an axial vector.
pub fn skew(a: &[f64; 3]) -> Mat3 {
[[0.0, -a[2], a[1]], [a[2], 0.0, -a[0]], [-a[1], a[0], 0.0]]
}
/// The polar rotation R and left stretch V (a 6-vector of tensor components) of the
/// Flanagan–Taylor integration.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Polar {
pub rotation: Mat3,
pub stretch: [f64; 6],
}
impl Default for Polar {
fn default() -> Polar {
Polar { rotation: IDENTITY, stretch: [1.0, 1.0, 1.0, 0.0, 0.0, 0.0] }
}
}
impl Polar {
/// Advance R and V over an increment with displacement increment gradient `g` (on the
/// midpoint configuration). V̇ = L V − V Ω is integrated with the midpoint rule (a predictor
/// for V at the midpoint, then Ω and the increment from it), so V and R are second-order
/// accurate in the increment. Returns the increment's rotation (R_{n+1} = q R_n); `None`
/// when the increment is not finite or V degenerates.
pub fn advance(&mut self, g: &Mat3) -> Option<Mat3> {
let v0 = tensor(&self.stretch);
let (dv, omega) = Polar::rate(g, &v0)?;
let mut mid = v0;
for i in 0..3 {
for j in 0..3 {
mid[i][j] += 0.5 * dv[i][j];
}
}
let (dv, omega_mid) = Polar::rate(g, &mid).unwrap_or((dv, omega));
let mut vn = v0;
for i in 0..3 {
for j in 0..3 {
vn[i][j] += dv[i][j];
}
}
let q = cayley(&omega_mid)?;
self.rotation = mul(&q, &self.rotation);
self.stretch = voigt(&vn);
self.stretch.iter().all(|s| s.is_finite()).then_some(q)
}
/// The increments of V (L V − V Ω, times Δt) and the polar spin Ω Δt at stretch `v`.
fn rate(g: &Mat3, v: &Mat3) -> Option<(Mat3, Mat3)> {
let (d6, w) = split(g);
let d = tensor(&[d6[0], d6[1], d6[2], 0.5 * d6[3], 0.5 * d6[4], 0.5 * d6[5]]);
let (dv, vd) = (mul(&d, v), mul(v, &d));
let mut rhs = [[0.0; 3]; 3];
for i in 0..3 {
for j in 0..3 {
rhs[i][j] = dv[i][j] - vd[i][j];
}
}
let trace = v[0][0] + v[1][1] + v[2][2];
let mut m = v.map(|r| r.map(|x| -x));
for (i, row) in m.iter_mut().enumerate() {
row[i] += trace;
}
let (inv, _) = inverse3(&m)?;
let xs = skew(&apply(&inv, &axial(&rhs)));
let mut omega = w;
for i in 0..3 {
for j in 0..3 {
omega[i][j] += xs[i][j];
}
}
let (lv, vo) = (mul(g, v), mul(v, &omega));
let mut rate = [[0.0; 3]; 3];
for i in 0..3 {
for j in 0..3 {
rate[i][j] = lv[i][j] - vo[i][j];
}
}
Some((rate, omega))
}
}
#[cfg(test)]
mod tests {
use super::*;
fn rotation_about(axis: [f64; 3], angle: f64) -> Mat3 {
let n = (axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]).sqrt();
let k = skew(&axis.map(|c| c / n));
let k2 = mul(&k, &k);
let mut r = IDENTITY;
for i in 0..3 {
for j in 0..3 {
r[i][j] += angle.sin() * k[i][j] + (1.0 - angle.cos()) * k2[i][j];
}
}
r
}
/// The midpoint gradient of the increment x → R x: 2 (R − I)(R + I)⁻¹.
fn rigid_increment(r: &Mat3) -> Mat3 {
let mut a = *r;
let mut b = *r;
for i in 0..3 {
a[i][i] -= 1.0;
b[i][i] += 1.0;
}
let (inv, _) = inverse3(&b).unwrap();
mul(&a, &inv).map(|row| row.map(|v| 2.0 * v))
}
#[test]
fn a_rigid_rotation_is_exact_with_both_rates() {
let r = rotation_about([0.3, -0.5, 0.8], 1.1);
let g = rigid_increment(&r);
let (d, w) = split(&g);
assert!(d.iter().all(|x| x.abs() < 1e-15), "{d:?}");
let q = cayley(&w).unwrap();
assert!(q.iter().flatten().zip(r.iter().flatten()).all(|(a, b)| (a - b).abs() < 1e-15), "{q:?} vs {r:?}");
let mut p = Polar::default();
p.advance(&g).unwrap();
assert!(p.rotation.iter().flatten().zip(r.iter().flatten()).all(|(a, b)| (a - b).abs() < 1e-14));
assert!(p.stretch.iter().zip([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]).all(|(a, b)| (a - b).abs() < 1e-14), "{:?}", p.stretch);
let s = [1.0, 2.0, 3.0, 0.4, 0.5, 0.6];
let rs = rotate(&q, &s);
// Invariants: trace and the norm.
let norm = |s: &[f64; 6]| s[..3].iter().map(|x| x * x).sum::<f64>() + 2.0 * s[3..].iter().map(|x| x * x).sum::<f64>();
assert!((rs[0] + rs[1] + rs[2] - 6.0).abs() < 1e-14 && (norm(&rs) - norm(&s)).abs() < 1e-13);
}
#[test]
fn polar_rotation_of_a_stretch_then_rotation() {
// Stretch along x by 1.5 in many small steps, then rotate 0.7 rad about z in many small
// steps: V = R diag(1.5, 1, 1) Rᵀ, R = rotation.
let mut p = Polar::default();
let n = 2000;
let step = 1.5f64.powf(1.0 / n as f64);
for _ in 0..n {
// x_{n+1} = diag(step) x_n; midpoint gradient 2 (s − 1)/(s + 1).
let mut g = [[0.0; 3]; 3];
g[0][0] = 2.0 * (step - 1.0) / (step + 1.0);
p.advance(&g).unwrap();
}
let r = rotation_about([0.0, 0.0, 1.0], 0.7 / n as f64);
for _ in 0..n {
p.advance(&rigid_increment(&r)).unwrap();
}
let want = rotation_about([0.0, 0.0, 1.0], 0.7);
let v = rotate(&want, &[1.5, 1.0, 1.0, 0.0, 0.0, 0.0]);
assert!(p.rotation.iter().flatten().zip(want.iter().flatten()).all(|(a, b)| (a - b).abs() < 1e-9), "{:?}", p.rotation);
assert!(p.stretch.iter().zip(v).all(|(a, b)| (a - b).abs() < 1e-6), "{:?} vs {v:?}", p.stretch);
}
}
crates/fem/src/uniform.rs
+105
View File
@@ -109,6 +109,65 @@ impl UniformStrain<8> {
}
}
/// The hexahedron's volume as a form in its corner positions (in the model's hex8 node order):
/// with C_IJK = ∫ ε_abc ∂_aN_I ∂_bN_J ∂_cN_K dξ (the determinant of the Jacobian expanded over
/// the corners), which is antisymmetric in I, J, K, V = Σ_{I<J<K} C_IJK x_I · (x_J × x_K). The 24
/// nonzero coefficients with I < J < K are ±1/12; this table holds them times 12 (the tests
/// integrate them again).
const VOLUME_FORM: [(usize, usize, usize, f64); 24] = [
(0, 1, 2, -1.0),
(0, 1, 3, -1.0),
(0, 1, 4, 1.0),
(0, 1, 5, 1.0),
(0, 2, 3, -1.0),
(0, 3, 4, -1.0),
(0, 3, 7, -1.0),
(0, 4, 5, -1.0),
(0, 4, 7, 1.0),
(1, 2, 3, -1.0),
(1, 2, 5, 1.0),
(1, 2, 6, 1.0),
(1, 4, 5, -1.0),
(1, 5, 6, -1.0),
(2, 3, 6, 1.0),
(2, 3, 7, 1.0),
(2, 5, 6, -1.0),
(2, 6, 7, -1.0),
(3, 4, 7, 1.0),
(3, 6, 7, -1.0),
(4, 5, 6, 1.0),
(4, 5, 7, 1.0),
(4, 6, 7, 1.0),
(5, 6, 7, 1.0),
];
impl UniformStrain<8> {
/// The integrated gradients and volume of a hexahedron in closed form: b_I = ∂V/∂x_I of the
/// volume form ([`VOLUME_FORM`]; D. P. Flanagan, T. Belytschko 1981 give it written out),
/// so b_I gathers C_IJK x_J × x_K, and V = ⅓ Σ_I b_I · x_I. A few hundred operations instead of a Gauss loop, for the
/// current configuration of large deformation. The element is refused when its volume is not
/// positive. `nodal_volume` is left at the mean V/8: an explicit solver keeps the lumped
/// masses of the undeformed element.
pub fn hex8_deformed(xyz: &[[f64; 3]; 8]) -> Result<Hex8> {
let mut b = [[0.0f64; 3]; 8];
let cross = |a: [f64; 3], b: [f64; 3]| [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
for (i, j, k, c) in VOLUME_FORM {
let (xi, xj, xk) = (xyz[i], xyz[j], xyz[k]);
for (target, v) in [(i, cross(xj, xk)), (j, cross(xk, xi)), (k, cross(xi, xj))] {
for d in 0..3 {
b[target][d] += c * v[d];
}
}
}
let b = b.map(|v| v.map(|c| c / 12.0));
let volume: f64 = b.iter().zip(xyz).map(|(bi, x)| bi[0] * x[0] + bi[1] * x[1] + bi[2] * x[2]).sum::<f64>() / 3.0;
if !(volume > 0.0 && volume.is_finite() && b.iter().flatten().all(|v| v.is_finite())) {
return Err(degenerate("hex8"));
}
Ok(UniformStrain { b, volume, nodal_volume: [volume / 8.0; 8] })
}
}
impl UniformStrain<4> {
/// The linear tetrahedron of corners `xyz` (the model's tet4 node order).
pub fn tet4(xyz: &[[f64; 3]; 4]) -> Result<Tet4> {
@@ -281,6 +340,52 @@ mod tests {
assert!(Hex8::hex8(&[[f64::NAN; 3]; 8]).is_err());
}
#[test]
fn the_volume_form_is_the_integrated_one() {
let mut c = [[[0.0f64; 8]; 8]; 8];
for p in Rule::Hex(2).points() {
let s = shape::hex8(p.xi);
let d = s.derivatives();
for i in 0..8 {
for j in 0..8 {
for k in 0..8 {
let m = [[d[i][0], d[j][0], d[k][0]], [d[i][1], d[j][1], d[k][1]], [d[i][2], d[j][2], d[k][2]]];
c[i][j][k] += p.weight * crate::dense::det3(&m);
}
}
}
}
let mut table = Vec::new();
for i in 0..8 {
for j in i + 1..8 {
for k in j + 1..8 {
if (12.0 * c[i][j][k]).abs() > 1e-9 {
assert!((12.0 * c[i][j][k] - (12.0 * c[i][j][k]).round()).abs() < 1e-12);
table.push((i, j, k, (12.0 * c[i][j][k]).round()));
}
}
}
}
assert_eq!(table, VOLUME_FORM.to_vec());
}
#[test]
fn closed_form_gradients_equal_the_integrated_ones() {
let hex = distorted();
let (a, b) = (Hex8::hex8(&hex).unwrap(), Hex8::hex8_deformed(&hex).unwrap());
assert!((a.volume - b.volume).abs() < 1e-14 * a.volume, "{} vs {}", a.volume, b.volume);
for (x, y) in a.b.iter().flatten().zip(b.b.iter().flatten()) {
assert!((x - y).abs() < 1e-14, "{x} vs {y}");
}
let mut inverted = hex;
inverted.swap(0, 1);
inverted.swap(2, 3);
inverted.swap(4, 5);
inverted.swap(6, 7);
assert!(Hex8::hex8_deformed(&inverted).is_err());
assert!(Hex8::hex8_deformed(&[[f64::NAN; 3]; 8]).is_err());
}
#[test]
fn hourglass_vectors_are_orthogonal_to_linear_fields() {
let hex = distorted();