Add pre-stressing

Author: finsberg

.cspell_dict.txt

@@ -4,6 +4,7 @@ allreduce
 argsort
 astype
 atol
+Barnafi
 Bartolucci
 basix
 bcast
@@ -17,9 +18,11 @@ cardiomyocytes
 cellname
 Chiara
 cofac
+cofnorm
 cpus
 Cristobal
 Culloch
+Davide
 deviatoric
 dirichletbc
 dofmap
@@ -40,6 +43,7 @@ Francesco
 functionspace
 geodir
 gradu
+Grice
 Guccione
 Holohan
 Holzapfel
@@ -90,6 +94,7 @@ Quarteroni
 regazzoni
 Regazzoni
 Remedios
+Riccobelli
 rique
 rtol
 scifem
@@ -99,10 +104,12 @@ Silvano
 Sorine
 Stefano
 subplus
+Taras
 tomek
 Tomek
 TRAME
 ureg
+Usyk
 varepsilon
 Venant
 Verlag

docs/refs.bib

@@ -67,3 +67,13 @@ @article{land2017model
   year={2017},
   publisher={Elsevier}
 }
+
+@article{barnafi2024reconstructing,
+  title={Reconstructing relaxed configurations in elastic bodies: Mathematical formulations and numerical methods for cardiac modeling},
+  author={Barnafi, NA and Regazzoni, Francesco and Riccobelli, Davide},
+  journal={Computer Methods in Applied Mechanics and Engineering},
+  volume={423},
+  pages={116845},
+  year={2024},
+  publisher={Elsevier}
+}

src/pulse/__init__.py

@@ -22,6 +22,7 @@
     kinematics,
     material_model,
     material_models,
+    prestress,
     problem,
     units,
     utils,
@@ -30,15 +31,23 @@
 from .active_stress import ActiveStress
 from .boundary_conditions import BoundaryConditions, NeumannBC, RobinBC
 from .cardiac_model import CardiacModel
-from .compressibility import Compressibility, Compressible, Incompressible
+from .compressibility import (
+    Compressibility,
+    Compressible,
+    Compressible2,
+    Compressible3,
+    Incompressible,
+)
 from .geometry import Geometry, HeartGeometry, Marker
 from .material_model import HyperElasticMaterial, Material
 from .material_models import (
     Guccione,
     HolzapfelOgden,
     NeoHookean,
     SaintVenantKirchhoff,
+    Usyk,
 )
+from .prestress import PrestressProblem
 from .problem import BaseBC, DynamicProblem, StaticProblem
 from .units import Variable, ureg
 from .viscoelasticity import NoneViscoElasticity, ViscoElasticity, Viscous
@@ -87,4 +96,9 @@
     "ViscoElasticity",
     "Viscous",
     "__version__",
+    "prestress",
+    "PrestressProblem",
+    "Usyk",
+    "Compressible2",
+    "Compressible3",
 ]

src/pulse/active_model.py

@@ -121,6 +121,21 @@ def S(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         """
         return 2.0 * ufl.diff(self.strain_energy(F), F)
 
+    def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+        """First Piola-Kirchhoff stress tensor for the active model.
+
+        Parameters
+        ----------
+        F : ufl.core.expr.Expr
+            The deformation gradient
+
+        Returns
+        -------
+        ufl.core.expr.Expr
+            The first Piola-Kirchhoff stress tensor
+        """
+        return ufl.diff(self.strain_energy(F), F)
+
 
 class Passive(ActiveModel):
     """Active model with no active component.

src/pulse/boundary_conditions.py

@@ -39,6 +39,7 @@ class RobinBC:
     value: Variable
     marker: int
     damping: bool = False
+    perpendicular: bool = False
 
     def __post_init__(self):
         if not isinstance(self.value, Variable):

src/pulse/cardiac_model.py

@@ -18,12 +18,16 @@ 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 P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
+
 
 class Compressibility(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 P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
+
     def is_compressible(self) -> bool: ...
 
     def register(self, p: dolfinx.fem.Function | None) -> None: ...
@@ -32,12 +36,16 @@ 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 S(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
 
 
 class ViscoElasticity(Protocol):
     def strain_energy(self, C_dot: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
 
+    def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
+
     def S(self, C_dot: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
 
 
@@ -74,3 +82,14 @@ def S(
         if C_dot is not None:
             S += self.viscoelasticity.S(C_dot)
         return S
+
+    def P(
+        self,
+        F: ufl.core.expr.Expr,
+        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)
+        if F_dot is not None:
+            P += self.viscoelasticity.P(F_dot)
+        return P

src/pulse/compressibility.py

@@ -42,6 +42,11 @@ def S(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         """Cauchy stress tensor for the compressibility model."""
         return 2.0 * ufl.diff(self.strain_energy(C), C)
 
+    def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+        """First Piola-Kirchhoff stress tensor for the compressibility model."""
+        C = F.T * F
+        return ufl.diff(self.strain_energy(C), F)
+
     @abc.abstractmethod
     def is_compressible(self) -> bool:
         """Returns True if the material model is compressible."""
@@ -126,16 +131,38 @@ class Compressible2(Compressible):
     Strain energy density function is given by
 
     .. math::
-        \Psi = \kappa (J^2 - 1 - 2 \ln(J))
+        \Psi = \kappa / 4 (J^2 - 1 - 2 \ln(J))
 
     """
 
     kappa: Variable = field(default_factory=lambda: Variable(1e6, "Pa"))
 
     def __str__(self) -> str:
-        return "\u03ba (J ** 2 - 1 - 2 ln(J))"
+        return "\u03ba / 4 (J ** 2 - 1 - 2 ln(J))"
 
     def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         J = ufl.sqrt(ufl.det(C))
         kappa = self.kappa.to_base_units()
         return 0.25 * kappa * (J**2 - 1 - 2 * ufl.ln(J))
+
+
+@dataclass(slots=True)
+class Compressible3(Compressible):
+    r"""Compressible material model used in Usyk et al. 2002
+
+    Strain energy density function is given by
+
+    .. math::
+        \Psi = \kappa / 2 (J - 1) * \ln(J)
+
+    """
+
+    kappa: Variable = field(default_factory=lambda: Variable(5e4, "Pa"))
+
+    def __str__(self) -> str:
+        return "\u03ba / 2 * (J - 1) * ln(J)"
+
+    def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
+        J = ufl.sqrt(ufl.det(C))
+        kappa = self.kappa.to_base_units()
+        return 0.5 * kappa * (J - 1) * ufl.ln(J)

src/pulse/geometry.py

@@ -81,6 +81,26 @@ def _set_measures(self) -> None:
             metadata=self.metadata,
         )
 
+    def deform(self, u: dolfinx.fem.Function) -> None:
+        """Deform the geometry by a displacement field.
+        Note that this will modify the mesh geometry in place.
+
+        Parameters
+        ----------
+        u : dolfinx.fem.Function
+            Displacement field to deform the geometry with.
+        """
+        if not isinstance(u, dolfinx.fem.Function):
+            raise TypeError("Displacement field must be a dolfinx.fem.Function")
+
+        import scifem
+
+        self.mesh.geometry.x[:] += scifem.evaluate_function(
+            u,
+            self.mesh.geometry.x,
+            broadcast=False,
+        )
+
     @classmethod
     def from_cardiac_geometries(
         cls,

src/pulse/material_model.py

@@ -94,8 +94,8 @@ def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         .. math::
             \mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}}
         """
-        F = ufl.variable(F)
-        return ufl.diff(self.strain_energy(F), F)
+        C = F.T * F
+        return ufl.diff(self.strain_energy(C), F)
 
     def sigma(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         return kinematics.InversePiolaTransform(self.P(F), F)

src/pulse/material_models/__init__.py

@@ -2,5 +2,6 @@
 from .holzapfelogden import HolzapfelOgden
 from .neo_hookean import NeoHookean
 from .saint_venant_kirchhoff import SaintVenantKirchhoff
+from .usyk import Usyk
 
-__all__ = ["NeoHookean", "HolzapfelOgden", "SaintVenantKirchhoff", "Guccione"]
+__all__ = ["NeoHookean", "HolzapfelOgden", "SaintVenantKirchhoff", "Guccione", "Usyk"]