carapace_core/model/materials/
plane.rs

1//! `PlaneMaterial`: the n-D material seam for 2D continuum elements.
2//!
3//! Same trial/commit design as the uniaxial `Material` (see its doc
4//! comment): a value holds only *committed* state, `trial_stress_tangent` is
5//! pure, and `commit` returns the new committed value. An element owns one
6//! independent `PlaneMaterial` per Gauss point, cloned from the prototype
7//! given to its constructor and committed individually, exactly as the fiber
8//! elements own one material per fiber. That is the seam where nonlinear
9//! plane materials (CSFM, damage, plasticity) plug in with per-point history.
10//!
11//! Conventions, fixed now so nonlinear materials later do not change them:
12//! the strain vector is `[eps_x, eps_y, gamma_xy]` with **engineering
13//! shear** `gamma_xy = 2 eps_xy`; the stress vector is `[sigma_x, sigma_y,
14//! tau_xy]`; `D` maps one to the other.
15
16use std::fmt;
17
18use nalgebra::{SMatrix, SVector};
19
20/// A plane strain or stress vector `[x, y, xy]` (engineering shear for strain).
21pub type PlaneVector = SVector<f64, 3>;
22/// The `3x3` plane constitutive tangent.
23pub type PlaneMatrix = SMatrix<f64, 3, 3>;
24
25#[derive(Debug, Clone, PartialEq)]
26pub enum PlaneMaterial {
27    /// A general symmetric `D`, given as the packed upper triangle row by
28    /// row: `[d11, d12, d13, d22, d23, d33]`, with respect to the
29    /// engineering-shear strain vector.
30    ElasticMatrix { d: [f64; 6] },
31    /// Orthotropic plane stress. `angle` is the counter-clockwise angle in
32    /// radians from global x to material axis 1; `nu_xy` is the major
33    /// Poisson ratio (strain in axis 2 per strain in axis 1, loading along 1).
34    Orthotropic {
35        ex: f64,
36        ey: f64,
37        nu_xy: f64,
38        g_xy: f64,
39        angle: f64,
40    },
41    /// Isotropic plane stress (`sigma_z = 0`).
42    PlaneStress { e: f64, nu: f64 },
43    /// Isotropic plane strain (`eps_z = 0`).
44    PlaneStrain { e: f64, nu: f64 },
45    /// Test-only stateful material: `D` of unit-modulus plane stress
46    /// (`nu = 0`) acting on the strain measured from the last committed
47    /// strain, so per-point history is observable.
48    #[cfg(test)]
49    Probe { last: [f64; 3] },
50}
51
52/// Why a `PlaneMaterial`'s constants are unusable.
53#[derive(Debug, Clone, Copy, PartialEq)]
54pub enum PlaneMaterialError {
55    /// A modulus is not finite and positive.
56    InvalidModulus,
57    /// A Poisson ratio is outside the range for which `D` is positive definite.
58    InvalidPoissonRatio,
59    /// The constitutive matrix is not finite, symmetric positive definite.
60    NotPositiveDefinite,
61}
62
63impl PlaneMaterialError {
64    pub fn message(&self) -> &'static str {
65        match self {
66            PlaneMaterialError::InvalidModulus => "a modulus is not finite and positive",
67            PlaneMaterialError::InvalidPoissonRatio => "Poisson ratio is out of range",
68            PlaneMaterialError::NotPositiveDefinite => {
69                "the constitutive matrix is not symmetric positive definite"
70            }
71        }
72    }
73}
74
75impl fmt::Display for PlaneMaterialError {
76    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
77        f.write_str(self.message())
78    }
79}
80
81impl std::error::Error for PlaneMaterialError {}
82
83impl PlaneMaterial {
84    /// Isotropic plane stress; `e > 0` and `-1 < nu < 1`.
85    pub fn plane_stress(e: f64, nu: f64) -> Result<Self, PlaneMaterialError> {
86        Self::checked(PlaneMaterial::PlaneStress { e, nu })
87    }
88
89    /// Isotropic plane strain; `e > 0` and `-1 < nu < 0.5`.
90    pub fn plane_strain(e: f64, nu: f64) -> Result<Self, PlaneMaterialError> {
91        Self::checked(PlaneMaterial::PlaneStrain { e, nu })
92    }
93
94    /// Orthotropic plane stress; positive moduli and `nu_xy^2 < ex / ey`.
95    pub fn orthotropic(
96        ex: f64,
97        ey: f64,
98        nu_xy: f64,
99        g_xy: f64,
100        angle: f64,
101    ) -> Result<Self, PlaneMaterialError> {
102        Self::checked(PlaneMaterial::Orthotropic {
103            ex,
104            ey,
105            nu_xy,
106            g_xy,
107            angle,
108        })
109    }
110
111    /// A general symmetric positive-definite `D` (packed upper triangle).
112    pub fn elastic_matrix(d: [f64; 6]) -> Result<Self, PlaneMaterialError> {
113        Self::checked(PlaneMaterial::ElasticMatrix { d })
114    }
115
116    fn checked(material: PlaneMaterial) -> Result<Self, PlaneMaterialError> {
117        material.validate()?;
118        Ok(material)
119    }
120
121    /// Checks the constants of a value built from the variants directly.
122    pub fn validate(&self) -> Result<(), PlaneMaterialError> {
123        let positive = |v: f64| {
124            if v.is_finite() && v > 0.0 {
125                Ok(())
126            } else {
127                Err(PlaneMaterialError::InvalidModulus)
128            }
129        };
130        let poisson = |nu: f64, upper: f64| {
131            if nu.is_finite() && nu > -1.0 && nu < upper {
132                Ok(())
133            } else {
134                Err(PlaneMaterialError::InvalidPoissonRatio)
135            }
136        };
137        match *self {
138            PlaneMaterial::PlaneStress { e, nu } => {
139                positive(e)?;
140                poisson(nu, 1.0)?;
141            }
142            PlaneMaterial::PlaneStrain { e, nu } => {
143                positive(e)?;
144                poisson(nu, 0.5)?;
145            }
146            PlaneMaterial::Orthotropic {
147                ex,
148                ey,
149                nu_xy,
150                g_xy,
151                angle,
152            } => {
153                positive(ex)?;
154                positive(ey)?;
155                positive(g_xy)?;
156                if !nu_xy.is_finite() || !angle.is_finite() || nu_xy * nu_xy >= ex / ey {
157                    return Err(PlaneMaterialError::InvalidPoissonRatio);
158                }
159            }
160            PlaneMaterial::ElasticMatrix { d } => {
161                let m = self.d_matrix();
162                if d.iter().any(|v| !v.is_finite()) || m.cholesky().is_none() {
163                    return Err(PlaneMaterialError::NotPositiveDefinite);
164                }
165            }
166            #[cfg(test)]
167            PlaneMaterial::Probe { .. } => {}
168        }
169        Ok(())
170    }
171
172    /// The constitutive matrix at the committed state (constant for every
173    /// elastic variant).
174    pub fn d_matrix(&self) -> PlaneMatrix {
175        match *self {
176            PlaneMaterial::ElasticMatrix { d } => PlaneMatrix::new(
177                d[0], d[1], d[2], //
178                d[1], d[3], d[4], //
179                d[2], d[4], d[5],
180            ),
181            PlaneMaterial::PlaneStress { e, nu } => {
182                let c = e / (1.0 - nu * nu);
183                PlaneMatrix::new(
184                    c,
185                    c * nu,
186                    0.0, //
187                    c * nu,
188                    c,
189                    0.0, //
190                    0.0,
191                    0.0,
192                    c * (1.0 - nu) / 2.0,
193                )
194            }
195            PlaneMaterial::PlaneStrain { e, nu } => {
196                let c = e / ((1.0 + nu) * (1.0 - 2.0 * nu));
197                PlaneMatrix::new(
198                    c * (1.0 - nu),
199                    c * nu,
200                    0.0, //
201                    c * nu,
202                    c * (1.0 - nu),
203                    0.0, //
204                    0.0,
205                    0.0,
206                    c * (1.0 - 2.0 * nu) / 2.0,
207                )
208            }
209            PlaneMaterial::Orthotropic {
210                ex,
211                ey,
212                nu_xy,
213                g_xy,
214                angle,
215            } => {
216                let nu_yx = nu_xy * ey / ex;
217                let k = 1.0 / (1.0 - nu_xy * nu_yx);
218                let local = PlaneMatrix::new(
219                    k * ex,
220                    k * nu_xy * ey,
221                    0.0, //
222                    k * nu_xy * ey,
223                    k * ey,
224                    0.0, //
225                    0.0,
226                    0.0,
227                    g_xy,
228                );
229                // Strain transformation (global -> material axes), engineering shear;
230                // stress transforms by its inverse transpose, so D = T^T D' T.
231                let (s, c) = angle.sin_cos();
232                let t = PlaneMatrix::new(
233                    c * c,
234                    s * s,
235                    c * s, //
236                    s * s,
237                    c * c,
238                    -c * s, //
239                    -2.0 * c * s,
240                    2.0 * c * s,
241                    c * c - s * s,
242                );
243                t.transpose() * local * t
244            }
245            #[cfg(test)]
246            PlaneMaterial::Probe { .. } => PlaneMatrix::identity(),
247        }
248    }
249
250    /// Stress and tangent at `strain`, always relative to the committed state.
251    pub fn trial_stress_tangent(&self, strain: &PlaneVector) -> (PlaneVector, PlaneMatrix) {
252        let d = self.d_matrix();
253        #[cfg(test)]
254        if let PlaneMaterial::Probe { last } = *self {
255            return (d * (strain - PlaneVector::from(last)), d);
256        }
257        (d * strain, d)
258    }
259
260    /// The new committed material after `strain` converged.
261    pub fn commit(&self, strain: &PlaneVector) -> Self {
262        #[cfg(test)]
263        if let PlaneMaterial::Probe { .. } = self {
264            return PlaneMaterial::Probe {
265                last: [strain[0], strain[1], strain[2]],
266            };
267        }
268        let _ = strain;
269        self.clone()
270    }
271
272    /// Whether the response is linear in strain and history-free (what the
273    /// enhanced `Quad4` requires).
274    pub fn is_linear(&self) -> bool {
275        match self {
276            PlaneMaterial::ElasticMatrix { .. }
277            | PlaneMaterial::Orthotropic { .. }
278            | PlaneMaterial::PlaneStress { .. }
279            | PlaneMaterial::PlaneStrain { .. } => true,
280            #[cfg(test)]
281            PlaneMaterial::Probe { .. } => false,
282        }
283    }
284}
285
286#[cfg(test)]
287mod tests {
288    use super::*;
289
290    fn close(a: &PlaneMatrix, b: &PlaneMatrix, tol: f64) -> bool {
291        (a - b).abs().max() <= tol * a.abs().max().max(1.0)
292    }
293
294    #[test]
295    fn d_is_symmetric_and_positive_definite_for_every_variant() {
296        let materials = [
297            PlaneMaterial::plane_stress(200e3, 0.3).unwrap(),
298            PlaneMaterial::plane_strain(200e3, 0.3).unwrap(),
299            PlaneMaterial::orthotropic(30e3, 10e3, 0.25, 4e3, 0.6).unwrap(),
300            PlaneMaterial::elastic_matrix([4.0, 1.0, 0.5, 3.0, 0.2, 2.0]).unwrap(),
301        ];
302        for m in materials {
303            let d = m.d_matrix();
304            assert!(close(&d, &d.transpose(), 1e-14), "{m:?}");
305            assert!(d.cholesky().is_some(), "{m:?}");
306        }
307    }
308
309    #[test]
310    fn plane_stress_with_effective_constants_equals_plane_strain() {
311        let (e, nu) = (200e3, 0.3);
312        let strain = PlaneMaterial::plane_strain(e, nu).unwrap();
313        let (e_eff, nu_eff) = (e / (1.0 - nu * nu), nu / (1.0 - nu));
314        let stress = PlaneMaterial::plane_stress(e_eff, nu_eff).unwrap();
315        assert!(close(&strain.d_matrix(), &stress.d_matrix(), 1e-13));
316    }
317
318    #[test]
319    fn orthotropic_with_equal_moduli_is_isotropic_for_any_angle() {
320        let (e, nu) = (100.0, 0.25);
321        let iso = PlaneMaterial::plane_stress(e, nu).unwrap().d_matrix();
322        for angle in [0.0, 0.3, 1.1, -2.0] {
323            let ortho =
324                PlaneMaterial::orthotropic(e, e, nu, e / (2.0 * (1.0 + nu)), angle).unwrap();
325            assert!(close(&ortho.d_matrix(), &iso, 1e-13), "angle {angle}");
326        }
327    }
328
329    #[test]
330    fn rotating_by_ninety_degrees_swaps_the_axes() {
331        let (ex, ey, nu_xy, g) = (30.0, 10.0, 0.2, 4.0);
332        let d0 = PlaneMaterial::orthotropic(ex, ey, nu_xy, g, 0.0)
333            .unwrap()
334            .d_matrix();
335        let d90 = PlaneMaterial::orthotropic(ex, ey, nu_xy, g, std::f64::consts::FRAC_PI_2)
336            .unwrap()
337            .d_matrix();
338        assert!((d0[(0, 0)] - d90[(1, 1)]).abs() < 1e-12);
339        assert!((d0[(1, 1)] - d90[(0, 0)]).abs() < 1e-12);
340        assert!((d0[(0, 1)] - d90[(0, 1)]).abs() < 1e-12);
341        assert!((d0[(2, 2)] - d90[(2, 2)]).abs() < 1e-12);
342        assert!(d90[(0, 2)].abs() < 1e-12 && d90[(1, 2)].abs() < 1e-12);
343    }
344
345    #[test]
346    fn rotation_by_pi_over_four_couples_normal_and_shear() {
347        // Uniaxial strain along the fibre direction at 45 degrees: sigma_x - sigma_y is zero by symmetry,
348        // and the rotated D has a nonzero normal-shear coupling for a strongly orthotropic ply.
349        let d = PlaneMaterial::orthotropic(100.0, 5.0, 0.3, 4.0, std::f64::consts::FRAC_PI_4)
350            .unwrap()
351            .d_matrix();
352        assert!((d[(0, 0)] - d[(1, 1)]).abs() < 1e-12);
353        assert!(d[(0, 2)].abs() > 1.0);
354    }
355
356    #[test]
357    fn elastic_matrix_packing_matches_the_equivalent_isotropic_d() {
358        let iso = PlaneMaterial::plane_stress(100.0, 0.25).unwrap().d_matrix();
359        let packed = [
360            iso[(0, 0)],
361            iso[(0, 1)],
362            iso[(0, 2)],
363            iso[(1, 1)],
364            iso[(1, 2)],
365            iso[(2, 2)],
366        ];
367        let m = PlaneMaterial::elastic_matrix(packed).unwrap();
368        assert!(close(&m.d_matrix(), &iso, 1e-15));
369        let strain = PlaneVector::new(1e-3, -2e-3, 3e-3);
370        assert!((m.trial_stress_tangent(&strain).0 - iso * strain).norm() < 1e-12);
371    }
372
373    #[test]
374    fn invalid_constants_are_errors_not_panics() {
375        assert_eq!(
376            PlaneMaterial::plane_stress(0.0, 0.3),
377            Err(PlaneMaterialError::InvalidModulus)
378        );
379        assert_eq!(
380            PlaneMaterial::plane_stress(-1.0, 0.3),
381            Err(PlaneMaterialError::InvalidModulus)
382        );
383        assert_eq!(
384            PlaneMaterial::plane_stress(f64::NAN, 0.3),
385            Err(PlaneMaterialError::InvalidModulus)
386        );
387        assert_eq!(
388            PlaneMaterial::plane_stress(1.0, 1.0),
389            Err(PlaneMaterialError::InvalidPoissonRatio)
390        );
391        assert_eq!(
392            PlaneMaterial::plane_strain(1.0, 0.5),
393            Err(PlaneMaterialError::InvalidPoissonRatio)
394        );
395        assert_eq!(
396            PlaneMaterial::plane_strain(1.0, -1.0),
397            Err(PlaneMaterialError::InvalidPoissonRatio)
398        );
399        assert_eq!(
400            PlaneMaterial::orthotropic(1.0, 10.0, 0.5, 1.0, 0.0),
401            Err(PlaneMaterialError::InvalidPoissonRatio)
402        );
403        assert_eq!(
404            PlaneMaterial::orthotropic(1.0, 1.0, 0.0, 0.0, 0.0),
405            Err(PlaneMaterialError::InvalidModulus)
406        );
407        assert_eq!(
408            PlaneMaterial::elastic_matrix([1.0, 2.0, 0.0, 1.0, 0.0, 1.0]),
409            Err(PlaneMaterialError::NotPositiveDefinite)
410        );
411        assert_eq!(
412            PlaneMaterial::elastic_matrix([f64::NAN, 0.0, 0.0, 1.0, 0.0, 1.0]),
413            Err(PlaneMaterialError::NotPositiveDefinite)
414        );
415        // A struct literal that bypasses the constructors is caught by `validate`.
416        assert!(PlaneMaterial::PlaneStress { e: -1.0, nu: 0.0 }
417            .validate()
418            .is_err());
419    }
420
421    #[test]
422    fn elastic_variants_are_linear_and_commit_changes_nothing() {
423        let m = PlaneMaterial::plane_stress(10.0, 0.2).unwrap();
424        assert!(m.is_linear());
425        assert_eq!(m.commit(&PlaneVector::new(1.0, 2.0, 3.0)), m);
426    }
427
428    #[test]
429    fn a_stateful_material_keeps_independent_history_per_copy() {
430        let proto = PlaneMaterial::Probe { last: [0.0; 3] };
431        let (mut a, mut b) = (proto.clone(), proto);
432        a = a.commit(&PlaneVector::new(1.0, 0.0, 0.0));
433        b = b.commit(&PlaneVector::new(0.0, 2.0, 0.0));
434        let probe = PlaneVector::new(1.0, 2.0, 0.0);
435        assert_eq!(
436            a.trial_stress_tangent(&probe).0,
437            PlaneVector::new(0.0, 2.0, 0.0)
438        );
439        assert_eq!(
440            b.trial_stress_tangent(&probe).0,
441            PlaneVector::new(1.0, 0.0, 0.0)
442        );
443        assert!(!a.is_linear());
444    }
445}
carapace_core/model/materials/plane.rs

carapace_core/model/materials/
plane.rs
