sparkl2d_core/dynamics/models/
plasticity_rankine.rs

1use crate::math::{Matrix, Real, DIM};
2use core::cmp::Ordering;
3
4#[cfg(not(feature = "std"))]
5use na::ComplexField;
6
7#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
8#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
9#[derive(Copy, Clone, Debug, PartialEq)]
10#[repr(C)]
11pub struct RankinePlasticity {
12    mu: Real,
13    lambda: Real,
14    tensile_strength: Real,
15    softening_rate: Real,
16}
17
18impl RankinePlasticity {
19    pub fn new(
20        young_modulus: Real,
21        poisson_ratio: Real,
22        tensile_strength: Real,
23        softening_rate: Real,
24    ) -> Self {
25        let (lambda, mu) = crate::utils::lame_lambda_mu(young_modulus, poisson_ratio);
26
27        Self {
28            lambda,
29            mu,
30            tensile_strength,
31            softening_rate,
32        }
33    }
34
35    pub fn update_particle(
36        &self,
37        particle_deformation_gradient: &mut Matrix<Real>,
38        particle_plastic_hardening: &mut Real,
39    ) {
40        let lambda = self.lambda;
41        let mu = self.mu;
42
43        let mut svd = particle_deformation_gradient.svd_unordered(true, true);
44        let mut eigv = svd.singular_values.map(|s| s.ln()); // Hencky strain eigenvalues.
45
46        let prev_eigv = eigv;
47        let mut idx = [0usize, 1, DIM - 1]; // We use DIM - 1 so that it works in 2D too.
48        idx.sort_unstable_by(|a, b| eigv[*a].partial_cmp(&eigv[*b]).unwrap_or(Ordering::Equal));
49        let [e3, e2, e1] = idx;
50        // TODO: The -1 is to account for the initial hardening value equal to 1.0.
51        //       Having to do this sounds error-prone (what if the user modifies that initial value?
52        //       Should we let the user modify that initial value?
53        let soft_tensile_strength = self.tensile_strength - (*particle_plastic_hardening - 1.0);
54
55        if lambda * eigv.sum() + 2.0 * mu * eigv[e1] <= soft_tensile_strength {
56            // No plastic deformation.
57            return;
58        } else if (2.0 * mu + lambda) * eigv[e2] + lambda * (eigv.sum() - eigv[e1])
59            <= soft_tensile_strength
60        {
61            let new_eigv =
62                (soft_tensile_strength - lambda * (eigv.sum() - eigv[e1])) / (2.0 * mu + lambda);
63            eigv[e1] = new_eigv;
64        } else if DIM == 3 && (2.0 * mu + 3.0 * lambda) * eigv[e3] <= soft_tensile_strength {
65            let new_eigv = (soft_tensile_strength - lambda * (eigv.sum() - eigv[e1] - eigv[e2]))
66                / (2.0 * mu + 2.0 * lambda);
67            eigv[e1] = new_eigv;
68            eigv[e2] = new_eigv;
69        } else {
70            let new_eigv = soft_tensile_strength / (2.0 * self.mu + 3.0 * self.lambda);
71            eigv.fill(new_eigv);
72        }
73
74        *particle_plastic_hardening += self.softening_rate * (prev_eigv - eigv).norm();
75        *particle_plastic_hardening = particle_plastic_hardening.min(self.tensile_strength);
76        svd.singular_values = eigv.map(|e| e.exp());
77        *particle_deformation_gradient = svd.recompose().unwrap();
78    }
79}
sparkl2d_core/dynamics/models/plasticity_rankine.rs

sparkl2d_core/dynamics/models/
plasticity_rankine.rs
