sparkl2d_core/dynamics/models/
elasticity_corotated_linear.rs

1use crate::dynamics::models::ActiveTimestepBounds;
2use crate::dynamics::timestep::ElasticitySoundSpeedTimestepBound;
3use crate::math::{Matrix, Real, Vector};
4
5#[cfg(not(feature = "std"))]
6use na::ComplexField;
7
8#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
9#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
10#[derive(Copy, Clone, Debug, PartialEq)]
11#[repr(C)]
12pub struct CorotatedLinearElasticity {
13    pub split_stress_on_failure: bool,
14    pub cfl_coeff: Real,
15    pub lambda: Real,
16    pub mu: Real,
17}
18
19impl CorotatedLinearElasticity {
20    pub fn new(young_modulus: Real, poisson_ratio: Real) -> Self {
21        let (lambda, mu) = crate::utils::lame_lambda_mu(young_modulus, poisson_ratio);
22
23        Self {
24            split_stress_on_failure: true,
25            cfl_coeff: 0.9,
26            lambda,
27            mu,
28        }
29    }
30
31    pub fn kirchhoff_stress(
32        &self,
33        phase: Real,
34        hardening: Real,
35        deformation_gradient: &Matrix<Real>,
36    ) -> Matrix<Real> {
37        let j = deformation_gradient.determinant();
38        let mut svd = deformation_gradient.svd_unordered(true, true);
39        svd.singular_values.apply(|e| *e = *e - 1.0);
40
41        if phase == 1.0 {
42            (2.0 * self.mu * hardening)
43                * svd.recompose().unwrap()
44                * deformation_gradient.transpose()
45                + (self.lambda * hardening * (j - 1.0) * j) * Matrix::identity()
46        } else {
47            let mut pos_def = svd;
48            pos_def.singular_values.apply(|e| *e = e.max(0.0));
49            let mut neg_def = svd;
50            neg_def.singular_values.apply(|e| *e = e.min(0.0));
51
52            let pos_dev_stress = 2.0
53                * self.mu
54                * hardening
55                * pos_def.recompose().unwrap()
56                * deformation_gradient.transpose();
57            let neg_dev_stress = 2.0
58                * self.mu
59                * hardening
60                * neg_def.recompose().unwrap()
61                * deformation_gradient.transpose();
62
63            let spherical_stress = (self.lambda * hardening * (j - 1.0) * j) * Matrix::identity();
64
65            let (pos_part, neg_part) = if j < 1.0 {
66                (pos_dev_stress, neg_dev_stress + spherical_stress)
67            } else {
68                (pos_dev_stress + spherical_stress, neg_dev_stress)
69            };
70
71            let phase_coeff = if self.split_stress_on_failure && phase == 0.0 {
72                0.0
73            } else {
74                1.0
75            };
76            pos_part * phase_coeff + neg_part
77        }
78    }
79
80    // General elastic density function: https://www.math.ucla.edu/~cffjiang/research/mpmcourse/mpmcourse.pdf#subsection.6.3 (49)
81    // With hardening: https://www.math.ucla.edu/~cffjiang/research/mpmcourse/mpmcourse.pdf#subsection.6.5 (87)
82    pub fn elastic_energy_density(
83        &self,
84        deformation_gradient: Matrix<Real>,
85        elastic_hardening: Real,
86    ) -> Real {
87        let singular_values = deformation_gradient
88            .svd_unordered(false, false)
89            .singular_values;
90        let determinant: Real = singular_values.iter().product();
91
92        let hardened_mu = self.mu * elastic_hardening;
93        let hardened_lambda = self.lambda * elastic_hardening;
94
95        hardened_mu
96            * singular_values
97                .iter()
98                .map(|sigma| (sigma - 1.).powi(2))
99                .sum::<Real>()
100            + hardened_lambda / 2. * (determinant - 1.).powi(2)
101    }
102
103    // Basically the same as [elastic_energy_density] but only the positive part
104    // See "Dynamic brittle fracture with eigenerosion enhanced material point method", 2.2 Elasticity degradation
105    pub fn pos_energy(&self, deformation_gradient: Matrix<Real>, elastic_hardening: Real) -> Real {
106        let j = deformation_gradient.determinant();
107        let mut pos_def = deformation_gradient.svd_unordered(false, false);
108        pos_def.singular_values.apply(|e| *e = (*e - 1.0).max(0.0));
109
110        let pos_dev_part = self.mu * elastic_hardening * pos_def.singular_values.norm_squared();
111        let spherical_part = self.lambda * elastic_hardening / 2.0 * (j - 1.0) * (j - 1.0);
112
113        if j < 1.0 {
114            pos_dev_part
115        } else {
116            pos_dev_part + spherical_part
117        }
118    }
119
120    pub fn active_timestep_bounds(&self) -> ActiveTimestepBounds {
121        ActiveTimestepBounds::CONSTITUTIVE_MODEL_BOUND
122            | ActiveTimestepBounds::PARTICLE_VELOCITY_BOUND
123            | ActiveTimestepBounds::DEFORMATION_GRADIENT_CHANGE_BOUND
124    }
125
126    pub fn timestep_bound(
127        &self,
128        particle_density0: Real,
129        particle_velocity: &Vector<Real>,
130        elastic_hardening: Real,
131        cell_width: Real,
132    ) -> Real {
133        let bulk_modulus = crate::utils::bulk_modulus_from_lame(self.lambda, self.mu);
134        let shear_modulus = crate::utils::shear_modulus_from_lame(self.lambda, self.mu);
135
136        let bound = ElasticitySoundSpeedTimestepBound {
137            alpha: self.cfl_coeff,
138            bulk_modulus: bulk_modulus * elastic_hardening,
139            shear_modulus: shear_modulus * elastic_hardening,
140        };
141        bound.timestep_bound(particle_density0, particle_velocity, cell_width)
142    }
143
144    pub fn is_fluid(&self) -> bool {
145        false
146    }
147}
sparkl2d_core/dynamics/models/elasticity_corotated_linear.rs

sparkl2d_core/dynamics/models/
elasticity_corotated_linear.rs
