Make sure we do deviatoric and volumetric split for compressible models

.cspell_dict.txt

@@ -6,6 +6,7 @@ atol
 basix
 cellname
 cofac
+deviatoric
 dirichletbc
 dofs
 dolfinx
@@ -22,6 +23,8 @@ hyperelastic
 isclose
 isochoric
 isscalar
+libgl
+libxrender
 lmbda
 mechanicsproblem
 meshtags
@@ -34,13 +37,11 @@ numpy
 Paraview
 petsc
 Piola
+PYVISTA
 rtol
 subplus
+TRAME
 varepsilon
 Venant
 XDMF
-PYVISTA
-TRAME
-libgl
-libxrender
 xvfb

_toc.yml

@@ -5,6 +5,7 @@ parts:
   - caption: Demos
     chapters:
     - file: "demo/unit_cube"
+    - file: "demo/lv_ellipsoid"
   - caption: Community
     chapters:
       - file: "CONTRIBUTING"

src/fenicsx_pulse/cardiac_model.py

@@ -13,6 +13,7 @@ def Fe(self, F) -> dolfinx.fem.Form: ...
 
 class Compressibility(Protocol):
     def strain_energy(self, J) -> dolfinx.fem.Form: ...
+    def is_compressible(self) -> bool: ...
     def register(self, p: dolfinx.fem.Function | None) -> None: ...
 
 
@@ -32,7 +33,14 @@ def strain_energy(self, F, p: dolfinx.fem.Function | None = None):
         # part of the deformation gradient
         Fe = self.active.Fe(F)
         J = kinematics.Jacobian(Fe)
-        Jm13 = J ** (-1 / 3)
+
+        # If model is compressible we need to to a
+        # deviatoric / volumetric split
+        if self.compressibility.is_compressible():
+            Jm13 = J ** (-1 / 3)
+        else:
+            Jm13 = 1.0
+
         return (
             self.material.strain_energy(Jm13 * Fe)
             + self.active.strain_energy(Jm13 * F)

src/fenicsx_pulse/compressibility.py

@@ -15,6 +15,10 @@ def strain_energy(self, J: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
     def register(self, *args, **kwargs) -> None:
         pass
 
+    @abc.abstractmethod
+    def is_compressible(self) -> bool:
+        pass
+
 
 @dataclass(slots=True)
 class Incompressible(Compressibility):
@@ -28,10 +32,16 @@ def strain_energy(self, J: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
             raise exceptions.MissingModelAttribute(attr="p", model=type(self).__name__)
         return self.p * (J - 1.0)
 
+    def is_compressible(self) -> bool:
+        return False
+
 
 @dataclass(slots=True)
 class Compressible(Compressibility):
     kappa: float | dolfinx.fem.Function | dolfinx.fem.Constant = 1e3
 
     def strain_energy(self, J: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         return self.kappa * (J * ufl.ln(J) - J + 1)
+
+    def is_compressible(self) -> bool:
+        return True

tests/test_cardiac_model.py

@@ -23,4 +23,4 @@ def test_CardiacModel(mesh, u):
     F = fenicsx_pulse.kinematics.DeformationGradient(u)
     psi = model.strain_energy(F)
     value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(psi * ufl.dx))
-    assert math.isclose(value, 103.22036041941614)
+    assert math.isclose(value, 49.57354795865461)

tests/test_mechanicsproblem.py

@@ -5,7 +5,7 @@
 import numpy as np
 
 
-def test_MechanicsProblem_and_boundary_conditions(mesh):
+def test_MechanicsProblemMixed_and_boundary_conditions(mesh):
     boundaries = [
         (1, 2, lambda x: np.isclose(x[0], 0)),
         (2, 2, lambda x: np.isclose(x[0], 1)),
@@ -59,7 +59,7 @@ def dirichlet_bc(
         body_force=(body_force,),
     )
 
-    problem = fenicsx_pulse.MechanicsProblem(model=model, geometry=geo, bcs=bcs)
+    problem = fenicsx_pulse.MechanicsProblemMixed(model=model, geometry=geo, bcs=bcs)
     problem.solve()
 
     u = problem.state.sub(0).collapse()