Source code for fenicsx_pulse.material_models.holzapfelogden

import logging
import typing
from dataclasses import dataclass, field
from functools import partial

import dolfinx
import numpy as np
import ufl

from .. import exceptions, functions, invariants
from ..material_model import HyperElasticMaterial
from ..units import Variable

Invariant = typing.Callable[[ufl.core.expr.Expr], ufl.core.expr.Expr]

logger = logging.getLogger(__name__)


def heaviside(x: ufl.core.expr.Expr, use_heaviside: bool) -> ufl.core.expr.Expr:
    if use_heaviside:
        return functions.heaviside(x)
    return ufl.as_ufl(1.0)


[docs] @dataclass(slots=True) class HolzapfelOgden(HyperElasticMaterial): r""" Orthotropic model by Holzapfel and Ogden Parameters ---------- f0: dolfinx.fem.Function | dolfinx.fem.Constant | None Function representing the direction of the fibers s0: dolfinx.fem.Function | dolfinx.fem.Constant | None Function representing the direction of the sheets a: float | dolfinx.fem.Function | dolfinx.fem.Constant Material parameter, by default 0.0 b: float | dolfinx.fem.Function | dolfinx.fem.Constant Material parameter, by default 0.0 a_f: float | dolfinx.fem.Function | dolfinx.fem.Constant Material parameter, by default 0.0 b_f: float | dolfinx.fem.Function | dolfinx.fem.Constant Material parameter, by default 0.0 a_s: float | dolfinx.fem.Function | dolfinx.fem.Constant Material parameter, by default 0.0 b_s: float | dolfinx.fem.Function | dolfinx.fem.Constant Material parameter, by default 0.0 a_fs: float | dolfinx.fem.Function | dolfinx.fem.Constant Material parameter, by default 0.0 b_fs: float | dolfinx.fem.Function | dolfinx.fem.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 Notes ----- Original model from Holzapfel and Ogden [1]_. The strain energy density function is given by .. math:: \Psi(I_1, I_{4\mathbf{f}_0}, I_{4\mathbf{s}_0}, I_{8\mathbf{f}_0\mathbf{s}_0}) = \frac{a}{2 b} \left( e^{ b (I_1 - 3)} -1 \right) + \frac{a_f}{2 b_f} \mathcal{H}(I_{4\mathbf{f}_0} - 1) \left( e^{ b_f (I_{4\mathbf{f}_0} - 1)_+^2} -1 \right) + \frac{a_s}{2 b_s} \mathcal{H}(I_{4\mathbf{s}_0} - 1) \left( e^{ b_s (I_{4\mathbf{s}_0} - 1)_+^2} -1 \right) + \frac{a_{fs}}{2 b_{fs}} \left( e^{ b_{fs} I_{8 \mathbf{f}_0 \mathbf{s}_0}^2} -1 \right) where .. math:: (x)_+ = \max\{x,0\} and .. math:: \mathcal{H}(x) = \begin{cases} 1, & \text{if $x > 0$} \\ 0, & \text{if $x \leq 0$} \end{cases} is the Heaviside function. .. [1] Holzapfel, Gerhard A., and Ray W. Ogden. "Constitutive modelling of passive myocardium: a structurally based framework for material characterization. "Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 367.1902 (2009): 3445-3475. """ f0: dolfinx.fem.Function | dolfinx.fem.Constant | None = None s0: dolfinx.fem.Function | dolfinx.fem.Constant | None = None a: Variable = field(default_factory=lambda: Variable(0.0, "kPa")) b: Variable = field(default_factory=lambda: Variable(0.0, "dimensionless")) a_f: Variable = field(default_factory=lambda: Variable(0.0, "kPa")) b_f: Variable = field(default_factory=lambda: Variable(0.0, "dimensionless")) a_s: Variable = field(default_factory=lambda: Variable(0.0, "kPa")) b_s: Variable = field(default_factory=lambda: Variable(0.0, "dimensionless")) a_fs: Variable = field(default_factory=lambda: Variable(0.0, "kPa")) 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) _W1_func: Invariant = field( init=False, repr=False, ) _W4f_func: Invariant = field( init=False, repr=False, ) _W4s_func: Invariant = field( init=False, repr=False, ) _W8fs_func: Invariant = field( init=False, repr=False, ) _I4f: Invariant = field( init=False, repr=False, ) _I4s: Invariant = field( init=False, repr=False, ) _I8fs: Invariant = field( init=False, repr=False, ) def __post_init__(self): # Check that all values are positive for attr in ["a", "b", "a_f", "b_f", "a_s", "a_s", "a_fs", "b_fs"]: p = getattr(self, attr) if not isinstance(p, Variable): unit = "dimensionless" if attr.startswith("b") else "kPa" if attr.startswith("a"): logger.warning("Setting %s to %s %s", attr, p, unit) p = Variable(p, unit) setattr(self, attr, p) if not exceptions.check_value_greater_than( getattr(self, attr).value, 0.0, inclusive=True, ): raise exceptions.InvalidRangeError( name=attr, expected_range=(0.0, np.inf), ) if self.f0 is not None: self._I4f = partial(invariants.I4, a0=self.f0) else: self._I4f = lambda x: 0.0 if self.s0 is not None: self._I4s = partial(invariants.I4, a0=self.s0) else: self._I4s = lambda x: 0.0 if self.f0 is not None and self.s0 is not None: self._I8fs = partial(invariants.I8, a0=self.f0, b0=self.s0) else: self._I8fs = lambda x: 0.0 a_f = self.a_f b_f = self.b_f a_s = self.a_s b_s = self.b_s self._W1_func = self._resolve_W1() self._W4f_func = self._resolve_W4(a=a_f, b=b_f, required_attr="f0") self._W4s_func = self._resolve_W4(a=a_s, b=b_s, required_attr="s0") self._W8fs_func = self._resolve_W8fs() def _resolve_W1(self) -> Invariant: a = self.a.to_base_units() b = self.b.to_base_units() if exceptions.check_value_greater_than(self.a.value, 1e-10): if exceptions.check_value_greater_than(self.b.value, 1e-10): return lambda I1: (a / (2.0 * b)) * (ufl.exp(b * (I1 - 3)) - 1.0) else: return lambda I1: (a / 2.0) * (I1 - 3) else: return lambda I1: 0.0 def _resolve_W4(self, a: Variable, b: Variable, required_attr: str) -> Invariant: subplus = functions.subplus if self.use_subplus else lambda x: x if exceptions.check_value_greater_than(a.value, 1e-10): a0 = getattr(self, required_attr) if a0 is None: raise exceptions.MissingModelAttribute( attr=required_attr, model=type(self).__name__, ) if exceptions.check_value_greater_than(b.value, 1e-10): return ( lambda I4: (a.to_base_units() / (2.0 * b.to_base_units())) * heaviside(I4 - 1, use_heaviside=self.use_heaviside) * (ufl.exp(b.to_base_units() * subplus(I4 - 1) ** 2) - 1.0) ) else: return ( lambda I4: (a.to_base_units() / 2.0) * heaviside(I4 - 1, use_heaviside=self.use_heaviside) * subplus(I4 - 1) ** 2 ) else: return lambda I4: 0.0 def _resolve_W8fs(self) -> Invariant: a_fs = self.a_fs.to_base_units() b_fs = self.b_fs.to_base_units() if exceptions.check_value_greater_than(self.a_fs.value, 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.value, 1e-10): return lambda I8: a_fs / (2.0 * b_fs) * (ufl.exp(b_fs * I8**2) - 1.0) else: return lambda I8: a_fs / 2.0 * I8**2 else: return lambda I8: 0.0
[docs] @staticmethod def transversely_isotropic_parameters() -> dict[str, Variable]: """ Material parameters for the Holzapfel Ogden model Taken from Table 1 row 3 in the main paper """ return { "a": Variable(2.280, "kPa"), "b": Variable(9.726, "dimensionless"), "a_f": Variable(1.685, "kPa"), "b_f": Variable(15.779, "dimensionless"), "a_s": Variable(0.0, "kPa"), "b_s": Variable(0.0, "dimensionless"), "a_fs": Variable(0.0, "kPa"), "b_fs": Variable(0.0, "dimensionless"), }
[docs] @staticmethod def partly_orthotropic_parameters() -> dict[str, Variable]: """ Material parameters for the Holzapfel Ogden model Taken from Table 1 row 1 in the main paper """ return { "a": Variable(0.057, "kPa"), "b": Variable(8.094, "dimensionless"), "a_f": Variable(21.503, "kPa"), "b_f": Variable(15.819, "dimensionless"), "a_s": Variable(6.841, "kPa"), "b_s": Variable(6.959, "dimensionless"), "a_fs": Variable(0.0, "kPa"), "b_fs": Variable(0.0, "dimensionless"), }
[docs] @staticmethod def orthotropic_parameters() -> dict[str, Variable]: """ Material parameters for the Holzapfel Ogden model Taken from Table 1 row 2 in the main paper """ return { "a": Variable(0.059, "kPa"), "b": Variable(8.023, "dimensionless"), "a_f": Variable(18.472, "kPa"), "b_f": Variable(16.026, "dimensionless"), "a_s": Variable(2.481, "kPa"), "b_s": Variable(11.120, "dimensionless"), "a_fs": Variable(0.216, "kPa"), "b_fs": Variable(11.436, "dimensionless"), }
def _W1(self, I1: ufl.core.expr.Expr) -> ufl.core.expr.Expr: return self._W1_func(I1) def _W4f(self, I4f: ufl.core.expr.Expr) -> ufl.core.expr.Expr: return self._W4f_func(I4f) def _W4s(self, I4s: ufl.core.expr.Expr) -> ufl.core.expr.Expr: return self._W4s_func(I4s) def _W8fs(self, I8fs: ufl.core.expr.Expr) -> ufl.core.expr.Expr: return self._W8fs_func(I8fs)
[docs] def strain_energy(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr: I1 = invariants.I1(F) I4f = self._I4f(F) I4s = self._I4s(F) I8fs = self._I8fs(F) return self._W1(I1) + self._W4f(I4f) + self._W4s(I4s) + self._W8fs(I8fs)
def __str__(self) -> str: return ( "a/2b (exp(b(I1 - 3)) - 1) + " "af/2bf H(I4f - 1) (exp(bf (I4f - 1)_+^2) - 1) + " "as/2bs H(I4s - 1) (exp(bs (I4s - 1)_+^2) - 1) + " "afs/2bfs (exp(bfs I8fs^2) - 1)" )