Change the waywe pass parameters to material model

_config.yml

@@ -4,7 +4,7 @@
 title: pulse2
 author: Henrik Finsberg
 logo: "docs/logo.png"
-copyright: "2022"
+copyright: "2023"
 only_build_toc_files: true
 
 # Force re-execution of notebooks on each build.

_toc.yml

@@ -4,7 +4,7 @@ root: README
 parts:
   - caption: Demos
     chapters:
-    - file: "demos/simple_pv_loop.py"
+    - file: "demos/simple_pvloop.py"
   - caption: Community
     chapters:
       - file: "CONTRIBUTING"

demos/simple_pvloop.py

@@ -37,7 +37,7 @@
 geo.mesh.coordinates()[:] /= 3
 
 material_params = pulse2.HolzapfelOgden.transversely_isotropic_parameters()
-material = pulse2.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)
+material = pulse2.HolzapfelOgden(f0=geo.f0, s0=geo.s0, parameters=material_params)
 
 Ta = dolfin.Constant(0.0)
 active_model = pulse2.ActiveStress(geo.f0, activation=Ta)

src/pulse2/holzapfelogden.py

@@ -1,3 +1,4 @@
+from __future__ import annotations
 import typing
 from dataclasses import dataclass
 from dataclasses import field
@@ -14,6 +15,21 @@
 from .material_model import HyperElasticMaterial
 
 
+class HolzapfelOgdenParameters(typing.TypedDict):
+    a: float | dolfin.Function | dolfin.Constant
+    b: float | dolfin.Function | dolfin.Constant
+    a_f: float | dolfin.Function | dolfin.Constant
+    b_f: float | dolfin.Function | dolfin.Constant
+    a_s: float | dolfin.Function | dolfin.Constant
+    b_s: float | dolfin.Function | dolfin.Constant
+    a_fs: float | dolfin.Function | dolfin.Constant
+    b_fs: float | dolfin.Function | dolfin.Constant
+
+
+def default_parameters() -> HolzapfelOgdenParameters:
+    return HolzapfelOgden.transversely_isotropic_parameters()
+
+
 @dataclass(slots=True)
 class HolzapfelOgden(HyperElasticMaterial):
     r"""
@@ -25,28 +41,31 @@ class HolzapfelOgden(HyperElasticMaterial):
         Function representing the direction of the fibers
     s0: dolfin.Function | dolfin.Constant | None
         Function representing the direction of the sheets
-    a: float | dolfin.Function | dolfin.Constant
-        Material parameter, by default 0.0
-    b: float | dolfin.Function | dolfin.Constant
-        Material parameter, by default 0.0
-    a_f: float | dolfin.Function | dolfin.Constant
-        Material parameter, by default 0.0
-    b_f: float | dolfin.Function | dolfin.Constant
-        Material parameter, by default 0.0
-    a_s: float | dolfin.Function | dolfin.Constant
-        Material parameter, by default 0.0
-    b_s: float | dolfin.Function | dolfin.Constant
-        Material parameter, by default 0.0
-    a_fs: float | dolfin.Function | dolfin.Constant
-        Material parameter, by default 0.0
-    b_fs: float | dolfin.Function | dolfin.Constant
-        Material parameter, by default 0.0
     use_subplus: bool
         Use subplus function when computing anisotropic contribution,
         by default True
     use_heaviside: bool
         Use heaviside function when computing anisotropic contribution,
         by default True
+    parameters : dict[str, float | dolfin.Function | dolfin.Constant]
+        Parameters in the model. Possible values are
+
+            a: float | dolfin.Function | dolfin.Constant
+                Material parameter, by default 0.0
+            b: float | dolfin.Function | dolfin.Constant
+                Material parameter, by default 0.0
+            a_f: float | dolfin.Function | dolfin.Constant
+                Material parameter, by default 0.0
+            b_f: float | dolfin.Function | dolfin.Constant
+                Material parameter, by default 0.0
+            a_s: float | dolfin.Function | dolfin.Constant
+                Material parameter, by default 0.0
+            b_s: float | dolfin.Function | dolfin.Constant
+                Material parameter, by default 0.0
+            a_fs: float | dolfin.Function | dolfin.Constant
+                Material parameter, by default 0.0
+            b_fs: float | dolfin.Function | dolfin.Constant
+                Material parameter, by default 0.0
 
     Notes
     -----
@@ -88,14 +107,7 @@ class HolzapfelOgden(HyperElasticMaterial):
     """
     f0: dolfin.Function | dolfin.Constant | None = None
     s0: dolfin.Function | dolfin.Constant | None = None
-    a: float | dolfin.Function | dolfin.Constant = 0.0
-    b: float | dolfin.Function | dolfin.Constant = 0.0
-    a_f: float | dolfin.Function | dolfin.Constant = 0.0
-    b_f: float | dolfin.Function | dolfin.Constant = 0.0
-    a_s: float | dolfin.Function | dolfin.Constant = 0.0
-    b_s: float | dolfin.Function | dolfin.Constant = 0.0
-    a_fs: float | dolfin.Function | dolfin.Constant = 0.0
-    b_fs: float | dolfin.Function | dolfin.Constant = 0.0
+    parameters: HolzapfelOgdenParameters = field(default_factory=default_parameters)
     use_subplus: bool = field(default=True, repr=False)
     use_heaviside: bool = field(default=True, repr=False)
 
@@ -129,19 +141,23 @@ class HolzapfelOgden(HyperElasticMaterial):
     )
 
     def __post_init__(self):
+        params = type(self).transversely_isotropic_parameters()
+        params.update(self.parameters)
+        self.parameters = params
         # Check that all values are positive
-        for attr in ["a", "b", "a_f", "b_f", "a_s", "a_s", "a_fs", "b_fs"]:
+
+        for name in self.parameters.keys():
             if not exceptions.check_value_greater_than(
-                getattr(self, attr),
+                self.parameters[name],
                 0.0,
                 inclusive=True,
             ):
                 raise exceptions.InvalidRangeError(
-                    name=attr,
+                    name=name,
                     expected_range=(0.0, np.inf),
                 )
-            if isinstance((v := getattr(self, attr)), float):
-                setattr(self, attr, dolfin.Constant(v))
+            if isinstance((v := self.parameters[name]), float):
+                self.parameters[name] = dolfin.Constant(v)
 
         if self.f0 is not None:
             self._I4f = partial(invariants.I4, a0=self.f0)
@@ -159,31 +175,22 @@ def __post_init__(self):
             self._I8fs = lambda x: 0.0
 
         self._W1_func = self._resolve_W1()
-        self._W4f_func = self._resolve_W4(a=self.a_f, b=self.b_f, required_attr="f0")
-        self._W4s_func = self._resolve_W4(a=self.a_s, b=self.b_s, required_attr="s0")
+        self._W4f_func = self._resolve_W4(
+            a=self.parameters["a_f"], b=self.parameters["b_f"], required_attr="f0"
+        )
+        self._W4s_func = self._resolve_W4(
+            a=self.parameters["a_s"], b=self.parameters["b_s"], required_attr="s0"
+        )
         self._W8fs_func = self._resolve_W8fs()
 
-    @property
-    def parameters(self) -> dict[str, float | dolfin.Constant | dolfin.Function]:
-        return {
-            "a": self.a,
-            "b": self.b,
-            "a_f": self.a_f,
-            "b_f": self.b_f,
-            "a_s": self.a_s,
-            "b_s": self.b_s,
-            "a_fs": self.a_fs,
-            "b_fs": self.b_fs,
-        }
-
     def _resolve_W1(self):
-        if exceptions.check_value_greater_than(self.a, 1e-10):
-            if exceptions.check_value_greater_than(self.b, 1e-10):
-                return lambda I1: (self.a / (2.0 * self.b)) * (
-                    ufl.exp(self.b * (I1 - 3)) - 1.0
-                )
+        if exceptions.check_value_greater_than(self.parameters["a"], 1e-10):
+            if exceptions.check_value_greater_than(self.parameters["b"], 1e-10):
+                return lambda I1: (
+                    self.parameters["a"] / (2.0 * self.parameters["b"])
+                ) * (ufl.exp(self.parameters["b"] * (I1 - 3)) - 1.0)
             else:
-                return lambda I1: (self.a / 2.0) * (I1 - 3)
+                return lambda I1: (self.parameters["a"] / 2.0) * (I1 - 3)
         else:
             return lambda I1: 0.0
 
@@ -209,25 +216,25 @@ def _resolve_W4(self, a, b, required_attr: str):
             return lambda I4: 0.0
 
     def _resolve_W8fs(self):
-        if exceptions.check_value_greater_than(self.a_fs, 1e-10):
+        if exceptions.check_value_greater_than(self.parameters["a_fs"], 1e-10):
             if self.f0 is None or self.s0 is None:
                 raise exceptions.MissingModelAttribute(
                     attr="f0 and/or s0",
                     model=type(self).__name__,
                 )
-            if exceptions.check_value_greater_than(self.b_fs, 1e-10):
+            if exceptions.check_value_greater_than(self.parameters["b_fs"], 1e-10):
                 return (
-                    lambda I8: self.a_fs
-                    / (2.0 * self.b_fs)
-                    * (ufl.exp(self.b_fs * I8**2) - 1.0)
+                    lambda I8: self.parameters["a_fs"]
+                    / (2.0 * self.parameters["b_fs"])
+                    * (ufl.exp(self.parameters["b_fs"] * I8**2) - 1.0)
                 )
             else:
-                return lambda I8: self.a_fs / 2.0 * I8**2
+                return lambda I8: self.parameters["a_fs"] / 2.0 * I8**2
         else:
             return lambda I8: 0.0
 
     @staticmethod
-    def transversely_isotropic_parameters() -> dict[str, float]:
+    def transversely_isotropic_parameters() -> HolzapfelOgdenParameters:
         """
         Material parameters for the Holzapfel Ogden model
         Taken from Table 1 row 3 in the main paper
@@ -245,7 +252,7 @@ def transversely_isotropic_parameters() -> dict[str, float]:
         }
 
     @staticmethod
-    def partly_orthotropic_parameters() -> dict[str, float]:
+    def partly_orthotropic_parameters() -> HolzapfelOgdenParameters:
         """
         Material parameters for the Holzapfel Ogden model
         Taken from Table 1 row 1 in the main paper
@@ -263,7 +270,7 @@ def partly_orthotropic_parameters() -> dict[str, float]:
         }
 
     @staticmethod
-    def orthotropic_parameters() -> dict[str, float]:
+    def orthotropic_parameters() -> HolzapfelOgdenParameters:
         """
         Material parameters for the Holzapfel Ogden model
         Taken from Table 1 row 2 in the main paper

src/pulse2/linear_elastic.py

@@ -1,15 +1,22 @@
+from __future__ import annotations
 from dataclasses import dataclass
 
 import dolfin
 import ufl_legacy as ufl
+from typing import TypedDict
 
 
 from . import exceptions
 from . import kinematics
 from .material_model import Material
 
 
-@dataclass(frozen=True, slots=True)
+class LinearElasticParameters(TypedDict):
+    E: float | dolfin.Function | dolfin.Constant
+    nu: float | dolfin.Function | dolfin.Constant
+
+
+@dataclass(slots=True)
 class LinearElastic(Material):
     """Linear elastic material
 
@@ -21,18 +28,13 @@ class LinearElastic(Material):
         Poisson's ratio
     """
 
-    E: float | dolfin.Function | dolfin.Constant
-    nu: float | dolfin.Function | dolfin.Constant
+    parameters: LinearElasticParameters
 
     def __post_init__(self):
         # The poisson ratio has to be between -1.0 and 0.5
-        if not exceptions.check_value_between(self.nu, -1.0, 0.5):
+        if not exceptions.check_value_between(self.parameters["nu"], -1.0, 0.5):
             raise exceptions.InvalidRangeError(name="mu", expected_range=(-1.0, 0.5))
 
-    @property
-    def parameters(self) -> dict[str, float | dolfin.Constant | dolfin.Function]:
-        return {"E": self.E, "nu": self.nu}
-
     def sigma(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         r"""Cauchy stress for linear elastic material
 
@@ -53,6 +55,7 @@ def sigma(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
         """
         e = kinematics.EngineeringStrain(F)
         I = kinematics.SecondOrderIdentity(F)
-        return (self.E / (1 + self.nu)) * (
-            e + (self.nu / (1 - 2 * self.nu)) * ufl.tr(e) * I
+        return (self.parameters["E"] / (1 + self.parameters["nu"])) * (
+            e
+            + (self.parameters["nu"] / (1 - 2 * self.parameters["nu"])) * ufl.tr(e) * I
         )

src/pulse2/material_model.py

@@ -1,6 +1,5 @@
 import abc
 
-import dolfin
 import ufl_legacy as ufl
 
 from . import kinematics
@@ -9,7 +8,7 @@
 class Material(abc.ABC):
     @property
     @abc.abstractmethod
-    def parameters(self) -> dict[str, dolfin.Constant | dolfin.Function]:
+    def parameters(self):
         ...
 
     @abc.abstractmethod

tests/test_cardiac_model.py

@@ -6,10 +6,9 @@
 
 
 def test_CardiacModel(mesh, u):
-    material_params = pulse2.HolzapfelOgden.transversely_isotropic_parameters()
     f0 = dolfin.Constant((1.0, 0.0, 0.0))
     s0 = dolfin.Constant((0.0, 1.0, 0.0))
-    material = pulse2.HolzapfelOgden(f0=f0, s0=s0, **material_params)
+    material = pulse2.HolzapfelOgden(f0=f0, s0=s0)
 
     active_model = pulse2.ActiveStress(f0)
     comp_model = pulse2.Compressible()

tests/test_itertarget.py

@@ -17,7 +17,7 @@ def problem(lvgeo):
     )
 
     material_params = pulse2.HolzapfelOgden.transversely_isotropic_parameters()
-    material = pulse2.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)
+    material = pulse2.HolzapfelOgden(f0=geo.f0, s0=geo.s0, parameters=material_params)
 
     Ta = dolfin.Constant(0.0)
     active_model = pulse2.ActiveStress(geo.f0, activation=Ta)