Try to let active stress depend on deviatoric strains

src/pulse/active_model.py

@@ -47,6 +47,8 @@ class for active models. Active models are used to incorporate
 import dolfinx
 import ufl
 
+from . import kinematics
+
 logger = logging.getLogger(__name__)
 
 
@@ -109,35 +111,48 @@ def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
             The active strain energy density function
         """
 
-    def S(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+    def S(self, C: ufl.core.expr.Expr, dev: bool = False) -> ufl.core.expr.Expr:
         """Cauchy stress tensor for the active model.
 
         Parameters
         ----------
-        F : ufl.core.expr.Expr
+        C : ufl.core.expr.Expr
             The right Cauchy-Green deformation tensor
+        dev : bool
+            Whether to compute the stress for the deviatoric part only
 
         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), Cdev)
 
-    def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+    def P(self, F: ufl.core.expr.Expr, dev: bool = False) -> ufl.core.expr.Expr:
         """First Piola-Kirchhoff stress tensor for the active model.
 
         Parameters
         ----------
         F : ufl.core.expr.Expr
             The deformation gradient
+        dev : bool
+            Whether to compute the stress for the deviatoric part only
 
         Returns
         -------
         ufl.core.expr.Expr
             The first Piola-Kirchhoff stress tensor
         """
-        return ufl.diff(self.strain_energy(F.T * F), F)
+        C = F.T * F
+        if dev:
+            Cdev = kinematics.Cdev(C)
+        else:
+            Cdev = C
+        return ufl.diff(self.strain_energy(Cdev), F)
 
 
 class Passive(ActiveModel):

src/pulse/active_stress.py

@@ -120,13 +120,15 @@ def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         else:
             raise NotImplementedError
 
-    def S(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+    def S(self, C: ufl.core.expr.Expr, dev: bool = False) -> ufl.core.expr.Expr:
         """Cauchy stress tensor for the active stress 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
 
         Returns
         -------

src/pulse/cardiac_model.py

@@ -20,9 +20,9 @@
 class ActiveModel(Protocol):
     def strain_energy(self, C: ufl.core.expr.Expr) -> 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: ...
 
-    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: ...
 
 
 class Compressibility(Protocol):
@@ -93,7 +93,7 @@ def strain_energy(
 
         psi = (
             self.material.strain_energy(Cdev)
-            + self.active.strain_energy(C)
+            + self.active.strain_energy(Cdev)
             + self.compressibility.strain_energy(C)
         )
         if C_dot is not None:
@@ -107,9 +107,8 @@ 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, dev=dev) + self.active.S(C) + self.compressibility.S(C)
+        S = self.material.S(C, dev=True) + self.active.S(C, dev=True) + self.compressibility.S(C)
         if C_dot is not None:
             S += self.viscoelasticity.S(C_dot)
         return S
@@ -120,9 +119,7 @@ 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()
-
-        P = self.material.P(F, dev=dev) + self.active.P(F) + self.compressibility.P(F)
+        P = self.material.P(F, dev=True) + self.active.P(F, dev=True) + self.compressibility.P(F)
         if F_dot is not None:
             P += self.viscoelasticity.P(F_dot)
         return P

tests/test_static_problem_Holzapfel.py

@@ -67,7 +67,11 @@ def dirichlet_bc(
 
     # With the HolzapfelOgden model the hydrostatic pressure
     # should equal the negative of the material parameter a
-    assert np.allclose(p.x.array, -material_params["a"].to_base_units())
+    # assert np.allclose(p.x.array, -material_params["a"].to_base_units())
+
+    # However when using the deviatoric strain only,
+    # the hydrostatic pressure should be zero
+    assert np.allclose(p.x.array, 0.0)
     # And with no external forces, there should be no displacement
     assert np.allclose(u.x.array, 0.0)