Remove deviatoric arguments from material models

Author: finsberg

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):
@@ -66,42 +67,6 @@ def __post_init__(self):
         logger.debug(f"  Compressibility: {type(self.compressibility).__name__}")
         logger.debug(f"  Viscoelasticity: {type(self.viscoelasticity).__name__}")
 
-    def Cdev(self, 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
-
-    def Fdev(self, 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 strain_energy(
         self,
         C: ufl.core.expr.Expr,
@@ -122,7 +87,7 @@ def strain_energy(
             The total strain energy density
         """
         if self.compressibility.is_compressible():
-            Cdev = self.Cdev(C)
+            Cdev = kinematics.Cdev(C)
         else:
             Cdev = C
 
@@ -142,12 +107,9 @@ def S(
         C_dot: ufl.core.expr.Expr | None = None,
     ) -> ufl.core.expr.Expr:
         """Cauchy stress for the cardiac model."""
-        if self.compressibility.is_compressible():
-            Cdev = self.Cdev(C)
-        else:
-            Cdev = C
+        dev = self.compressibility.is_compressible()
 
-        S = self.material.S(Cdev) + 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
@@ -158,13 +120,9 @@ def P(
         F_dot: ufl.core.expr.Expr | None = None,
     ) -> ufl.core.expr.Expr:
         """First Piola-Kirchhoff stress for the cardiac model."""
+        dev = self.compressibility.is_compressible()
 
-        if self.compressibility.is_compressible():
-            Fdev = self.Fdev(F)
-        else:
-            Fdev = F
-
-        P = self.material.P(Fdev) + self.active.P(F) + self.compressibility.P(F)
+        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)
@@ -66,7 +65,7 @@ def test_CardiacModel_HolzapfelOgden(comp_model_cls, isotropy, mesh, u):
 
         # 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
+        # + 1e6 * (1.1 ** 3 * math.log(1.1 ** 3) - 1.1 ** 3 + 1) = 49573.547958669194
         assert math.isclose(value, 49573.547958669194)
 
 
@@ -78,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()
@@ -99,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,))
@@ -112,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(
@@ -130,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)