sparkl3d_core/dynamics/models/
elasticity_neo_hookean.rs

1use crate::dynamics::models::ActiveTimestepBounds;
2use crate::dynamics::timestep::ElasticitySoundSpeedTimestepBound;
3use crate::math::{Matrix, Real, Vector, DIM};
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 NeoHookeanElasticity {
12    pub cfl_coeff: Real,
13    pub lambda: Real,
14    pub mu: Real,
15}
16
17impl NeoHookeanElasticity {
18    pub fn new(young_modulus: Real, poisson_ratio: Real) -> Self {
19        let (lambda, mu) = crate::utils::lame_lambda_mu(young_modulus, poisson_ratio);
20
21        Self {
22            cfl_coeff: 0.5,
23            lambda,
24            mu,
25        }
26    }
27
28    pub fn fpk_stress(&self, deformation_gradient: &Matrix<Real>, hardening: Real) -> Matrix<Real> {
29        let j = deformation_gradient.determinant();
30        let k = 2.0 / 3.0 * self.mu * hardening + self.lambda * hardening;
31        let inv_f_tr = deformation_gradient
32            .transpose()
33            .try_inverse()
34            .unwrap_or_else(Matrix::identity);
35        let cauchy_green_strain = deformation_gradient * deformation_gradient.transpose();
36        let deviatoric_stress = (self.mu * hardening * j.powf(-2.0 / (DIM as Real)))
37            * crate::utils::deviatoric_part(&cauchy_green_strain)
38            * inv_f_tr;
39        let volumetric_stress = j * k / 2.0 * (j - 1.0 / j) * inv_f_tr;
40        deviatoric_stress + volumetric_stress
41    }
42
43    pub fn kirchhoff_stress(
44        &self,
45        particle_phase: Real,
46        particle_elastic_hardening: Real,
47        particle_deformation_gradient: &Matrix<Real>,
48    ) -> Matrix<Real> {
49        let particle_phase_coeff = Self::phase_coeff(particle_phase);
50        let j = particle_deformation_gradient.determinant();
51        let k = 2.0 / 3.0 * self.mu * particle_elastic_hardening
52            + self.lambda * particle_elastic_hardening;
53        let cauchy_green_strain =
54            particle_deformation_gradient * particle_deformation_gradient.transpose();
55        let deviatoric_stress = (self.mu * particle_elastic_hardening)
56            * j.powf(-2.0 / (DIM as Real))
57            * crate::utils::deviatoric_part(&cauchy_green_strain);
58        let volumetric_stress = k / 2.0 * (j * j - 1.0) * Matrix::identity();
59
60        let (pos_part, neg_part) = if j >= 1.0 {
61            (deviatoric_stress + volumetric_stress, Matrix::zeros())
62        } else {
63            (deviatoric_stress, volumetric_stress)
64        };
65
66        pos_part * particle_phase_coeff + neg_part
67    }
68
69    // https://dl.acm.org/doi/pdf/10.1145/3306346.3322949#subsection.3.2
70    pub fn phase_coeff(phase: Real) -> Real {
71        const R: Real = 0.001;
72        (1.0 - R) * phase * phase + R
73    }
74
75    pub fn is_fluid(&self) -> bool {
76        false
77    }
78
79    // General elastic density function: https://dl.acm.org/doi/pdf/10.1145/3306346.3322949#section.5 (8) and (9)
80    // With degradation: https://dl.acm.org/doi/pdf/10.1145/3306346.3322949#subsection.3.2 (3)
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        particle_phase: Real,                // aka. c
85        elastic_hardening: Real,             // aka. e^(xi (1 - J_P))
86        deformation_gradient: &Matrix<Real>, // aka. F
87    ) -> Real {
88        let hardened_mu = self.mu * elastic_hardening;
89        let hardened_lambda = self.lambda * elastic_hardening;
90
91        // aka. bulk modulus
92        let kappa = 2.0 / 3.0 * hardened_mu + hardened_lambda;
93
94        let j = deformation_gradient.determinant();
95        let alpha = -1. / DIM as Real;
96
97        // aka. Psi_mu
98        let deviatoric_part = hardened_mu / 2.
99            * (((j.powf(alpha) * deformation_gradient).transpose()
100                * (j.powf(alpha) * deformation_gradient))
101                .trace()
102                - DIM as Real);
103
104        // aka. Psi_kappa
105        let volumetric_part = kappa / 2. * ((j.powi(2) - 1.) / 2. - j.ln());
106
107        // aka. (Psi_plus, Psi_minus)
108        let (tensile_part, compressive_part) = if j >= 1. {
109            (volumetric_part + deviatoric_part, 0.)
110        } else {
111            (volumetric_part, deviatoric_part)
112        };
113
114        Self::phase_coeff(particle_phase) * tensile_part + compressive_part
115    }
116
117    // Basically the same as [elastic_energy_density] but only the positive part
118    pub fn pos_energy(
119        &self,
120        particle_phase: Real,
121        particle_elastic_hardening: Real,
122        particle_deformation_gradient: &Matrix<Real>,
123    ) -> Real {
124        let phase_coeff = Self::phase_coeff(particle_phase);
125        let hardening = particle_elastic_hardening;
126
127        let j = particle_deformation_gradient.determinant();
128        let k = 2.0 / 3.0 * self.mu * hardening + self.lambda * hardening;
129        let deviatoric_part = hardening * self.mu / 2.0
130            * ((particle_deformation_gradient * particle_deformation_gradient.transpose()).trace()
131                * j.powf(-2.0 / (DIM as Real))
132                - DIM as Real);
133        let volumetric_part = k / 2.0 * ((j * j - 1.0) / 2.0 - j.ln());
134
135        // TODO: shouldn't this be free of any phase dependency?
136        if j < 1.0 {
137            deviatoric_part * phase_coeff
138        } else {
139            (deviatoric_part + volumetric_part) * particle_phase
140        }
141    }
142
143    pub fn active_timestep_bounds(&self) -> ActiveTimestepBounds {
144        ActiveTimestepBounds::CONSTITUTIVE_MODEL_BOUND
145            | ActiveTimestepBounds::PARTICLE_VELOCITY_BOUND
146            | ActiveTimestepBounds::DEFORMATION_GRADIENT_CHANGE_BOUND
147    }
148
149    pub fn timestep_bound(
150        &self,
151        particle_density0: Real,
152        particle_velocity: &Vector<Real>,
153        particle_elastic_hardening: Real,
154        cell_width: Real,
155    ) -> Real {
156        let bulk_modulus = crate::utils::bulk_modulus_from_lame(self.lambda, self.mu);
157        let shear_modulus = crate::utils::shear_modulus_from_lame(self.lambda, self.mu);
158
159        let bound = ElasticitySoundSpeedTimestepBound {
160            alpha: self.cfl_coeff,
161            bulk_modulus: bulk_modulus * particle_elastic_hardening,
162            shear_modulus: shear_modulus * particle_elastic_hardening,
163        };
164        bound.timestep_bound(particle_density0, particle_velocity, cell_width)
165    }
166}
sparkl3d_core/dynamics/models/elasticity_neo_hookean.rs

sparkl3d_core/dynamics/models/
elasticity_neo_hookean.rs
