feat: Add elastic energy density function and refs

src/dynamics/models/constitutive_model.rs

@@ -4,6 +4,11 @@ use crate::math::{Matrix, Real};
 
 pub trait ConstitutiveModel: Send + Sync {
     fn is_fluid(&self) -> bool;
+
+    fn elastic_energy_density(&self, _deformation_gradient: Matrix<Real>) -> Real {
+        0.0
+    }
+
     fn pos_energy(&self, _particle: &Particle) -> Real {
         0.0
     }

src/dynamics/models/elasticity_corotated_linear.rs

@@ -17,6 +17,21 @@ impl ConstitutiveModel for CorotatedLinearElasticity {
         )
     }
 
+    // https://www.math.ucla.edu/~cffjiang/research/mpmcourse/mpmcourse.pdf#subsection.6.3
+    fn elastic_energy_density(&self, deformation_gradient: Matrix<Real>) -> Real {
+        let singular_values = deformation_gradient
+            .svd_unordered(false, false)
+            .singular_values;
+        let determinant: Real = singular_values.iter().product();
+
+        self.mu
+            * singular_values
+                .iter()
+                .map(|sigma| (sigma - 1.).powi(2))
+                .sum::<Real>()
+            + self.lambda / 2. * (determinant - 1.).powi(2)
+    }
+
     fn pos_energy(&self, particle: &Particle) -> Real {
         self.pos_energy(particle.deformation_gradient, particle.elastic_hardening)
     }

src/dynamics/models/elasticity_neo_hookean.rs

@@ -2,7 +2,7 @@ use crate::dynamics::models::{
     ActiveTimestepBounds, ConstitutiveModel, CoreConstitutiveModel, NeoHookeanElasticity,
 };
 use crate::dynamics::Particle;
-use crate::math::{Matrix, Real};
+use crate::math::{Matrix, Real, DIM};
 
 impl ConstitutiveModel for NeoHookeanElasticity {
     fn is_fluid(&self) -> bool {
@@ -17,6 +17,15 @@ impl ConstitutiveModel for NeoHookeanElasticity {
         )
     }
 
+    // https://www.math.ucla.edu/~cffjiang/research/mpmcourse/mpmcourse.pdf#subsection.6.2
+    fn elastic_energy_density(&self, deformation_gradient: Matrix<Real>) -> Real {
+        let determinant_log = deformation_gradient.determinant().ln();
+        self.mu / 2.
+            * ((deformation_gradient.transpose() * deformation_gradient).trace() - DIM as Real)
+            - self.mu * determinant_log
+            + self.lambda / 2. * determinant_log.powi(2)
+    }
+
     fn pos_energy(&self, particle: &Particle) -> Real {
         self.pos_energy(
             particle.phase,

src_core/dynamics/models/elasticity_corotated_linear.rs

@@ -76,7 +76,7 @@ impl CorotatedLinearElasticity {
 
     pub fn pos_energy(&self, deformation_gradient: Matrix<Real>, elastic_hardening: Real) -> Real {
         let j = deformation_gradient.determinant();
-        let mut pos_def = deformation_gradient.svd_unordered(true, true);
+        let mut pos_def = deformation_gradient.svd_unordered(true, true); // TODO: why compute U and V?
         pos_def.singular_values.apply(|e| *e = (*e - 1.0).max(0.0));
 
         let pos_dev_part = self.mu * elastic_hardening * pos_def.singular_values.norm_squared();

src_core/utils/physics.rs

@@ -14,6 +14,7 @@ pub fn inv_exact(e: Real) -> Real {
 }
 
 /// Computes the Lamé parameters (lambda, mu) from the young modulus and poisson ratio.
+/// https://encyclopediaofmath.org/wiki/Lam%C3%A9_constants
 pub fn lame_lambda_mu(young_modulus: Real, poisson_ratio: Real) -> (Real, Real) {
     (
         young_modulus * poisson_ratio / ((1.0 + poisson_ratio) * (1.0 - 2.0 * poisson_ratio)),