Correct stress rate implementation

python/comfe/cdm.py

@@ -9,7 +9,7 @@
 from petsc4py import PETSc
 from pydantic import BaseModel
 
-from .comfe import jaumann_rotation  # , jaumann_rotation_expensive
+from .comfe import jaumann_rotation, jaumann_rotation_expensive
 from .helpers import LogMixin, QuadratureEvaluator, QuadratureRule, diagonal_mass, set_mesh_coordinates
 from .laws import ConstitutiveModel, RustConstitutiveModel
 
@@ -153,7 +153,7 @@ def __init__(
         f_int_form = df.fem.form(f_int_ufl)
         try:
             L_evaluator = QuadratureEvaluator(
-                self._as_3d_tensor(ufl.nabla_grad(v)),
+                self._as_3d_tensor(ufl.grad(v)),
                 function_space.mesh,
                 quadrature_rule,
             )

python/comfe/laws.py

@@ -4,7 +4,7 @@
 import numpy as np
 from pydantic import BaseModel
 
-from .comfe import PyGradientJH23D, PyJH23D, PyLinearElastic3D, PyLinElas3D
+from .comfe import PyGradientJH23D, PyHypoelasticity3D, PyJH23D, PyLinearElastic3D, PyLinElas3D
 from .helpers import QuadratureRule
 
 __all__ = [
@@ -15,8 +15,9 @@
     "PyJH23D",
     "PyGradientJH23D",
     "PyLinearElastic3D",
+    "PyHypoelasticity3D",
 ]
-RustConstitutiveModel = PyLinElas3D | PyJH23D | PyLinearElastic3D | PyGradientJH23D
+RustConstitutiveModel = PyLinElas3D | PyJH23D | PyLinearElastic3D | PyGradientJH23D | PyHypoelasticity3D
 
 
 class QuadratureModel(ABC):

src/hypoelasticity.rs

@@ -0,0 +1,71 @@
+use crate::interfaces::{ConstitutiveModel, QDim, QValueInput, QValueOutput, Q};
+use crate::stress_strain::{
+    mandel_decomposition, mandel_rate_from_velocity_gradient, volumetric, MANDEL_IDENTITY,
+};
+
+use std::collections::HashMap;
+#[derive(Debug)]
+pub struct Hypoelasticity3D {
+    mu: f64,
+    lambda: f64,
+}
+
+impl ConstitutiveModel for Hypoelasticity3D {
+    fn new(parameters: &HashMap<String, f64>) -> Option<Self> {
+        Some(Self {
+            mu: *parameters.get("mu").unwrap(),
+            lambda: *parameters.get("lambda").unwrap(),
+        })
+    }
+    fn evaluate_ip(&self, ip: usize, del_t: f64, input: &QValueInput, output: &mut QValueOutput) {
+        let velocity_gradient = input
+            .get_tensor::<{ Q::VelocityGradient.dim() }, { Q::VelocityGradient.size() }>(
+                Q::VelocityGradient,
+                ip,
+            );
+        let d_eps = mandel_rate_from_velocity_gradient(&velocity_gradient);
+        let trace_d_eps = volumetric(&d_eps) * 3.0;
+        let sigma_0 = input.get_vector::<{ Q::MandelStress.size() }>(Q::MandelStress, ip);
+
+        let sigma_1 = sigma_0
+            + 2. * self.mu * d_eps * del_t
+            + self.lambda * trace_d_eps * del_t * MANDEL_IDENTITY;
+
+        output.set_vector(Q::MandelStress, ip, sigma_1);
+    }
+
+    /// Returns the physical quantities that are required as input for the
+    /// constitutive model together with their dimensions.
+    fn define_input(&self) -> HashMap<Q, QDim> {
+        HashMap::from([
+            (Q::VelocityGradient, QDim::SquareTensor(3)),
+        ])
+    }
+
+    /// Returns the physical quantities that are needed as internal variables
+    /// for the constitutive model together with their dimensions. These Variables are
+    /// stored both in in the input and the output.
+    fn define_history(&self) -> HashMap<Q, QDim> {
+        HashMap::from([
+            (Q::MandelStress, QDim::Vector(6)),
+        ])
+    }
+
+    /// Returns the physical quantities that are needed as output, but are not
+    /// necessarily needed in oredr to calculate the constitutive model. An example is
+    /// the consistent tangent which is not needed for the calculation of the stresses
+    /// and is therefore purely an output quantity.
+    fn define_output(&self) -> HashMap<Q, QDim> {
+        HashMap::from([(Q::MandelStress, QDim::Vector(6))])
+    }
+
+    /// Returns the physical quantities that are optional output of the constitutive
+    /// model. These quantities are not needed for the calculation of the stresses
+    /// but can be useful for postprocessing.
+    fn define_optional_output(&self) -> HashMap<Q, QDim> {
+        HashMap::from([])
+    }
+    fn define_optional_history(&self) -> HashMap<Q, QDim> {
+        HashMap::from([])
+    }
+}

src/lib.rs

@@ -7,6 +7,7 @@ use crate::jh2::JH23D;
 use crate::generic_jh2::GenericJH23D;
 use crate::gradient_jh2::GradientJH23D;
 use crate::smallstrain::linear_elastic::LinearElastic3D;
+use crate::hypoelasticity::Hypoelasticity3D;
 //use crate::stress_strain;
 use nalgebra::{Const, DVectorView, DVectorViewMut, Dyn, SMatrix};
 use numpy::{PyReadonlyArray1, PyReadwriteArray1};
@@ -19,6 +20,7 @@ pub mod jh2;
 //pub mod jh_concrete;
 pub mod generic_jh2;
 pub mod gradient_jh2;
+pub mod hypoelasticity;
 pub mod smallstrain;
 pub mod stress_strain;
 
@@ -398,6 +400,7 @@ fn comfe(_py: Python, m: &PyModule) -> PyResult<()> {
     impl_constitutive_model!(PyJH23D, JH23D, m);
     impl_constitutive_model!(PyGradientJH23D, GradientJH23D, m);
     impl_constitutive_model!(PyLinElas3D, LinearElastic3D, m);
+    impl_constitutive_model!(PyHypoelasticity3D, Hypoelasticity3D, m);
     m.add_function(wrap_pyfunction!(py_jaumann_rotation, m)?)?;
     m.add_function(wrap_pyfunction!(py_jaumann_rotation_expensive, m)?)?;
     Ok(())