diff --git a/.cspell_dict.txt b/.cspell_dict.txt index 78f9a9c..795a5ff 100644 --- a/.cspell_dict.txt +++ b/.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 diff --git a/demo/geometries/biv_ellipsoid.py b/demo/geometries/biv_ellipsoid.py index 4579de2..74ad8ce 100644 --- a/demo/geometries/biv_ellipsoid.py +++ b/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))) diff --git a/src/pulse/cardiac_model.py b/src/pulse/cardiac_model.py index 212cf91..d641f06 100644 --- a/src/pulse/cardiac_model.py +++ b/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 diff --git a/src/pulse/kinematics.py b/src/pulse/kinematics.py index 13f323f..a939a78 100644 --- a/src/pulse/kinematics.py +++ b/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""" diff --git a/src/pulse/material_model.py b/src/pulse/material_model.py index f7be436..63e968a 100644 --- a/src/pulse/material_model.py +++ b/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) diff --git a/src/pulse/material_models/guccione.py b/src/pulse/material_models/guccione.py index 889993c..e146931 100644 --- a/src/pulse/material_models/guccione.py +++ b/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() diff --git a/src/pulse/material_models/holzapfelogden.py b/src/pulse/material_models/holzapfelogden.py index 15396b6..75d8d74 100644 --- a/src/pulse/material_models/holzapfelogden.py +++ b/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) diff --git a/src/pulse/material_models/neo_hookean.py b/src/pulse/material_models/neo_hookean.py index d536be0..5ed5880 100644 --- a/src/pulse/material_models/neo_hookean.py +++ b/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) diff --git a/src/pulse/material_models/usyk.py b/src/pulse/material_models/usyk.py index c6422da..efff7b0 100644 --- a/src/pulse/material_models/usyk.py +++ b/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() diff --git a/tests/test_cardiac_model.py b/tests/test_cardiac_model.py index decf4e0..d4ab96d 100644 --- a/tests/test_cardiac_model.py +++ b/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,7 +51,7 @@ 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) @@ -60,10 +59,14 @@ def test_CardiacModel_HolzapfelOgden(comp_model_cls, isotropy, mesh, u): # 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) diff --git a/tests/test_material_models.py b/tests/test_material_models.py index 14dd19d..62778d8 100644 --- a/tests/test_material_models.py +++ b/tests/test_material_models.py @@ -22,7 +22,7 @@ def test_holzapfel_ogden(params_func, expected_value, mesh, u) -> None: params = params_func() f0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0)) s0 = dolfinx.fem.Constant(mesh, (0.0, 1.0, 0.0)) - model = pulse.HolzapfelOgden(f0=f0, s0=s0, **params, deviatoric=False) + model = pulse.HolzapfelOgden(f0=f0, s0=s0, **params) u.interpolate(lambda x: x / 10) F = pulse.kinematics.DeformationGradient(u) @@ -57,7 +57,7 @@ def test_holzapfel_ogden_raises_MissingModelAttribute(params, attr): def test_holzapfel_ogden_neohookean(u): - model = pulse.HolzapfelOgden(a=1.0, deviatoric=False) + model = pulse.HolzapfelOgden(a=1.0) u.interpolate(lambda x: x / 10) F = pulse.kinematics.DeformationGradient(u) C = F.T * F @@ -72,7 +72,7 @@ def test_holzapfel_ogden_neohookean(u): def test_holzapfel_ogden_pure_fiber(u, mesh): f0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0)) - model = pulse.HolzapfelOgden(a_f=1.0, f0=f0, deviatoric=False) + model = pulse.HolzapfelOgden(a_f=1.0, f0=f0) u.interpolate(lambda x: x / 10) F = pulse.kinematics.DeformationGradient(u) C = F.T * F @@ -85,7 +85,7 @@ def test_holzapfel_ogden_pure_fiber(u, mesh): def test_holzapfel_ogden_pure_fiber_sheets(u, mesh): f0 = dolfinx.fem.Constant(mesh, (1.0, 0.0, 0.0)) - model = pulse.HolzapfelOgden(a_fs=1.0, f0=f0, s0=f0, deviatoric=False) + model = pulse.HolzapfelOgden(a_fs=1.0, f0=f0, s0=f0) u.interpolate(lambda x: x / 10) F = pulse.kinematics.DeformationGradient(u) C = F.T * F @@ -97,7 +97,7 @@ def test_holzapfel_ogden_pure_fiber_sheets(u, mesh): def test_neo_hookean(u, mesh): - model = pulse.NeoHookean(mu=1.0, deviatoric=False) + model = pulse.NeoHookean(mu=1.0) u.interpolate(lambda x: x / 10) F = pulse.kinematics.DeformationGradient(u) C = F.T * F @@ -128,7 +128,7 @@ def test_saint_venant_kirchhoff(u): def test_guccione_isotropic(u): C = 10.0 bf = bt = bfs = 1.0 - model = pulse.Guccione(C=C, bf=bf, bt=bt, bfs=bfs, deviatoric=False) + model = pulse.Guccione(C=C, bf=bf, bt=bt, bfs=bfs) assert model.is_isotropic() u.interpolate(lambda x: x / 10) @@ -151,7 +151,7 @@ def test_guccione_anisotropic(u, mesh): bf = 1.0 bt = 2.0 bfs = 3.0 - model = pulse.Guccione(C=C, bf=bf, bt=bt, bfs=bfs, f0=f0, s0=s0, n0=n0, deviatoric=False) + model = pulse.Guccione(C=C, bf=bf, bt=bt, bfs=bfs, f0=f0, s0=s0, n0=n0) assert not model.is_isotropic() u.interpolate(lambda x: x / 10) diff --git a/tests/test_static_problem_Holzapfel.py b/tests/test_static_problem_Holzapfel.py index a8e085d..6dd6584 100644 --- a/tests/test_static_problem_Holzapfel.py +++ b/tests/test_static_problem_Holzapfel.py @@ -112,10 +112,7 @@ def test_CompressibleProblem_and_boundary_conditions(mesh): ) f0 = dolfinx.fem.Constant(mesh, PETSc.ScalarType((1.0, 0.0, 0.0))) - material = pulse.NeoHookean( - mu=dolfinx.fem.Constant(mesh, PETSc.ScalarType(15.0)), - deviatoric=True, - ) + material = pulse.NeoHookean(mu=dolfinx.fem.Constant(mesh, PETSc.ScalarType(15.0))) Ta = dolfinx.fem.Constant(mesh, PETSc.ScalarType(0.0)) active_model = pulse.ActiveStress(f0, activation=Ta)