Make sure we use deviatoric strains when using compressible models

.cspell_dict.txt

@@ -19,6 +19,7 @@ Bucelli
 Bursi
 caplog
 cardiomyocytes
+Cdev
 cellname
 Chiara
 CMAME
@@ -41,6 +42,7 @@ Elife
 Elsevier
 endo
 Erlicher
+Fdev
 Fedele
 fenics
 fenicsx

demo/geometries/biv_ellipsoid.py

@@ -157,7 +157,6 @@ def update_marker(old_name, new_id):
 # \Psi_{NH} = \frac{\mu}{2} (I_1 - 3)
 # $$
 #
-# `deviatoric=True` (default) means we use the isochoric invariant $\bar{I}_1 = J^{-2/3} I_1$.
 
 material = pulse.NeoHookean(mu=dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(15.0)))
 

src/pulse/cardiac_model.py

@@ -11,6 +11,7 @@
 import dolfinx
 import ufl
 
+from . import kinematics
 from .viscoelasticity import NoneViscoElasticity
 
 logger = logging.getLogger(__name__)
@@ -39,9 +40,9 @@ def register(self, p: dolfinx.fem.Function | None) -> None: ...
 class HyperElasticMaterial(Protocol):
     def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
 
-    def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
+    def P(self, F: ufl.core.expr.Expr, dev: bool) -> ufl.core.expr.Expr: ...
 
-    def S(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
+    def S(self, C: ufl.core.expr.Expr, dev: bool) -> ufl.core.expr.Expr: ...
 
 
 class ViscoElasticity(Protocol):
@@ -71,8 +72,27 @@ def strain_energy(
         C: ufl.core.expr.Expr,
         C_dot: ufl.core.expr.Expr | None = None,
     ) -> ufl.core.expr.Expr:
+        """Total strain energy for the cardiac model.
+
+        Parameters
+        ----------
+        C : ufl.core.expr.Expr
+            Right Cauchy-Green deformation tensor
+        C_dot : ufl.core.expr.Expr | None, optional
+            Time derivative of the right Cauchy-Green deformation tensor, by default None
+
+        Returns
+        -------
+        ufl.core.expr.Expr
+            The total strain energy density
+        """
+        if self.compressibility.is_compressible():
+            Cdev = kinematics.Cdev(C)
+        else:
+            Cdev = C
+
         psi = (
-            self.material.strain_energy(C)
+            self.material.strain_energy(Cdev)
             + self.active.strain_energy(C)
             + self.compressibility.strain_energy(C)
         )
@@ -87,8 +107,9 @@ def S(
         C_dot: ufl.core.expr.Expr | None = None,
     ) -> ufl.core.expr.Expr:
         """Cauchy stress for the cardiac model."""
+        dev = self.compressibility.is_compressible()
 
-        S = self.material.S(C) + self.active.S(C) + self.compressibility.S(C)
+        S = self.material.S(C, dev=dev) + self.active.S(C) + self.compressibility.S(C)
         if C_dot is not None:
             S += self.viscoelasticity.S(C_dot)
         return S
@@ -99,7 +120,9 @@ def P(
         F_dot: ufl.core.expr.Expr | None = None,
     ) -> ufl.core.expr.Expr:
         """First Piola-Kirchhoff stress for the cardiac model."""
-        P = self.material.P(F) + self.active.P(F) + self.compressibility.P(F)
+        dev = self.compressibility.is_compressible()
+
+        P = self.material.P(F, dev=dev) + self.active.P(F) + self.compressibility.P(F)
         if F_dot is not None:
             P += self.viscoelasticity.P(F_dot)
         return P

src/pulse/kinematics.py

@@ -2,6 +2,44 @@
 import ufl
 
 
+def Fdev(F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+    """Deviatoric part of the deformation gradient.
+
+    Parameters
+    ----------
+    F : ufl.core.expr.Expr
+        Deformation gradient
+
+    Returns
+    -------
+    ufl.core.expr.Expr
+        Deviatoric part of the deformation gradient
+    """
+    J = ufl.det(F)
+    dim = F.ufl_shape[0]
+    Fdev = J ** (-1.0 / dim) * F
+    return Fdev
+
+
+def Cdev(C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+    """Deviatoric part of the right Cauchy-Green deformation tensor.
+
+    Parameters
+    ----------
+    C : ufl.core.expr.Expr
+        Right Cauchy-Green deformation tensor
+
+    Returns
+    -------
+    ufl.core.expr.Expr
+        Deviatoric part of the right Cauchy-Green deformation tensor
+    """
+    J = ufl.sqrt(ufl.det(C))
+    dim = C.ufl_shape[0]
+    Cdev = J ** (-2.0 / dim) * C
+    return Cdev
+
+
 # Second order identity tensor
 def SecondOrderIdentity(F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
     """Return identity with same dimension as input"""

src/pulse/material_model.py

@@ -72,13 +72,16 @@ def strain_energy(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
             The strain energy density
         """
 
-    def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+    def P(self, F: ufl.core.expr.Expr, dev: bool = True) -> ufl.core.expr.Expr:
         r"""First Piola-Kirchhoff stress tensor
 
         Parameters
         ----------
         F : ufl.core.expr.Expr
            The deformation gradient
+        dev : bool
+            Whether to compute the stress for the deviatoric part only
+            This should be True for compressible materials
 
         Returns
         -------
@@ -95,22 +98,34 @@ def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
             \mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}}
         """
         C = F.T * F
-        return ufl.diff(self.strain_energy(C), F)
+        if dev:
+            Cdev = kinematics.Cdev(C)
+        else:
+            Cdev = C
+        return ufl.diff(self.strain_energy(Cdev), F)
 
     def sigma(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         return kinematics.InversePiolaTransform(self.P(F), F)
 
-    def S(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+    def S(self, C: ufl.core.expr.Expr, dev: bool = True) -> ufl.core.expr.Expr:
         """Cauchy stress tensor for the hyperelastic material model.
 
         Parameters
         ----------
         C : ufl.core.expr.Expr
             The right Cauchy-Green deformation tensor
+        dev : bool
+            Whether to compute the stress for the deviatoric part only
+            This should be True for compressible materials
 
         Returns
         -------
         ufl.core.expr.Expr
             The Cauchy stress tensor
+
         """
-        return 2.0 * ufl.diff(self.strain_energy(C), C)
+        if dev:
+            Cdev = kinematics.Cdev(C)
+        else:
+            Cdev = C
+        return 2.0 * ufl.diff(self.strain_energy(Cdev), C)

src/pulse/material_models/guccione.py

@@ -71,7 +71,6 @@ class Guccione(HyperElasticMaterial):
     bf: Variable = field(default_factory=lambda: Variable(8.0, "dimensionless"))
     bt: Variable = field(default_factory=lambda: Variable(2.0, "dimensionless"))
     bfs: Variable = field(default_factory=lambda: Variable(4.0, "dimensionless"))
-    deviatoric: bool = True
 
     def __post_init__(self):
         # Check that all values are positive
@@ -130,9 +129,6 @@ def is_isotropic(self) -> bool:
 
     def _Q(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         dim = C.ufl_shape[0]
-        if self.deviatoric:
-            Jm23 = pow(ufl.det(C), -1.0 / dim)
-            C *= Jm23  # Make C deviatoric
         E = 0.5 * (C - ufl.Identity(dim))
 
         bt = self.bt.to_base_units()

src/pulse/material_models/holzapfelogden.py

@@ -118,7 +118,6 @@ class HolzapfelOgden(HyperElasticMaterial):
     b_fs: Variable = field(default_factory=lambda: Variable(0.0, "dimensionless"))
     use_subplus: bool = field(default=True, repr=False)
     use_heaviside: bool = field(default=True, repr=False)
-    deviatoric: bool = field(default=True, repr=False)
 
     _W1_func: Invariant = field(
         init=False,
@@ -331,15 +330,20 @@ def _W8fs(self, I8fs: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         return self._W8fs_func(I8fs)
 
     def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
-        dim = C.ufl_shape[0]
+        r"""Strain energy density function
 
-        I1 = invariants.I1(C)
+        Parameters
+        ----------
+        C : ufl.core.expr.Expr
+            Right Cauchy-Green deformation tensor
+
+        Returns
+        -------
+        ufl.core.expr.Expr
+            The strain energy density
+        """
 
-        if self.deviatoric:
-            Jm23 = pow(ufl.det(C), -1.0 / dim)
-            C *= Jm23  # Make C deviatoric
-            Jm23 = pow(invariants.I3(C), -1 / dim)
-            I1 *= Jm23  # Convert to deviatoric I1 - see https://arxiv.org/pdf/2009.08754
+        I1 = invariants.I1(C)
 
         I4f = self._I4f(C)
         I4s = self._I4s(C)

src/pulse/material_models/neo_hookean.py

@@ -34,7 +34,6 @@ class NeoHookean(HyperElasticMaterial):
     """
 
     mu: Variable = field(default_factory=lambda: Variable(15.0, "kPa"))
-    deviatoric: bool = field(default=True)
 
     def __post_init__(self):
         if not isinstance(self.mu, Variable):
@@ -61,10 +60,6 @@ def parameters(self) -> dict[str, Variable]:
     def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         I1 = invariants.I1(C)
         dim = C.ufl_shape[0]
-        if self.deviatoric:
-            # Deviatoric strain energy
-            Jm23 = pow(ufl.det(C), -1.0 / dim)
-            I1 *= Jm23
         mu = self.mu.to_base_units()
         return 0.5 * mu * (I1 - dim)
 

src/pulse/material_models/usyk.py

@@ -66,7 +66,6 @@ class Usyk(HyperElasticMaterial):
     bfs: Variable = field(default_factory=lambda: Variable(12.0, "dimensionless"))
     bfn: Variable = field(default_factory=lambda: Variable(3.0, "dimensionless"))
     bsn: Variable = field(default_factory=lambda: Variable(3.0, "dimensionless"))
-    deviatoric: bool = True
 
     def __post_init__(self):
         # Check that all values are positive
@@ -141,9 +140,6 @@ def is_isotropic(self) -> bool:
 
     def _Q(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         dim = C.ufl_shape[0]
-        if self.deviatoric:
-            Jm23 = pow(ufl.det(C), -1.0 / dim)
-            C *= Jm23  # Make C deviatoric
         E = 0.5 * (C - ufl.Identity(dim))
 
         bf = self.bf.to_base_units()

tests/test_cardiac_model.py

@@ -27,7 +27,6 @@ def test_CardiacModel_HolzapfelOgden(comp_model_cls, isotropy, mesh, u):
         b_f=0.0,
         a_fs=0.0,
         b_fs=0.0,
-        deviatoric=False,
     )
     comp_model = comp_model_cls()
     active_model = pulse.ActiveStress(f0, isotropy=isotropy)
@@ -52,18 +51,22 @@ def test_CardiacModel_HolzapfelOgden(comp_model_cls, isotropy, mesh, u):
     # )
 
     # F = I + 0.1 I, C = 1.21 I, I4f = 1.21
-    # J = det(F) = 1.1 ** 3
+    # J = det(F) = 1.1 ** 3, J^{-2/3} = 1.1 ** -2
 
     if isinstance(comp_model, pulse.compressibility.Incompressible):
         # psi = 0.5 * a * (I1 - 3) + 0.5 * a_f * (I4f - 1)**2 + p (J - 1)
         # psi = 0.5 * 1000*1 * (3 * 1.21 - 3) + 0.5 * 1000.0 * 1.0 * (1.21 - 1)**2 +
         # 1000.0 * (1.1 ** 3 - 1) = 668.05
         assert math.isclose(value, 668.05)
     else:
-        # psi = 0.5 * a * (I1 - 3) + 0.5 * a_f * (I4f - 1)**2 + kappa * (J * ln(J) - J + 1)
-        # psi = 0.5 * 1000*1 * (3 * 1.21 - 3) + 0.5 * 1000.0 * 1.0 * (1.21 - 1)**2 +
-        # 1e6 * (1.1 ** 3 * math.log(1.1 ** 3) - 1.1 ** 3 + 1) = 49910.5979586692
-        assert math.isclose(value, 49910.5979586692)
+        # J^{-2/3} = 1.1 ** -2
+        # psi = 0.5 * a * (J^{-2/3} * I1 - 3) + 0.5 * a_f * (^{-2/3} * I4f - 1)**2
+        # + kappa * (J * ln(J) - J + 1)
+
+        # psi = 0.5 * 1000*1 * (1.1 ** -2  * 3 * 1.21 - 3)
+        # + 0.5 * 1000.0 * 1.0 * (1.1 ** -2  * 1.21 - 1)**2
+        # + 1e6 * (1.1 ** 3 * math.log(1.1 ** 3) - 1.1 ** 3 + 1) = 49573.547958669194
+        assert math.isclose(value, 49573.547958669194)
 
 
 @pytest.mark.parametrize("isotropy", (pulse.active_stress.ActiveStressModels.transversely,))
@@ -74,7 +77,6 @@ def test_CardiacModel_HolzapfelOgden(comp_model_cls, isotropy, mesh, u):
 def test_CardiacModel_NeoHookean(comp_model_cls, isotropy, mesh, u):
     material = pulse.NeoHookean(
         mu=dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(15.0)),
-        deviatoric=False,
     )
     f0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0))
     comp_model = comp_model_cls()
@@ -95,7 +97,7 @@ def test_CardiacModel_NeoHookean(comp_model_cls, isotropy, mesh, u):
     if isinstance(comp_model, pulse.compressibility.Incompressible):
         assert math.isclose(value, 4725.331000000082)
     else:
-        assert math.isclose(value, 54298.54795867078)
+        assert math.isclose(value, 49573.54795867355)
 
 
 @pytest.mark.parametrize("isotropy", (pulse.active_stress.ActiveStressModels.transversely,))
@@ -108,7 +110,7 @@ def test_CardiacModel_Guccione(comp_model_cls, isotropy, mesh, u):
     s0 = dolfinx.fem.Constant(mesh, (0.0, 1.0, 0.0))
     n0 = dolfinx.fem.Constant(mesh, (0.0, 0.0, 1.0))
     material_params = pulse.Guccione.default_parameters()
-    material = pulse.Guccione(f0=f0, s0=s0, n0=n0, **material_params, deviatoric=False)
+    material = pulse.Guccione(f0=f0, s0=s0, n0=n0, **material_params)
     active_model = pulse.ActiveStress(f0, isotropy=isotropy)
     comp_model = comp_model_cls()
     model = pulse.CardiacModel(
@@ -126,4 +128,4 @@ def test_CardiacModel_Guccione(comp_model_cls, isotropy, mesh, u):
     if isinstance(comp_model, pulse.compressibility.Incompressible):
         assert math.isclose(value, 141.78170311802802)
     else:
-        assert math.isclose(value, 49714.998661786354)
+        assert math.isclose(value, 49573.54795866912)