API Reference

Contents

API Reference#

kinematics#

pulse.kinematics.Cdev(C)[source]#

Deviatoric part of the right Cauchy-Green deformation tensor.

Parameters:

C (ufl.core.expr.Expr) – Right Cauchy-Green deformation tensor

Returns:

Deviatoric part of the right Cauchy-Green deformation tensor

Return type:

ufl.core.expr.Expr

pulse.kinematics.DeformationGradient(u)[source]#

Return deformation gradient from displacement.

Parameters:

  • u (dolfinx.fem.Function) – The displacement field

  • isochoric (bool, optional) – If true return the isochoric deformation gradient, by default False

Returns:

The deformation gradient

Return type:

ufl.core.expr.Expr

Notes

Given a displacement field \(\mathbf{u}\), the deformation gradient is given by

\[\mathbf{F} = \mathbf{I} + \nabla \mathbf{u}\]
pulse.kinematics.EngineeringStrain(F)[source]#

Engineering strain

\[\mathbf{\varepsilon} = \frac{1}{2}\left( \nabla u + (\nabla u)^T \right)\]
Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

Engineering strain

Return type:

ufl.core.expr.Expr

pulse.kinematics.Fdev(F)[source]#

Deviatoric part of the deformation gradient.

Parameters:

F (ufl.core.expr.Expr) – Deformation gradient

Returns:

Deviatoric part of the deformation gradient

Return type:

ufl.core.expr.Expr

pulse.kinematics.GreenLagrangeStrain(F)[source]#

Green-Lagrange strain tensor

\[\mathbf{E} = \frac{1}{2}\left( \mathbf{C} - \mathbf{I}\right)\]
Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

Green-Lagrange strain tensor

Return type:

ufl.core.expr.Expr

pulse.kinematics.InversePiolaTransform(A, F)[source]#

Push-forward of a rank two-tensor from the reference to the current configuration

Parameters:
Returns:

The push-forward

Return type:

ufl.core.expr.Expr

Notes

A push-forward is a transformation of a rank-2 tensor from the reference configuration to the current configuration. A common example is the push-forward of the First Piola-Kirchhoff stress tensor to the current configuration which yields the Cauchy stress tensor, i.e

\[\sigma = \frac{1}{J} \mathbf{P} \mathbf{F}^T\]
pulse.kinematics.IsochoricDeformationGradient(F)[source]#

Return the isochoric deformation gradient

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The isochoric deformation gradient

Return type:

ufl.core.expr.Expr

Notes

We can decompose the the deformation gradient multiplicatively into the volumetric (\(\mathbf{F}_{\mathrm{vol}}\)) and isochoric (\(\mathbf{F}_{\mathrm{iso}}\)) components

\[\mathbf{F} = \mathbf{F}_{\mathrm{iso}} \cdot \mathbf{F}_{\mathrm{vol}}\]

such that \(\mathrm{det}(\mathbf{F}_{\mathrm{iso}}) = 1\). In this case, we can work out that

\[\mathbf{F}_{\mathrm{vol}} = J^{1/3}\mathbf{I}\]

and consequently

\[\mathbf{F}_{\mathrm{vol}} = J^{-1/3}\mathbf{F}\]

with \(J = \mathrm{det}(\mathbf{F})\).

pulse.kinematics.Jacobian(F)[source]#

Determinant of the deformation gradient

\[J = \mathrm{det}(\mathbf{F})\]
Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The jacobian

Return type:

ufl.core.expr.Expr

pulse.kinematics.LeftCauchyGreen(F)[source]#

Left Cauchy-Green deformation tensor

\[\mathbf{C} = \mathbf{F}\mathbf{F}^T\]
Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

Left Cauchy-Green deformation tensor

Return type:

ufl.core.expr.Expr

pulse.kinematics.PiolaTransform(A, F)[source]#

Pull-back of a two-tensor from the current to the reference configuration

Parameters:
Returns:

The pull-back

Return type:

ufl.core.expr.Expr

Notes

A pull-back is a transformation of a rank-2 tensor from the current configuration to the reference configuration. A common example is the pull-back of the Cauchy stress tensor to the reference configuration which yields the First Piola-Kirchhoff stress tensor, i.e

\[\mathbf{P} = J \sigma \mathbf{F}^{-T}\]
pulse.kinematics.RightCauchyGreen(F)[source]#

Right Cauchy-Green deformation tensor

\[\mathbf{C} = \mathbf{F}^T\mathbf{F}\]
Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

Right Cauchy-Green deformation tensor

Return type:

ufl.core.expr.Expr

pulse.kinematics.SecondOrderIdentity(F)[source]#

Return identity with same dimension as input

Parameters:

F (Expr)

Return type:

Expr

invariants#

This module defines invariants of the deformation gradient.

The deformation gradient \(\mathbf{F}\) is a tensor field that describes the deformation of a material. The invariants of the deformation gradient are scalar quantities that are invariant under rigid body motion.

The invariants are defined as follows:

  • First principal invariant

    \(I_1 = \text{tr}(\mathbf{C})\)

  • Second principal invariant

    \(I_2 = \left( I_1^2 - \text{tr}(\mathbf{C}\cdot\mathbf{C})\right)\)

  • Third principal invariant

    \(I_3 = \text{det}(\mathbf{C})\)

  • Fourth quasi invariant

    \(I_{4\mathbf{a_0}} = \mathbf{C}\mathbf{a}_0 \cdot \mathbf{a}_0\)

  • Fifth quasi invariant

    \(I_{5\mathbf{a_0}} = \mathbf{C}\mathbf{a}_0 \cdot \mathbf{C}\mathbf{a}_0\)

  • Eight quasi invariant

    \(I_{8\mathbf{a_0}\mathbf{b_0}} = \mathbf{F}\mathbf{a}_0 \cdot \mathbf{F}\mathbf{b}_0\)

where \(\mathbf{C} = \mathbf{F}^T\mathbf{F}\) is the right Cauchy-Green tensor.

pulse.invariants.I1(C)[source]#

First principal invariant

\[I_1 = \mathrm{tr}(\mathbf{C})\]
Parameters:

C (ufl.core.expr.Expr) – The deformation gradient

Returns:

First principal invariant

Return type:

ufl.core.expr.Expr

pulse.invariants.I2(C)[source]#

Second principal invariant

\[I_2 = \left( I_1^2 - \mathrm{tr}(\mathbf{C}\cdot\mathbf{C})\right)\]
Parameters:
Returns:

Second principal invariant

Return type:

ufl.core.expr.Expr

pulse.invariants.I3(C)[source]#

Third principal invariant

\[I_3 = \mathrm{det}(\mathbf{C})\]
Parameters:

C (ufl.core.expr.Expr) – The deformation gradient

Returns:

Third principal invariant

Return type:

ufl.core.expr.Expr

pulse.invariants.I4(C, a0)[source]#

Fourth quasi invariant

\[I_{4\mathbf{a_0}} = \mathbf{C}\mathbf{a}_0 \cdot \mathbf{a}_0\]
Parameters:
Returns:

Fourth quasi invariant in the direction a0

Return type:

ufl.core.expr.Expr

pulse.invariants.I5(C, a0)[source]#

Fifth quasi invariant

\[I_{5\mathbf{a_0}} = \mathbf{C}\mathbf{a}_0 \cdot \mathbf{C}\mathbf{a}_0\]
Parameters:
Returns:

Fifth quasi invariant in the direction a0

Return type:

ufl.core.expr.Expr

pulse.invariants.I8(C, a0, b0)[source]#

Eight quasi invariant

\[I_{8\mathbf{a_0}\mathbf{b_0}} = \mathbf{F}\mathbf{a}_0 \cdot \mathbf{F}\mathbf{b}_0\]
Parameters:

  • C (ufl.core.expr.Expr) – The deformation gradient

  • a0 (FuncConst) – Some direction

  • b0 (FuncConst) – Another direction

Returns:

Eight quasi invariant in the direction a0

Return type:

ufl.core.expr.Expr

material_model#

This module defines the material model interface and some common material models.

The material model describes the mechanical behavior of a material. The material model is used to compute the stress tensor given the deformation gradient.

The material model interface defines two methods:

  • sigma(F): The Cauchy stress tensor

  • P(F): The first Piola-Kirchhoff stress tensor

The sigma method computes the Cauchy stress tensor given the deformation gradient. The P method computes the first Piola-Kirchhoff stress tensor given the deformation gradient.

The HyperElasticMaterial class is a base class for hyperelastic material models. Hyperelastic materials are materials that have a strain energy density function that depends only on the deformation gradient. The strain_energy method computes the strain energy density function given the deformation gradient.

class pulse.material_model.HyperElasticMaterial[source]#
P(F, dev=True)[source]#

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:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

Notes

For a hyperelastic material model with strain energy density function \(\Psi = \Psi(\mathbf{F})\), the first Piola-Kirchhoff stress tensor is given by

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}}\]
S(C, dev=True)[source]#

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:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

sigma(F)[source]#

Cauchy stress tensor

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

abstractmethod strain_energy(F)[source]#

Strain energy density function

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The strain energy density

Return type:

ufl.core.expr.Expr

class pulse.material_model.Material[source]#
P(F)[source]#

First Piola-Kirchhoff stress tensor

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

abstractmethod sigma(F)[source]#

Cauchy stress tensor

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

class pulse.material_models.holzapfelogden.HolzapfelOgden(f0=None, s0=None, a=<factory>, b=<factory>, a_f=<factory>, b_f=<factory>, a_s=<factory>, b_s=<factory>, a_fs=<factory>, b_fs=<factory>, use_subplus=True, use_heaviside=True)[source]#

Orthotropic model by Holzapfel and Ogden

Parameters:

Notes

Original model from Holzapfel and Ogden [HO09]. The strain energy density function is given by

\[\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

\[(x)_+ = \max\{x,0\}\]

and

\[\begin{split}\mathcal{H}(x) = \begin{cases} 1, & \text{if $x > 0$} \\ 0, & \text{if $x \leq 0$} \end{cases}\end{split}\]

is the Heaviside function.

P(F, dev=True)#

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:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

Notes

For a hyperelastic material model with strain energy density function \(\Psi = \Psi(\mathbf{F})\), the first Piola-Kirchhoff stress tensor is given by

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}}\]
S(C, dev=True)#

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:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

static orthotropic_parameters()[source]#

Material parameters for the Holzapfel Ogden model Taken from Table 1 row 2 in the main paper

Return type:

HolzapfelOgdenParameters

static partly_orthotropic_parameters()[source]#

Material parameters for the Holzapfel Ogden model Taken from Table 1 row 1 in the main paper

Return type:

HolzapfelOgdenParameters

sigma(F)#

Cauchy stress tensor

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

strain_energy(C)[source]#

Strain energy density function

Parameters:

C (ufl.core.expr.Expr) – Right Cauchy-Green deformation tensor

Returns:

The strain energy density

Return type:

ufl.core.expr.Expr

static transversely_isotropic_parameters()[source]#

Material parameters for the Holzapfel Ogden model Taken from Table 1 row 3 in the main paper

Return type:

HolzapfelOgdenParameters

class pulse.material_models.holzapfelogden.HolzapfelOgdenParameters[source]#
clear() → None.  Remove all items from D.#
copy() → a shallow copy of D#
classmethod fromkeys(iterable, value=None, /)#

Create a new dictionary with keys from iterable and values set to value.

get(key, default=None, /)#

Return the value for key if key is in the dictionary, else default.

items() → a set-like object providing a view on D's items#
keys() → a set-like object providing a view on D's keys#
pop(k[, d]) → v, remove specified key and return the corresponding value.#

If the key is not found, return the default if given; otherwise, raise a KeyError.

popitem()#

Remove and return a (key, value) pair as a 2-tuple.

Pairs are returned in LIFO (last-in, first-out) order. Raises KeyError if the dict is empty.

setdefault(key, default=None, /)#

Insert key with a value of default if key is not in the dictionary.

Return the value for key if key is in the dictionary, else default.

update([E, ]**F) → None.  Update D from dict/iterable E and F.#

If E is present and has a .keys() method, then does: for k in E: D[k] = E[k] If E is present and lacks a .keys() method, then does: for k, v in E: D[k] = v In either case, this is followed by: for k in F: D[k] = F[k]

values() → an object providing a view on D's values#
class pulse.material_models.guccione.Guccione(f0=None, s0=None, n0=None, C=<factory>, bf=<factory>, bt=<factory>, bfs=<factory>)[source]#

Transversely isotropic exponential model by Guccione et al.

Parameters:

Notes

Original model from Guccione et al. [GMW91]. The strain energy density function is given by

\[\Psi = \frac{C}{2} \left( \mathrm{exp}^{Q} - 1 \right)\]

where

\[Q = b_f E_{11}^2 + b_t \left( E_{22}^2 + E_{33}^2 + E_{23}^2 + E_{32}^2 \right) + b_{fs} \left( E_{12}^2 + E_{21}^2 + E_{13}^2 + E_{31}^2 \right)\]
P(F, dev=True)#

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:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

Notes

For a hyperelastic material model with strain energy density function \(\Psi = \Psi(\mathbf{F})\), the first Piola-Kirchhoff stress tensor is given by

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}}\]
S(C, dev=True)#

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:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

is_isotropic()[source]#

Return True if the material is isotropic.

Return type:

bool

sigma(F)#

Cauchy stress tensor

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

strain_energy(C)[source]#

Strain energy density function

Parameters:
Returns:

The strain energy density

Return type:

ufl.core.expr.Expr

class pulse.material_models.guccione.GuccioneParameters[source]#
clear() → None.  Remove all items from D.#
copy() → a shallow copy of D#
classmethod fromkeys(iterable, value=None, /)#

Create a new dictionary with keys from iterable and values set to value.

get(key, default=None, /)#

Return the value for key if key is in the dictionary, else default.

items() → a set-like object providing a view on D's items#
keys() → a set-like object providing a view on D's keys#
pop(k[, d]) → v, remove specified key and return the corresponding value.#

If the key is not found, return the default if given; otherwise, raise a KeyError.

popitem()#

Remove and return a (key, value) pair as a 2-tuple.

Pairs are returned in LIFO (last-in, first-out) order. Raises KeyError if the dict is empty.

setdefault(key, default=None, /)#

Insert key with a value of default if key is not in the dictionary.

Return the value for key if key is in the dictionary, else default.

update([E, ]**F) → None.  Update D from dict/iterable E and F.#

If E is present and has a .keys() method, then does: for k in E: D[k] = E[k] If E is present and lacks a .keys() method, then does: for k, v in E: D[k] = v In either case, this is followed by: for k in F: D[k] = F[k]

values() → an object providing a view on D's values#
class pulse.material_models.saint_venant_kirchhoff.SaintVenantKirchhoff(mu, lmbda)[source]#

Class for Saint Venant-Kirchhoff material

Let

\[\epsilon = \frac{1}{2} ( \nabla u + \nabla u^T + \nabla u + \nabla u^T )\]

Then

Parameters:
P(F, dev=True)#

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:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

Notes

For a hyperelastic material model with strain energy density function \(\Psi = \Psi(\mathbf{F})\), the first Piola-Kirchhoff stress tensor is given by

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}}\]
S(C, dev=True)#

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:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

sigma(F)#

Cauchy stress tensor

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

strain_energy(C)[source]#

Strain energy density function

Parameters:
Returns:

The strain energy density

Return type:

ufl.core.expr.Expr

active_model#

This module defines the ActiveModel class which is an abstract class for active models. Active models are used to incorporate active stress or active strain in the material model.

The ActiveModel class defines two methods:

  • Fe(F): Transforming the deformation gradient to an active deformation gradient

  • strain_energy(F): Active strain energy density function

The Fe method transforms the deformation gradient to an active deformation gradient. For example in the active strain approach we perform a multiplicative decomposition of the deformation gradient into an elastic and an active part, i.e

\[\mathbf{F} = \mathbf{F}_e \mathbf{F}_a\]

In which case the active model can be incorporated by transforming the full deformation gradient into a pure elastic component.

The strain_energy method defines the active strain energy density function. For example in the active stress approach, the active stress is added as an extra stress component

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}} + \mathbf{P}_a\]

where \(\mathbf{P}_a\). Now we can instead rewrite this as a total strain energy by considering the following equation

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}} = \frac{\partial \Psi_p}{\partial \mathbf{F}} + \frac{\partial \Psi_a}{\partial \mathbf{F}}\]

where \(\Psi_p\) is the passive (classical) strain energy density function and \(\Psi_a\) is the corresponding active strain energy density function.

The Passive class is a simple active model with no active component. This model could for example be used if you want to use a pure passive model.

class pulse.active_model.ActiveModel[source]#
abstractmethod Fe(F)[source]#

Method to transforming the deformation gradient to an an active deformation gradient. For example in the active strain approach we perform a multiplicative decomposition of the deformation gradient into an elastic and an active part, i.e

\[\mathbf{F} = \mathbf{F}_e \mathbf{F}_a\]

In which case the active model can be incorporated by transforming the full deformation gradient into a pure elastic component

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The elastic deformation gradient

Return type:

ufl.core.expr.Expr

P(F)[source]#

First Piola-Kirchhoff stress tensor for the active model.

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

S(C)[source]#

Second Piola-Kirchhoff stress tensor for the active model.

The active model is evaluated on the full right Cauchy-Green tensor, never on its isochoric part: see the note in pulse.cardiac_model.CardiacModel.strain_energy().

Parameters:

C (ufl.core.expr.Expr) – The right Cauchy-Green deformation tensor

Returns:

The second Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

evaluate_at_end_of_step: bool = False#

Evaluate this model’s stress at the end of the time step rather than at the generalized-alpha \(\alpha_f\) point.

Set this (as a class attribute on a subclass, or as an instance attribute, which overrides it) when the stress depends on state advanced over the step – a stretch rate, or ODE states integrated using it – rather than on the instantaneous deformation alone. Assembling such a model at the alpha_f-interpolated configuration would feed it a blended stretch and a rate scaled by \(1 - \alpha_f\), not the ones it actually advanced with. pulse.problem.DynamicProblem reads this flag and, when it is set, assembles the active stress at the true end-of-step displacement instead – the same treatment it already gives the cavity constraint and \(J - 1\). A pulse.problem.StaticProblem has no alpha_f point to differ from, so the flag has no effect there.

register(u)[source]#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Parameters:

u (Function)

Return type:

None

abstractmethod strain_energy(C)[source]#

Active strain energy function. For example in the active stress approach, the active stress is added as an extra stress component

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}} + \mathbf{P}_a\]

where \(\mathbf{P}_a\). Now we can instead rewrite this as a total strain energy by considering the following equation

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}} = \frac{\partial \Psi_p}{\partial \mathbf{F}} + \frac{\partial \Psi_a}{\partial \mathbf{F}}\]

where \(\Psi_p\) is the passive (classical) strain energy density function and \(\Psi_a\) is the corresponding active strain energy density function.

Parameters:

C (ufl.core.expr.Expr) – The right Cauchy-Green deformation tensor

Returns:

The active strain energy density function

Return type:

ufl.core.expr.Expr

class pulse.active_model.Passive[source]#

Active model with no active component. This model could for example be used if you want to use a pure passive model.

Fe(F)[source]#

Method to transforming the deformation gradient to an an active deformation gradient. For example in the active strain approach we perform a multiplicative decomposition of the deformation gradient into an elastic and an active part, i.e

\[\mathbf{F} = \mathbf{F}_e \mathbf{F}_a\]

In which case the active model can be incorporated by transforming the full deformation gradient into a pure elastic component

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The elastic deformation gradient

Return type:

ufl.core.expr.Expr

P(F)#

First Piola-Kirchhoff stress tensor for the active model.

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

S(C)#

Second Piola-Kirchhoff stress tensor for the active model.

The active model is evaluated on the full right Cauchy-Green tensor, never on its isochoric part: see the note in pulse.cardiac_model.CardiacModel.strain_energy().

Parameters:

C (ufl.core.expr.Expr) – The right Cauchy-Green deformation tensor

Returns:

The second Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

evaluate_at_end_of_step: bool = False#

Evaluate this model’s stress at the end of the time step rather than at the generalized-alpha \(\alpha_f\) point.

Set this (as a class attribute on a subclass, or as an instance attribute, which overrides it) when the stress depends on state advanced over the step – a stretch rate, or ODE states integrated using it – rather than on the instantaneous deformation alone. Assembling such a model at the alpha_f-interpolated configuration would feed it a blended stretch and a rate scaled by \(1 - \alpha_f\), not the ones it actually advanced with. pulse.problem.DynamicProblem reads this flag and, when it is set, assembles the active stress at the true end-of-step displacement instead – the same treatment it already gives the cavity constraint and \(J - 1\). A pulse.problem.StaticProblem has no alpha_f point to differ from, so the flag has no effect there.

register(u)#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Parameters:

u (Function)

Return type:

None

strain_energy(C)[source]#

Active strain energy function. For example in the active stress approach, the active stress is added as an extra stress component

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}} + \mathbf{P}_a\]

where \(\mathbf{P}_a\). Now we can instead rewrite this as a total strain energy by considering the following equation

\[\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}} = \frac{\partial \Psi_p}{\partial \mathbf{F}} + \frac{\partial \Psi_a}{\partial \mathbf{F}}\]

where \(\Psi_p\) is the passive (classical) strain energy density function and \(\Psi_a\) is the corresponding active strain energy density function.

Parameters:

C (ufl.core.expr.Expr) – The right Cauchy-Green deformation tensor

Returns:

The active strain energy density function

Return type:

ufl.core.expr.Expr

This module contains the active stress model for the cardiac mechanics problem. The active stress model is used to describe the active contraction of the heart. The active stress model is used to compute the active stress given the deformation gradient.

class pulse.active_stress.ActiveStress(f0, activation=<factory>, s0=None, n0=None, T_ref=1.0, eta=0.0, isotropy=ActiveStressModels.transversely, formulation=ActiveStressFormulation.invariant)[source]#

Active stress model

f0: dolfinx.fem.Function | dolfinx.fem.Constant

The cardiac fiber direction

activation: dolfinx.fem.Function | dolfinx.fem.Constant | None

A function or constant representing the activation. If not provided a constant will be created.

s0: dolfinx.fem.Function | dolfinx.fem.Constant | None

The sheets orientation. Only needed for orthotropic active stress models

n0: dolfinx.fem.Function | dolfinx.fem.Constant | None

The sheet-normal orientation. Only needed for orthotropic active stress models

T_ref: float = 1.0

Reference active stress, by default 1.0

eta: float = 0.0

Amount of transverse active stress, by default 0.0. A value of zero means that all active stress is along the fiber direction. If the value is 1.0 then all active stress will be in the transverse direction.

isotropy: ActiveStressModels

What kind of active stress model to use, by default ‘transversely’

formulation: ActiveStressFormulation

Which active-stress convention to use, by default ‘invariant’. The two differ by a factor of the fiber stretch, so this changes results rather than just their derivation – see ActiveStressFormulation for how to choose.

Parameters:
Fe(F)[source]#

Method to transforming the deformation gradient to an an active deformation gradient. For example in the active strain approach we perform a multiplicative decomposition of the deformation gradient into an elastic and an active part, i.e

\[\mathbf{F} = \mathbf{F}_e \mathbf{F}_a\]

In which case the active model can be incorporated by transforming the full deformation gradient into a pure elastic component

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The elastic deformation gradient

Return type:

ufl.core.expr.Expr

P(F)#

First Piola-Kirchhoff stress tensor for the active model.

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

S(C)[source]#

Second Piola-Kirchhoff stress tensor for the active stress model.

Parameters:

C (ufl.core.expr.Expr) – The right Cauchy-Green deformation tensor

Returns:

The second Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

property Ta: Expr#

The active stress

evaluate_at_end_of_step: bool = False#

Evaluate this model’s stress at the end of the time step rather than at the generalized-alpha \(\alpha_f\) point.

Set this (as a class attribute on a subclass, or as an instance attribute, which overrides it) when the stress depends on state advanced over the step – a stretch rate, or ODE states integrated using it – rather than on the instantaneous deformation alone. Assembling such a model at the alpha_f-interpolated configuration would feed it a blended stretch and a rate scaled by \(1 - \alpha_f\), not the ones it actually advanced with. pulse.problem.DynamicProblem reads this flag and, when it is set, assembles the active stress at the true end-of-step displacement instead – the same treatment it already gives the cavity constraint and \(J - 1\). A pulse.problem.StaticProblem has no alpha_f point to differ from, so the flag has no effect there.

register(u)#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Parameters:

u (Function)

Return type:

None

strain_energy(C)[source]#

Active strain energy density

Parameters:

C (ufl.core.expr.Expr) – The right Cauchy-Green deformation tensor

Returns:

The active strain energy density

Return type:

ufl.core.expr.Expr

Raises:

NotImplementedError – _description_

class pulse.active_stress.ActiveStressFormulation(*values)[source]#

Which power of the fiber stretch the active energy is linear in.

Both formulations put the active stress along the fiber, and differ only in what \(T_a\) is taken to multiply – by a factor of the fiber stretch \(\lambda\). Both appear in the literature, so the choice must be explicit rather than implied.

invariant

\(\Psi_a = \frac{1}{2} T_a (I_{4f} - 1)\), giving \(\mathbf{S}_a = T_a\, f_0 \otimes f_0\) and \(\mathbf{P}_a = T_a\, \mathbf{F} f_0 \otimes f_0\). Adding \(T_a f_0 \otimes f_0\) to the second Piola stress is the most widespread form of the active stress approach in cardiac mechanics, and it remains the default here for that reason and for backwards compatibility. Note the fibre traction it delivers scales with the stretch: \(|\mathbf{P}_a f_0| = T_a \lambda\).

stretch

\(\Psi_a = T_a (\lambda - 1)\) with \(\lambda = \sqrt{I_{4f}}\), giving \(\mathbf{S}_a = \frac{T_a}{\lambda} f_0 \otimes f_0\) and \(\mathbf{P}_a = T_a \frac{\mathbf{F} f_0 \otimes f_0} {|\mathbf{F} f_0|}\). This is the convention used by Regazzoni & Quarteroni [RQ21], and the one StabilizedActiveStress is derived in. Choose it when \(T_a\) comes from a force-generation model whose active stiffness is defined as \(\partial \dot{T_a}/\partial \dot\lambda\), so that tension and stiffness refer to the same kinematic variable.

capitalize()#

Return a capitalized version of the string.

More specifically, make the first character have upper case and the rest lower case.

casefold()#

Return a version of the string suitable for caseless comparisons.

center(width, fillchar=' ', /)#

Return a centered string of length width.

Padding is done using the specified fill character (default is a space).

count(sub[, start[, end]]) → int#

Return the number of non-overlapping occurrences of substring sub in string S[start:end]. Optional arguments start and end are interpreted as in slice notation.

encode(encoding='utf-8', errors='strict')#

Encode the string using the codec registered for encoding.

encoding

The encoding in which to encode the string.

errors

The error handling scheme to use for encoding errors. The default is ‘strict’ meaning that encoding errors raise a UnicodeEncodeError. Other possible values are ‘ignore’, ‘replace’ and ‘xmlcharrefreplace’ as well as any other name registered with codecs.register_error that can handle UnicodeEncodeErrors.

endswith(suffix[, start[, end]]) → bool#

Return True if S ends with the specified suffix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. suffix can also be a tuple of strings to try.

expandtabs(tabsize=8)#

Return a copy where all tab characters are expanded using spaces.

If tabsize is not given, a tab size of 8 characters is assumed.

find(sub[, start[, end]]) → int#

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

format(*args, **kwargs) → str#

Return a formatted version of S, using substitutions from args and kwargs. The substitutions are identified by braces (‘{’ and ‘}’).

format_map(mapping) → str#

Return a formatted version of S, using substitutions from mapping. The substitutions are identified by braces (‘{’ and ‘}’).

index(sub[, start[, end]]) → int#

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Raises ValueError when the substring is not found.

isalnum()#

Return True if the string is an alpha-numeric string, False otherwise.

A string is alpha-numeric if all characters in the string are alpha-numeric and there is at least one character in the string.

isalpha()#

Return True if the string is an alphabetic string, False otherwise.

A string is alphabetic if all characters in the string are alphabetic and there is at least one character in the string.

isascii()#

Return True if all characters in the string are ASCII, False otherwise.

ASCII characters have code points in the range U+0000-U+007F. Empty string is ASCII too.

isdecimal()#

Return True if the string is a decimal string, False otherwise.

A string is a decimal string if all characters in the string are decimal and there is at least one character in the string.

isdigit()#

Return True if the string is a digit string, False otherwise.

A string is a digit string if all characters in the string are digits and there is at least one character in the string.

isidentifier()#

Return True if the string is a valid Python identifier, False otherwise.

Call keyword.iskeyword(s) to test whether string s is a reserved identifier, such as “def” or “class”.

islower()#

Return True if the string is a lowercase string, False otherwise.

A string is lowercase if all cased characters in the string are lowercase and there is at least one cased character in the string.

isnumeric()#

Return True if the string is a numeric string, False otherwise.

A string is numeric if all characters in the string are numeric and there is at least one character in the string.

isprintable()#

Return True if the string is printable, False otherwise.

A string is printable if all of its characters are considered printable in repr() or if it is empty.

isspace()#

Return True if the string is a whitespace string, False otherwise.

A string is whitespace if all characters in the string are whitespace and there is at least one character in the string.

istitle()#

Return True if the string is a title-cased string, False otherwise.

In a title-cased string, upper- and title-case characters may only follow uncased characters and lowercase characters only cased ones.

isupper()#

Return True if the string is an uppercase string, False otherwise.

A string is uppercase if all cased characters in the string are uppercase and there is at least one cased character in the string.

join(iterable, /)#

Concatenate any number of strings.

The string whose method is called is inserted in between each given string. The result is returned as a new string.

Example: ‘.’.join([‘ab’, ‘pq’, ‘rs’]) -> ‘ab.pq.rs’

ljust(width, fillchar=' ', /)#

Return a left-justified string of length width.

Padding is done using the specified fill character (default is a space).

lower()#

Return a copy of the string converted to lowercase.

lstrip(chars=None, /)#

Return a copy of the string with leading whitespace removed.

If chars is given and not None, remove characters in chars instead.

static maketrans()#

Return a translation table usable for str.translate().

If there is only one argument, it must be a dictionary mapping Unicode ordinals (integers) or characters to Unicode ordinals, strings or None. Character keys will be then converted to ordinals. If there are two arguments, they must be strings of equal length, and in the resulting dictionary, each character in x will be mapped to the character at the same position in y. If there is a third argument, it must be a string, whose characters will be mapped to None in the result.

partition(sep, /)#

Partition the string into three parts using the given separator.

This will search for the separator in the string. If the separator is found, returns a 3-tuple containing the part before the separator, the separator itself, and the part after it.

If the separator is not found, returns a 3-tuple containing the original string and two empty strings.

removeprefix(prefix, /)#

Return a str with the given prefix string removed if present.

If the string starts with the prefix string, return string[len(prefix):]. Otherwise, return a copy of the original string.

removesuffix(suffix, /)#

Return a str with the given suffix string removed if present.

If the string ends with the suffix string and that suffix is not empty, return string[:-len(suffix)]. Otherwise, return a copy of the original string.

replace(old, new, count=-1, /)#

Return a copy with all occurrences of substring old replaced by new.

count

Maximum number of occurrences to replace. -1 (the default value) means replace all occurrences.

If the optional argument count is given, only the first count occurrences are replaced.

rfind(sub[, start[, end]]) → int#

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

rindex(sub[, start[, end]]) → int#

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Raises ValueError when the substring is not found.

rjust(width, fillchar=' ', /)#

Return a right-justified string of length width.

Padding is done using the specified fill character (default is a space).

rpartition(sep, /)#

Partition the string into three parts using the given separator.

This will search for the separator in the string, starting at the end. If the separator is found, returns a 3-tuple containing the part before the separator, the separator itself, and the part after it.

If the separator is not found, returns a 3-tuple containing two empty strings and the original string.

rsplit(sep=None, maxsplit=-1)#

Return a list of the substrings in the string, using sep as the separator string.

sep

The separator used to split the string.

When set to None (the default value), will split on any whitespace character (including n r t f and spaces) and will discard empty strings from the result.

maxsplit

Maximum number of splits. -1 (the default value) means no limit.

Splitting starts at the end of the string and works to the front.

rstrip(chars=None, /)#

Return a copy of the string with trailing whitespace removed.

If chars is given and not None, remove characters in chars instead.

split(sep=None, maxsplit=-1)#

Return a list of the substrings in the string, using sep as the separator string.

sep

The separator used to split the string.

When set to None (the default value), will split on any whitespace character (including n r t f and spaces) and will discard empty strings from the result.

maxsplit

Maximum number of splits. -1 (the default value) means no limit.

Splitting starts at the front of the string and works to the end.

Note, str.split() is mainly useful for data that has been intentionally delimited. With natural text that includes punctuation, consider using the regular expression module.

splitlines(keepends=False)#

Return a list of the lines in the string, breaking at line boundaries.

Line breaks are not included in the resulting list unless keepends is given and true.

startswith(prefix[, start[, end]]) → bool#

Return True if S starts with the specified prefix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. prefix can also be a tuple of strings to try.

strip(chars=None, /)#

Return a copy of the string with leading and trailing whitespace removed.

If chars is given and not None, remove characters in chars instead.

swapcase()#

Convert uppercase characters to lowercase and lowercase characters to uppercase.

title()#

Return a version of the string where each word is titlecased.

More specifically, words start with uppercased characters and all remaining cased characters have lower case.

translate(table, /)#

Replace each character in the string using the given translation table.

table

Translation table, which must be a mapping of Unicode ordinals to Unicode ordinals, strings, or None.

The table must implement lookup/indexing via __getitem__, for instance a dictionary or list. If this operation raises LookupError, the character is left untouched. Characters mapped to None are deleted.

upper()#

Return a copy of the string converted to uppercase.

zfill(width, /)#

Pad a numeric string with zeros on the left, to fill a field of the given width.

The string is never truncated.

class pulse.active_stress.ActiveStressModels(*values)[source]#
capitalize()#

Return a capitalized version of the string.

More specifically, make the first character have upper case and the rest lower case.

casefold()#

Return a version of the string suitable for caseless comparisons.

center(width, fillchar=' ', /)#

Return a centered string of length width.

Padding is done using the specified fill character (default is a space).

count(sub[, start[, end]]) → int#

Return the number of non-overlapping occurrences of substring sub in string S[start:end]. Optional arguments start and end are interpreted as in slice notation.

encode(encoding='utf-8', errors='strict')#

Encode the string using the codec registered for encoding.

encoding

The encoding in which to encode the string.

errors

The error handling scheme to use for encoding errors. The default is ‘strict’ meaning that encoding errors raise a UnicodeEncodeError. Other possible values are ‘ignore’, ‘replace’ and ‘xmlcharrefreplace’ as well as any other name registered with codecs.register_error that can handle UnicodeEncodeErrors.

endswith(suffix[, start[, end]]) → bool#

Return True if S ends with the specified suffix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. suffix can also be a tuple of strings to try.

expandtabs(tabsize=8)#

Return a copy where all tab characters are expanded using spaces.

If tabsize is not given, a tab size of 8 characters is assumed.

find(sub[, start[, end]]) → int#

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

format(*args, **kwargs) → str#

Return a formatted version of S, using substitutions from args and kwargs. The substitutions are identified by braces (‘{’ and ‘}’).

format_map(mapping) → str#

Return a formatted version of S, using substitutions from mapping. The substitutions are identified by braces (‘{’ and ‘}’).

index(sub[, start[, end]]) → int#

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Raises ValueError when the substring is not found.

isalnum()#

Return True if the string is an alpha-numeric string, False otherwise.

A string is alpha-numeric if all characters in the string are alpha-numeric and there is at least one character in the string.

isalpha()#

Return True if the string is an alphabetic string, False otherwise.

A string is alphabetic if all characters in the string are alphabetic and there is at least one character in the string.

isascii()#

Return True if all characters in the string are ASCII, False otherwise.

ASCII characters have code points in the range U+0000-U+007F. Empty string is ASCII too.

isdecimal()#

Return True if the string is a decimal string, False otherwise.

A string is a decimal string if all characters in the string are decimal and there is at least one character in the string.

isdigit()#

Return True if the string is a digit string, False otherwise.

A string is a digit string if all characters in the string are digits and there is at least one character in the string.

isidentifier()#

Return True if the string is a valid Python identifier, False otherwise.

Call keyword.iskeyword(s) to test whether string s is a reserved identifier, such as “def” or “class”.

islower()#

Return True if the string is a lowercase string, False otherwise.

A string is lowercase if all cased characters in the string are lowercase and there is at least one cased character in the string.

isnumeric()#

Return True if the string is a numeric string, False otherwise.

A string is numeric if all characters in the string are numeric and there is at least one character in the string.

isprintable()#

Return True if the string is printable, False otherwise.

A string is printable if all of its characters are considered printable in repr() or if it is empty.

isspace()#

Return True if the string is a whitespace string, False otherwise.

A string is whitespace if all characters in the string are whitespace and there is at least one character in the string.

istitle()#

Return True if the string is a title-cased string, False otherwise.

In a title-cased string, upper- and title-case characters may only follow uncased characters and lowercase characters only cased ones.

isupper()#

Return True if the string is an uppercase string, False otherwise.

A string is uppercase if all cased characters in the string are uppercase and there is at least one cased character in the string.

join(iterable, /)#

Concatenate any number of strings.

The string whose method is called is inserted in between each given string. The result is returned as a new string.

Example: ‘.’.join([‘ab’, ‘pq’, ‘rs’]) -> ‘ab.pq.rs’

ljust(width, fillchar=' ', /)#

Return a left-justified string of length width.

Padding is done using the specified fill character (default is a space).

lower()#

Return a copy of the string converted to lowercase.

lstrip(chars=None, /)#

Return a copy of the string with leading whitespace removed.

If chars is given and not None, remove characters in chars instead.

static maketrans()#

Return a translation table usable for str.translate().

If there is only one argument, it must be a dictionary mapping Unicode ordinals (integers) or characters to Unicode ordinals, strings or None. Character keys will be then converted to ordinals. If there are two arguments, they must be strings of equal length, and in the resulting dictionary, each character in x will be mapped to the character at the same position in y. If there is a third argument, it must be a string, whose characters will be mapped to None in the result.

partition(sep, /)#

Partition the string into three parts using the given separator.

This will search for the separator in the string. If the separator is found, returns a 3-tuple containing the part before the separator, the separator itself, and the part after it.

If the separator is not found, returns a 3-tuple containing the original string and two empty strings.

removeprefix(prefix, /)#

Return a str with the given prefix string removed if present.

If the string starts with the prefix string, return string[len(prefix):]. Otherwise, return a copy of the original string.

removesuffix(suffix, /)#

Return a str with the given suffix string removed if present.

If the string ends with the suffix string and that suffix is not empty, return string[:-len(suffix)]. Otherwise, return a copy of the original string.

replace(old, new, count=-1, /)#

Return a copy with all occurrences of substring old replaced by new.

count

Maximum number of occurrences to replace. -1 (the default value) means replace all occurrences.

If the optional argument count is given, only the first count occurrences are replaced.

rfind(sub[, start[, end]]) → int#

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

rindex(sub[, start[, end]]) → int#

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Raises ValueError when the substring is not found.

rjust(width, fillchar=' ', /)#

Return a right-justified string of length width.

Padding is done using the specified fill character (default is a space).

rpartition(sep, /)#

Partition the string into three parts using the given separator.

This will search for the separator in the string, starting at the end. If the separator is found, returns a 3-tuple containing the part before the separator, the separator itself, and the part after it.

If the separator is not found, returns a 3-tuple containing two empty strings and the original string.

rsplit(sep=None, maxsplit=-1)#

Return a list of the substrings in the string, using sep as the separator string.

sep

The separator used to split the string.

When set to None (the default value), will split on any whitespace character (including n r t f and spaces) and will discard empty strings from the result.

maxsplit

Maximum number of splits. -1 (the default value) means no limit.

Splitting starts at the end of the string and works to the front.

rstrip(chars=None, /)#

Return a copy of the string with trailing whitespace removed.

If chars is given and not None, remove characters in chars instead.

split(sep=None, maxsplit=-1)#

Return a list of the substrings in the string, using sep as the separator string.

sep

The separator used to split the string.

When set to None (the default value), will split on any whitespace character (including n r t f and spaces) and will discard empty strings from the result.

maxsplit

Maximum number of splits. -1 (the default value) means no limit.

Splitting starts at the front of the string and works to the end.

Note, str.split() is mainly useful for data that has been intentionally delimited. With natural text that includes punctuation, consider using the regular expression module.

splitlines(keepends=False)#

Return a list of the lines in the string, breaking at line boundaries.

Line breaks are not included in the resulting list unless keepends is given and true.

startswith(prefix[, start[, end]]) → bool#

Return True if S starts with the specified prefix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. prefix can also be a tuple of strings to try.

strip(chars=None, /)#

Return a copy of the string with leading and trailing whitespace removed.

If chars is given and not None, remove characters in chars instead.

swapcase()#

Convert uppercase characters to lowercase and lowercase characters to uppercase.

title()#

Return a version of the string where each word is titlecased.

More specifically, words start with uppercased characters and all remaining cased characters have lower case.

translate(table, /)#

Replace each character in the string using the given translation table.

table

Translation table, which must be a mapping of Unicode ordinals to Unicode ordinals, strings, or None.

The table must implement lookup/indexing via __getitem__, for instance a dictionary or list. If this operation raises LookupError, the character is left untouched. Characters mapped to None are deleted.

upper()#

Return a copy of the string converted to uppercase.

zfill(width, /)#

Pad a numeric string with zeros on the left, to fill a field of the given width.

The string is never truncated.

class pulse.active_stress.FrankStarlingActiveStress(f0, activation=<factory>, s0=None, n0=None, T_ref=1.0, eta=0.0, isotropy=ActiveStressModels.transversely, formulation=ActiveStressFormulation.invariant, amp_min=0.0, amp_max=1.0, stretch_threshold=0.85, stretch_optimal=1.15)[source]#

Active stress model incorporating the Frank-Starling mechanism. Multiplies the baseline time-dependent activation by a stretch-dependent factor.

Parameters:

  • amp_min (float, optional) – The minimum amplification factor, by default 0.0.

  • amp_max (float, optional) – The maximum amplification factor, by default 1.0.

  • stretch_threshold (float, optional) – The stretch ratio below which the active force remains at its minimum, by default 0.85.

  • stretch_optimal (float, optional) – The optimal stretch ratio where the active force reaches its maximum plateau, by default 1.15.

  • f0 (Function | Constant)

  • activation (Variable)

  • s0 (Function | Constant | None)

  • n0 (Function | Constant | None)

  • T_ref (Constant | float)

  • eta (Constant | float)

  • isotropy (ActiveStressModels)

  • formulation (ActiveStressFormulation)

Fe(F)#

Method to transforming the deformation gradient to an an active deformation gradient. For example in the active strain approach we perform a multiplicative decomposition of the deformation gradient into an elastic and an active part, i.e

\[\mathbf{F} = \mathbf{F}_e \mathbf{F}_a\]

In which case the active model can be incorporated by transforming the full deformation gradient into a pure elastic component

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The elastic deformation gradient

Return type:

ufl.core.expr.Expr

P(F)#

First Piola-Kirchhoff stress tensor for the active model.

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

S(C)#

Second Piola-Kirchhoff stress tensor for the active stress model.

Parameters:

C (ufl.core.expr.Expr) – The right Cauchy-Green deformation tensor

Returns:

The second Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

property Ta: Expr#

Overrides the base active tension property from ActiveStress. The parent class methods (like S and stress_tensor) will automatically use this dynamically scaled active tension.

Returns:

The total active tension (baseline activation * multiplier)

Return type:

ufl.core.expr.Expr

evaluate_at_end_of_step: bool = False#

Evaluate this model’s stress at the end of the time step rather than at the generalized-alpha \(\alpha_f\) point.

Set this (as a class attribute on a subclass, or as an instance attribute, which overrides it) when the stress depends on state advanced over the step – a stretch rate, or ODE states integrated using it – rather than on the instantaneous deformation alone. Assembling such a model at the alpha_f-interpolated configuration would feed it a blended stretch and a rate scaled by \(1 - \alpha_f\), not the ones it actually advanced with. pulse.problem.DynamicProblem reads this flag and, when it is set, assembles the active stress at the true end-of-step displacement instead – the same treatment it already gives the cavity constraint and \(J - 1\). A pulse.problem.StaticProblem has no alpha_f point to differ from, so the flag has no effect there.

frank_starling_multiplier()[source]#

Class method wrapper that evaluates the standalone Frank-Starling multiplier function using the registered displacement and material properties.

Returns:

A symbolic UFL expression of the multiplier.

Return type:

ufl.core.expr.Expr

Raises:

ValueError – If the displacement field u has not been registered yet.

register(u)[source]#

Registers the displacement field into the material model. This must be called before the active stress is evaluated so the model can calculate the dynamic stretch.

Parameters:

u (ufl.core.expr.Expr) – The displacement vector field to register.

strain_energy(C)#

Active strain energy density

Parameters:

C (ufl.core.expr.Expr) – The right Cauchy-Green deformation tensor

Returns:

The active strain energy density

Return type:

ufl.core.expr.Expr

Raises:

NotImplementedError – _description_

class pulse.active_stress.StabilizedActiveStress(f0, activation=<factory>, active_stiffness=<factory>, lmbda_prev=None)[source]#

Active stress with the consistent stabilization term of Regazzoni & Quarteroni, for use when \(T_a\) comes from an external force-generation solver.

Why you probably want this#

The usual way to drive ActiveStress is to advance some cell-level contraction model, write its tension into activation, and solve mechanics with that value held fixed. This is a segregated (staggered) scheme, and it has a failure mode that is easy to hit and hard to diagnose: whenever the active stiffness of the tissue exceeds its passive stiffness – routine in contracting myocardium – the scheme develops non-physical oscillations in \(T_a\) and \(\lambda\). Regazzoni & Quarteroni [RQ21] show it is then not merely inaccurate but not convergent, its amplification factor tending to \(-K_a/K_p < -1\) as \(\Delta t \to 0\). Reducing the time step makes it worse, so the problem cannot be tuned away.

The cause is that a staggered scheme treats active tension as a dead load over the mechanics solve, when physically it is a population of crossbridges behaving as springs. Restoring that gives

\[\mathbf{P}_{act} = \left[T_a + K_a(\lambda - \lambda_{prev})\right] \frac{\mathbf{F} f_0 \otimes f_0}{|\mathbf{F} f_0|}\]

which is the gradient of

\[\Psi_a = T_a (\lambda - \lambda_{prev}) + \tfrac{1}{2} K_a (\lambda - \lambda_{prev})^2\]

and is what this class implements. The extra term is \(\mathcal{O}(\Delta t)\) and vanishes in the limit, so the scheme remains consistent with the same continuous problem – it is a numerical device, not a change of model – while becoming unconditionally stable.

Usage#

Each time step, in this order:

  1. advance the force-generation model using \(\lambda_{prev}\),

  2. assign the resulting tension and stiffness to activation and active_stiffness,

  3. solve mechanics,

  4. call update_prev() with the new displacement.

Step 4 matters: \(\lambda_{prev}\) must be the same stretch that was fed to the force-generation model in step 1. If the two drift apart the added term is no longer a consistent perturbation and can itself destabilize the solve.

param f0:

The cardiac fiber direction

type f0:

dolfinx.fem.Function | dolfinx.fem.Constant

param activation:

The active tension \(T_a\), from the force-generation model

type activation:

Variable

param active_stiffness:

The active stiffness \(K_a = \partial \dot{T_a} / \partial \dot\lambda\), from the same model, in the same units as activation and per unit of the same stretch variable. Setting it to zero recovers the plain staggered scheme – useful for demonstrating the instability, not for production.

type active_stiffness:

Variable

param lmbda_prev:

Fiber stretch at the previous time step. Defaults to a constant 1.0, i.e. the reference configuration. Pass a Function (and use update_prev()) for anything beyond a single step.

type lmbda_prev:

dolfinx.fem.Function | dolfinx.fem.Constant | None

Notes

Unlike ActiveStress this takes no T_ref or eta. A reference scaling applied to \(T_a\) but not \(K_a\) would silently break the consistency of the stabilization, and there is no accepted transverse generalization of an energy written in \(\lambda\). Scale both quantities before assigning them instead.

Fe(F)[source]#

Method to transforming the deformation gradient to an an active deformation gradient. For example in the active strain approach we perform a multiplicative decomposition of the deformation gradient into an elastic and an active part, i.e

\[\mathbf{F} = \mathbf{F}_e \mathbf{F}_a\]

In which case the active model can be incorporated by transforming the full deformation gradient into a pure elastic component

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The elastic deformation gradient

Return type:

ufl.core.expr.Expr

property Ka: Expr#

Active stiffness, in base units.

P(F)#

First Piola-Kirchhoff stress tensor for the active model.

Parameters:

F (ufl.core.expr.Expr) – The deformation gradient

Returns:

The first Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

S(C)[source]#

Second Piola-Kirchhoff stress,

\[\mathbf{S} = \frac{T_a + K_a \Delta\lambda}{\lambda} f_0 \otimes f_0\]

Given in closed form rather than by differentiating strain_energy(); the two agree, which test_stabilized_active_stress.py checks.

Parameters:

C (Expr)

Return type:

Expr

property Ta: Expr#

Active tension, in base units.

dlmbda(C)[source]#

Increment in fiber stretch since the previous time step.

Parameters:

C (Expr)

Return type:

Expr

evaluate_at_end_of_step: bool = False#

Evaluate this model’s stress at the end of the time step rather than at the generalized-alpha \(\alpha_f\) point.

Set this (as a class attribute on a subclass, or as an instance attribute, which overrides it) when the stress depends on state advanced over the step – a stretch rate, or ODE states integrated using it – rather than on the instantaneous deformation alone. Assembling such a model at the alpha_f-interpolated configuration would feed it a blended stretch and a rate scaled by \(1 - \alpha_f\), not the ones it actually advanced with. pulse.problem.DynamicProblem reads this flag and, when it is set, assembles the active stress at the true end-of-step displacement instead – the same treatment it already gives the cavity constraint and \(J - 1\). A pulse.problem.StaticProblem has no alpha_f point to differ from, so the flag has no effect there.

register(u)#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Parameters:

u (Function)

Return type:

None

strain_energy(C)[source]#

\(\Psi_a = T_a \Delta\lambda + \frac{1}{2} K_a \Delta\lambda^2\)

Parameters:

C (Expr)

Return type:

Expr

update_prev(u)[source]#

Record the fiber stretch of displacement u as \(\lambda_{prev}\).

Call once per time step, after the mechanics solve. Requires lmbda_prev to be a Function; a Constant cannot hold a spatially varying stretch.

Parameters:

u (Function)

Return type:

None

Parameters:
pulse.active_stress.compute_frank_starling_multiplier(u, f0, amp_min, amp_max, stretch_threshold, stretch_optimal)[source]#

Computes a stretch-dependent scalar multiplier for active tension to model the Frank-Starling mechanism using a piecewise linear ascending limb.

Parameters:

  • u (ufl.core.expr.Expr) – The macroscopic displacement vector field.

  • f0 (ufl.core.expr.Expr) – The reference fiber direction vector field.

  • amp_min (float) – The minimum amplification factor (used for tissue at resting or compressed lengths).

  • amp_max (float) – The maximum amplification factor (used for tissue at or beyond the optimal stretch length).

  • stretch_threshold (float) – The stretch ratio below which the active force remains at its minimum.

  • stretch_optimal (float) – The optimal stretch ratio where the active force reaches its maximum plateau.

Returns:

A symbolic UFL expression representing the spatial multiplier field g(lambda).

Return type:

ufl.core.expr.Expr

Notes

Mathematical Formulation:

Let the right Cauchy-Green deformation tensor be \(\mathbf{C} = \mathbf{F}^T \mathbf{F}\) where \(\mathbf{F} = \mathbf{I} + \nabla \mathbf{u}\) is the deformation gradient.

The local fiber stretch \(\lambda\) is computed as:

\[\lambda = \sqrt{\mathbf{f}_0 \cdot (\mathbf{C} \mathbf{f}_0)}\]

The multiplier \(g(\lambda)\) is defined as a piecewise function:

\[\begin{split}g(\lambda) = \begin{cases} a_{\min} & \text{if } \lambda \le \lambda_{\text{threshold}} \\ a_{\min} + m (\lambda - \lambda_{\text{threshold}}) & \text{if } \lambda_{\text{threshold}} < \lambda \le \lambda_{\text{opt}} \\ a_{\max} & \text{if } \lambda > \lambda_{\text{opt}} \end{cases}\end{split}\]

where the slope \(m\) is calculated as:

\[m = \frac{a_{\max} - a_{\min}}{\lambda_{\text{opt}} - \lambda_{\text{threshold}}}\]
pulse.active_stress.fiber_stretch(C, f0)[source]#

Stretch along the fiber direction, \(\lambda = \sqrt{f_0 \cdot C f_0}\).

Equal to \(|\mathbf{F} f_0|\), and to 1 in the reference configuration.

Parameters:
Return type:

Expr

pulse.active_stress.stretch_active_stress(Ta, C, f0)[source]#

Second Piola-Kirchhoff stress for stretch_active_stress_strain_energy(),

\[\mathbf{S} = \frac{T_a}{\lambda} f_0 \otimes f_0\]
Parameters:
pulse.active_stress.stretch_active_stress_strain_energy(Ta, C, f0)[source]#

Active strain energy that is linear in the fiber stretch,

\[W = T_a (\lambda - 1), \qquad \lambda = \sqrt{I_{4f}}\]

whose second Piola-Kirchhoff stress is \(T_a f_0 \otimes f_0 / \lambda\) and whose first Piola-Kirchhoff stress is therefore \(T_a \mathbf{F} f_0 \otimes f_0 / |\mathbf{F} f_0|\) – the normalization used by Regazzoni & Quarteroni [RQ21].

Compare transversely_active_stress_strain_energy(), which is linear in \(I_{4f} = \lambda^2\) instead and hence differs by a factor of \(\lambda\) in the resulting stress.

Parameters:
pulse.active_stress.transversely_active_stress(Ta, f0, eta=0.0)[source]#

Return the Cauchy stress tensor for the active stress model when activation is only working along the fibers, with a possible transverse component defined by \(\eta\) with \(\eta = 0\) meaning that all active stress is along the fiber and \(\eta = 1\) meaning that all active stress is in the transverse direction. The Cauchy stress tensor is given by

\[\sigma = T_a \left( I_{4f} - 1 + \eta ((I_1 - 3) - (I_{4f} - 1)) \right) f_0\]
Parameters:

  • Ta (dolfinx.fem.Function or dolfinx.fem.Constant) – A scalar function representing the magnitude of the active stress. With this (invariant) formulation the resulting fibre traction is \(|\mathbf{P}_a f_0| = T_a \lambda\); see ActiveStressFormulation.

  • f0 (dolfinx.fem.Function) – A vector function representing the direction of the active stress

  • eta (float) – Amount of active stress in the transverse direction (relative to f0)

pulse.active_stress.transversely_active_stress_strain_energy(Ta, C, f0, eta=0.0)[source]#

Return active strain energy when activation is only working along the fibers, with a possible transverse component defined by \(\eta\) with \(\eta = 0\) meaning that all active stress is along the fiber and \(\eta = 1\) meaning that all active stress is in the transverse direction. The active strain energy is given by

\[W = \frac{1}{2} T_a \left( I_{4f} - 1 + \eta ((I_1 - 3) - (I_{4f} - 1)) \right)\]
Parameters:

  • Ta (dolfinx.fem.Function or dolfinx.fem.Constant) – A scalar function representing the magnitude of the active stress. Note that with this (invariant) formulation the resulting fibre traction is \(|\mathbf{P}_a f_0| = T_a \lambda\), not \(T_a\) – it is ActiveStressFormulation stretch under which \(T_a\) is itself the first Piola fibre traction. Which one to supply depends on what your activation model’s \(T_a\) was calibrated to mean; see ActiveStressFormulation.

  • C (ufl.Form) – The right Cauchy-Green deformation tensor

  • f0 (dolfinx.fem.Function) – A vector function representing the direction of the active stress

  • eta (float) – Amount of active stress in the transverse direction (relative to f0)

compressibility#

This module defines compressibility models for the material models. We define two compressibility models: Incompressible and Compressible.

An incompressible material is a material that does not change its volume under deformation. The volume change is described by the Jacobian \(J = \det(F)\), where \(F\) is the deformation gradient.

The Incompressible model is defined by the strain energy density function \(\Psi = p (J - 1)\), where \(p\) is a function representing the Lagrange multiplier. The Compressible model is defined by the strain energy density function \(\Psi = \kappa (J \ln(J) - J + 1)\), where \(\kappa\) is a material parameter representing the bulk modulus. Higher values of \(\kappa\) correspond to more incompressible material.

class pulse.compressibility.Compressibility[source]#

Base class for compressibility models.

P(F)[source]#

First Piola-Kirchhoff stress tensor for the compressibility model.

Parameters:

F (Expr)

Return type:

Expr

S(C)[source]#

Cauchy stress tensor for the compressibility model.

Parameters:

C (Expr)

Return type:

Expr

abstractmethod is_compressible()[source]#

Returns True if the material model is compressible.

Return type:

bool

register(*args, **kwargs)[source]#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Return type:

None

abstractmethod strain_energy(C)[source]#

Strain energy density function

Parameters:

C (Expr)

Return type:

Expr

class pulse.compressibility.Compressible(kappa=<factory>)[source]#

Compressible material model

Strain energy density function is given by

\[\Psi = \kappa (J \ln(J) - J + 1)\]
Parameters:

kappa (Variable)

P(F)#

First Piola-Kirchhoff stress tensor for the compressibility model.

Parameters:

F (Expr)

Return type:

Expr

S(C)#

Cauchy stress tensor for the compressibility model.

Parameters:

C (Expr)

Return type:

Expr

is_compressible()[source]#

Returns True if the material model is compressible.

Return type:

bool

register(*args, **kwargs)#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Return type:

None

strain_energy(C)[source]#

Strain energy density function

Parameters:

C (Expr)

Return type:

Expr

class pulse.compressibility.Compressible2(kappa=<factory>)[source]#

Compressible material model

Strain energy density function is given by

\[\Psi = \kappa / 4 (J^2 - 1 - 2 \ln(J))\]
Parameters:

kappa (Variable)

P(F)#

First Piola-Kirchhoff stress tensor for the compressibility model.

Parameters:

F (Expr)

Return type:

Expr

S(C)#

Cauchy stress tensor for the compressibility model.

Parameters:

C (Expr)

Return type:

Expr

is_compressible()#

Returns True if the material model is compressible.

Return type:

bool

register(*args, **kwargs)#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Return type:

None

strain_energy(C)[source]#

Strain energy density function

Parameters:

C (Expr)

Return type:

Expr

class pulse.compressibility.Compressible3(kappa=<factory>)[source]#

Compressible material model used in Usyk et al. [ULM02]

Strain energy density function is given by

\[\Psi = \kappa / 2 (J - 1) * \ln(J)\]
Parameters:

kappa (Variable)

P(F)#

First Piola-Kirchhoff stress tensor for the compressibility model.

Parameters:

F (Expr)

Return type:

Expr

S(C)#

Cauchy stress tensor for the compressibility model.

Parameters:

C (Expr)

Return type:

Expr

is_compressible()#

Returns True if the material model is compressible.

Return type:

bool

register(*args, **kwargs)#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Return type:

None

strain_energy(C)[source]#

Strain energy density function

Parameters:

C (Expr)

Return type:

Expr

class pulse.compressibility.Incompressible[source]#

Incompressible material model

Strain energy density function is given by

\[\Psi = p (J - 1)\]
P(F)#

First Piola-Kirchhoff stress tensor for the compressibility model.

Parameters:

F (Expr)

Return type:

Expr

S(C)#

Cauchy stress tensor for the compressibility model.

Parameters:

C (Expr)

Return type:

Expr

is_compressible()[source]#

Returns True if the material model is compressible.

Return type:

bool

register(p)[source]#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Parameters:

p (Function | None)

Return type:

None

strain_energy(C)[source]#

Strain energy density function

Parameters:

C (Expr)

Return type:

Expr

viscoelasticity#

class pulse.viscoelasticity.NoneViscoElasticity[source]#
P(F_dot=None)#

First Piola-Kirchhoff stress for the viscoelasticity model.

Parameters:

F_dot (Expr | None)

Return type:

Expr

S(C_dot)#

Cauchy stress tensor for the viscoelasticity model.

Parameters:

C_dot (ufl.core.expr.Expr) – The time derivative of the deformation gradient (strain rate)

Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

strain_energy(C_dot)[source]#

Strain energy density function.

Parameters:

C_dot (ufl.core.expr.Expr) – The time derivative of the deformation gradient (strain rate)

Returns:

The strain energy density function

Return type:

ufl.Form

class pulse.viscoelasticity.ViscoElasticity[source]#
P(F_dot=None)[source]#

First Piola-Kirchhoff stress for the viscoelasticity model.

Parameters:

F_dot (Expr | None)

Return type:

Expr

S(C_dot)[source]#

Cauchy stress tensor for the viscoelasticity model.

Parameters:

C_dot (ufl.core.expr.Expr) – The time derivative of the deformation gradient (strain rate)

Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

abstractmethod strain_energy(C_dot)[source]#

Strain energy density function.

Parameters:

C_dot (ufl.core.expr.Expr) – The time derivative of the deformation gradient (strain rate)

Returns:

The strain energy density function

Return type:

ufl.Form

class pulse.viscoelasticity.Viscous(eta: pulse.units.Variable = <factory>)[source]#
Parameters:

eta (Variable)

P(F_dot=None)#

First Piola-Kirchhoff stress for the viscoelasticity model.

Parameters:

F_dot (Expr | None)

Return type:

Expr

S(C_dot)#

Cauchy stress tensor for the viscoelasticity model.

Parameters:

C_dot (ufl.core.expr.Expr) – The time derivative of the deformation gradient (strain rate)

Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

strain_energy(C_dot)[source]#

Strain energy density function.

Parameters:

C_dot (ufl.core.expr.Expr) – The time derivative of the deformation gradient (strain rate)

Returns:

The strain energy density function

Return type:

ufl.Form

cardiac_model#

This module defines the cardiac model.

The cardiac model is a combination of a material model, an active model, and a compressibility model.

class pulse.cardiac_model.ActiveModel(*args, **kwargs)[source]#

An active model contributes an extra stress, on the full deformation.

Unlike HyperElasticMaterial, there is no dev argument: active tension along a fiber is not an isochoric quantity, and the volumetric response is the compressibility model’s job. See CardiacModel.strain_energy().

register(u)[source]#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Parameters:

u (Function)

Return type:

None

class pulse.cardiac_model.CardiacModel(material: pulse.cardiac_model.HyperElasticMaterial, active: pulse.cardiac_model.ActiveModel, compressibility: pulse.cardiac_model.Compressibility, viscoelasticity: pulse.cardiac_model.ViscoElasticity = <factory>)[source]#
Parameters:
P(F, F_dot=None)[source]#

First Piola-Kirchhoff stress for the cardiac model.

Parameters:
Return type:

Expr

S(C, C_dot=None, *, active=True)[source]#

Second Piola-Kirchhoff stress for the cardiac model.

Parameters:
Returns:

The second Piola-Kirchhoff stress tensor

Return type:

ufl.core.expr.Expr

sigma(F, F_dot=None)[source]#

Cauchy stress tensor

Parameters:
Returns:

The Cauchy stress tensor

Return type:

ufl.core.expr.Expr

strain_energy(C, C_dot=None)[source]#

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:

The total strain energy density

Return type:

ufl.core.expr.Expr

Notes

The isochoric split is applied to the passive material only. The active model sees the full \(\mathbf{C}\), so that a prescribed active tension means what it says along the fiber rather than being reduced by a factor \(J^{-1/3}\) and picking up a spurious isotropic term from the chain rule through \(\bar{\mathbf{C}}\).

class pulse.cardiac_model.Compressibility(*args, **kwargs)[source]#
register(p)[source]#

Register a virtual subclass of an ABC.

Returns the subclass, to allow usage as a class decorator.

Parameters:

p (Function | None)

Return type:

None

class pulse.cardiac_model.HyperElasticMaterial(*args, **kwargs)[source]#
class pulse.cardiac_model.ViscoElasticity(*args, **kwargs)[source]#

geometry#

class pulse.geometry.CardiacGeometriesObject(*args, **kwargs)[source]#
pulse.geometry.DEFAULT_QUADRATURE_DEGREE = 4#

Quadrature degree of dx and ds when metadata names none, as the CLI’s geometry.quadrature_degree. Without one, UFL estimates the degree from each integrand, which for the exponentials, logarithms and inverses of cardiac models is very high: on a 4x4x4 cube a 4-iteration Newton solve took 98 s instead of 0.15 s. Pass metadata={"quadrature_degree": -1} to get that estimate anyway.

class pulse.geometry.Geometry(*, mesh: dolfinx.mesh.Mesh, boundaries: Sequence[pulse.geometry.Marker] = (), metadata: dict[str, typing.Any]=<factory>, facet_tags: dolfinx.mesh.MeshTags | None = None, markers: dict[str, tuple[int, int]]=<factory>)[source]#
Parameters:

  • mesh (Mesh)

  • boundaries (Sequence[Marker])

  • metadata (dict[str, Any])

  • facet_tags (MeshTags | None)

  • markers (dict[str, tuple[int, int]])

deform(u)[source]#

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.

Return type:

None

class pulse.geometry.HeartGeometry(*, mesh: dolfinx.mesh.Mesh, boundaries: Sequence[pulse.geometry.Marker] = (), metadata: dict[str, typing.Any]=<factory>, facet_tags: dolfinx.mesh.MeshTags | None = None, markers: dict[str, tuple[int, int]]=<factory>)[source]#
Parameters:

  • mesh (Mesh)

  • boundaries (Sequence[Marker])

  • metadata (dict[str, Any])

  • facet_tags (MeshTags | None)

  • markers (dict[str, tuple[int, int]])

base_center(base='BASE', u=None, dtype=<class 'numpy.float64'>)[source]#

Return the normal of the base

Parameters:

  • base (str, optional) – Marker for the base, by default “BASE”

  • u (dolfinx.fem.Function | None, optional) – Displacement field, by default None

Returns:

Normal of the base

Return type:

npt.NDArray[np.float64]

base_center_form(base='BASE', u=None)[source]#

Return the normal of the base

Parameters:

  • base (str, optional) – Marker for the base, by default “BASE”

  • u (dolfinx.fem.Function | None, optional) – Displacement field, by default None

Returns:

Normal of the base

Return type:

npt.NDArray[np.float64]

deform(u)#

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.

Return type:

None

volume(marker, u=None)[source]#

Return the volume of the cavity for a given marker

Parameters:

  • marker (str) – Marker for the surface of the cavity

  • u (dolfinx.fem.Function | None, optional) – Optional displacement field, by default None

Returns:

Volume of the cavity

Return type:

float

Raises:

exceptions.MarkerNotFoundError – If the marker is not found in the geometry

volume_form(u=None)[source]#

Return the form for the volume of the cavity for a given marker

Parameters:

u (dolfinx.fem.Function | None, optional) – Optional displacement field, by default None

Returns:

The form for the volume of the cavity

Return type:

dolfinx.fem.forms.Form

Raises:

exceptions.MarkerNotFoundError – If the marker is not found in the geometry

class pulse.geometry.Marker(name, marker, dim, locator)[source]#
Parameters:

  • name (str)

  • marker (int)

  • dim (int)

  • locator (Callable[[ndarray[tuple[Any, ...], dtype[float64]]], bool])

count(value, /)#

Return number of occurrences of value.

dim: int#

Alias for field number 2

index(value, start=0, stop=9223372036854775807, /)#

Return first index of value.

Raises ValueError if the value is not present.

locator: Callable[[ndarray[tuple[Any, ...], dtype[float64]]], bool]#

Alias for field number 3

marker: int#

Alias for field number 1

name: str#

Alias for field number 0

problem#

pulse.problem.BACKWARD_EULER_STENCIL = (1.0, -1.0, 0.0)#

Coefficients of (y, y_old, y_prev) in dt * dy/dt; see StaticProblem._circulation_form().

class pulse.problem.BaseBC(*values)[source]#

Base boundary condition

capitalize()#

Return a capitalized version of the string.

More specifically, make the first character have upper case and the rest lower case.

casefold()#

Return a version of the string suitable for caseless comparisons.

center(width, fillchar=' ', /)#

Return a centered string of length width.

Padding is done using the specified fill character (default is a space).

count(sub[, start[, end]]) → int#

Return the number of non-overlapping occurrences of substring sub in string S[start:end]. Optional arguments start and end are interpreted as in slice notation.

encode(encoding='utf-8', errors='strict')#

Encode the string using the codec registered for encoding.

encoding

The encoding in which to encode the string.

errors

The error handling scheme to use for encoding errors. The default is ‘strict’ meaning that encoding errors raise a UnicodeEncodeError. Other possible values are ‘ignore’, ‘replace’ and ‘xmlcharrefreplace’ as well as any other name registered with codecs.register_error that can handle UnicodeEncodeErrors.

endswith(suffix[, start[, end]]) → bool#

Return True if S ends with the specified suffix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. suffix can also be a tuple of strings to try.

expandtabs(tabsize=8)#

Return a copy where all tab characters are expanded using spaces.

If tabsize is not given, a tab size of 8 characters is assumed.

find(sub[, start[, end]]) → int#

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

format(*args, **kwargs) → str#

Return a formatted version of S, using substitutions from args and kwargs. The substitutions are identified by braces (‘{’ and ‘}’).

format_map(mapping) → str#

Return a formatted version of S, using substitutions from mapping. The substitutions are identified by braces (‘{’ and ‘}’).

index(sub[, start[, end]]) → int#

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Raises ValueError when the substring is not found.

isalnum()#

Return True if the string is an alpha-numeric string, False otherwise.

A string is alpha-numeric if all characters in the string are alpha-numeric and there is at least one character in the string.

isalpha()#

Return True if the string is an alphabetic string, False otherwise.

A string is alphabetic if all characters in the string are alphabetic and there is at least one character in the string.

isascii()#

Return True if all characters in the string are ASCII, False otherwise.

ASCII characters have code points in the range U+0000-U+007F. Empty string is ASCII too.

isdecimal()#

Return True if the string is a decimal string, False otherwise.

A string is a decimal string if all characters in the string are decimal and there is at least one character in the string.

isdigit()#

Return True if the string is a digit string, False otherwise.

A string is a digit string if all characters in the string are digits and there is at least one character in the string.

isidentifier()#

Return True if the string is a valid Python identifier, False otherwise.

Call keyword.iskeyword(s) to test whether string s is a reserved identifier, such as “def” or “class”.

islower()#

Return True if the string is a lowercase string, False otherwise.

A string is lowercase if all cased characters in the string are lowercase and there is at least one cased character in the string.

isnumeric()#

Return True if the string is a numeric string, False otherwise.

A string is numeric if all characters in the string are numeric and there is at least one character in the string.

isprintable()#

Return True if the string is printable, False otherwise.

A string is printable if all of its characters are considered printable in repr() or if it is empty.

isspace()#

Return True if the string is a whitespace string, False otherwise.

A string is whitespace if all characters in the string are whitespace and there is at least one character in the string.

istitle()#

Return True if the string is a title-cased string, False otherwise.

In a title-cased string, upper- and title-case characters may only follow uncased characters and lowercase characters only cased ones.

isupper()#

Return True if the string is an uppercase string, False otherwise.

A string is uppercase if all cased characters in the string are uppercase and there is at least one cased character in the string.

join(iterable, /)#

Concatenate any number of strings.

The string whose method is called is inserted in between each given string. The result is returned as a new string.

Example: ‘.’.join([‘ab’, ‘pq’, ‘rs’]) -> ‘ab.pq.rs’

ljust(width, fillchar=' ', /)#

Return a left-justified string of length width.

Padding is done using the specified fill character (default is a space).

lower()#

Return a copy of the string converted to lowercase.

lstrip(chars=None, /)#

Return a copy of the string with leading whitespace removed.

If chars is given and not None, remove characters in chars instead.

static maketrans()#

Return a translation table usable for str.translate().

If there is only one argument, it must be a dictionary mapping Unicode ordinals (integers) or characters to Unicode ordinals, strings or None. Character keys will be then converted to ordinals. If there are two arguments, they must be strings of equal length, and in the resulting dictionary, each character in x will be mapped to the character at the same position in y. If there is a third argument, it must be a string, whose characters will be mapped to None in the result.

partition(sep, /)#

Partition the string into three parts using the given separator.

This will search for the separator in the string. If the separator is found, returns a 3-tuple containing the part before the separator, the separator itself, and the part after it.

If the separator is not found, returns a 3-tuple containing the original string and two empty strings.

removeprefix(prefix, /)#

Return a str with the given prefix string removed if present.

If the string starts with the prefix string, return string[len(prefix):]. Otherwise, return a copy of the original string.

removesuffix(suffix, /)#

Return a str with the given suffix string removed if present.

If the string ends with the suffix string and that suffix is not empty, return string[:-len(suffix)]. Otherwise, return a copy of the original string.

replace(old, new, count=-1, /)#

Return a copy with all occurrences of substring old replaced by new.

count

Maximum number of occurrences to replace. -1 (the default value) means replace all occurrences.

If the optional argument count is given, only the first count occurrences are replaced.

rfind(sub[, start[, end]]) → int#

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

rindex(sub[, start[, end]]) → int#

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Raises ValueError when the substring is not found.

rjust(width, fillchar=' ', /)#

Return a right-justified string of length width.

Padding is done using the specified fill character (default is a space).

rpartition(sep, /)#

Partition the string into three parts using the given separator.

This will search for the separator in the string, starting at the end. If the separator is found, returns a 3-tuple containing the part before the separator, the separator itself, and the part after it.

If the separator is not found, returns a 3-tuple containing two empty strings and the original string.

rsplit(sep=None, maxsplit=-1)#

Return a list of the substrings in the string, using sep as the separator string.

sep

The separator used to split the string.

When set to None (the default value), will split on any whitespace character (including n r t f and spaces) and will discard empty strings from the result.

maxsplit

Maximum number of splits. -1 (the default value) means no limit.

Splitting starts at the end of the string and works to the front.

rstrip(chars=None, /)#

Return a copy of the string with trailing whitespace removed.

If chars is given and not None, remove characters in chars instead.

split(sep=None, maxsplit=-1)#

Return a list of the substrings in the string, using sep as the separator string.

sep

The separator used to split the string.

When set to None (the default value), will split on any whitespace character (including n r t f and spaces) and will discard empty strings from the result.

maxsplit

Maximum number of splits. -1 (the default value) means no limit.

Splitting starts at the front of the string and works to the end.

Note, str.split() is mainly useful for data that has been intentionally delimited. With natural text that includes punctuation, consider using the regular expression module.

splitlines(keepends=False)#

Return a list of the lines in the string, breaking at line boundaries.

Line breaks are not included in the resulting list unless keepends is given and true.

startswith(prefix[, start[, end]]) → bool#

Return True if S starts with the specified prefix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. prefix can also be a tuple of strings to try.

strip(chars=None, /)#

Return a copy of the string with leading and trailing whitespace removed.

If chars is given and not None, remove characters in chars instead.

swapcase()#

Convert uppercase characters to lowercase and lowercase characters to uppercase.

title()#

Return a version of the string where each word is titlecased.

More specifically, words start with uppercased characters and all remaining cased characters have lower case.

translate(table, /)#

Replace each character in the string using the given translation table.

table

Translation table, which must be a mapping of Unicode ordinals to Unicode ordinals, strings, or None.

The table must implement lookup/indexing via __getitem__, for instance a dictionary or list. If this operation raises LookupError, the character is left untouched. Characters mapped to None are deleted.

upper()#

Return a copy of the string converted to uppercase.

zfill(width, /)#

Pad a numeric string with zeros on the left, to fill a field of the given width.

The string is never truncated.

pulse.problem.CONTROLLED_VOLUME_SCALE = 1000000.0#

The rows of a controlled cavity measure the volume in mL and the pressure in kPa, so that both modes have residuals of order one for a heart and Newton’s tolerance means the same thing whichever mode is active.

class pulse.problem.Cavity(marker, volume=None, control=None)[source]#

A chamber whose pressure is an unknown of the problem, and its constraint.

Exactly one of these constrains it:

  • volume: the volume is held at that value, with the pressure as its Lagrange multiplier. It is usually a Constant you set each step.

  • control: a CavityControl, whose constraint (volume, pressure, or pressure affine in the volume) can be switched at run time.

  • neither, when the marker is coupled to a chamber of a 0D circulation model (see pulse.circulation): the problem then points volume at that chamber’s volume state, so the volume is itself an unknown and the constraint couples the two.

Parameters:
control: CavityControl | None#

Alias for field number 2

count(value, /)#

Return number of occurrences of value.

index(value, start=0, stop=9223372036854775807, /)#

Return first index of value.

Raises ValueError if the value is not present.

marker: str#

Alias for field number 0

volume: Constant | Function | Expr | None#

Alias for field number 1

class pulse.problem.CavityControl(mesh)[source]#

Which constraint a cavity’s pressure unknown satisfies, set at run time.

A controlled cavity has one pressure unknown p, as a prescribed-volume cavity does, but its equation is chosen by the constants held here:

  • volume mode (mode == 1): V(u) = V_target;

  • pressure mode (mode == 0): p = A + B V(u). B = 0 prescribes the pressure, and B != 0 ties the pressure to the volume inside Newton, as a Windkessel does during ejection.

Every one of these is a Constant read when the form is assembled, so switching between them does not rebuild the problem. A new control starts in pressure mode at zero pressure. Values are in SI units: V_target in m^3, A in Pa and B in Pa/m^3 – which only means what it says when the problem’s own V(u) is in cubic metres, so a controlled cavity requires parameters["mesh_unit"] == "m"; StaticProblem refuses one otherwise.

Do not also put a Neumann pressure on the cavity’s marker: the load on the wall comes from p.

Parameters:

mesh (Mesh)

set_affine_pressure(A, B)[source]#

Make the cavity pressure A + B V(u), with A in Pa and B in Pa/m^3.

Parameters:

  • A (float)

  • B (float)

Return type:

None

set_pressure(P)[source]#

Hold the cavity pressure at P (Pa).

Parameters:

P (float)

Return type:

None

set_volume(V)[source]#

Hold the cavity volume at V (m^3).

Parameters:

V (float)

Return type:

None

class pulse.problem.DynamicProblem(model, geometry, parameters=<factory>, bcs=<factory>, cavities=<factory>, circulation=None, chambers=<factory>, circulation_missing=<factory>, NonlinearProblem=<class 'dolfinx.fem.petsc.NonlinearProblem'>, Function=<class 'dolfinx.fem.function.Function'>, monitor=<factory>)[source]#

Second-order elastodynamics, via the generalized-\(\alpha\) method.

Every term of the residual is evaluated at one of three points:

  • at the \(\alpha_f\) point (interpolate(u_old, u, alpha_f), built in R): the material, compressibility and viscous stress, the Robin and Neumann loads, and the body force – the genuine second-order dynamics, for which alpha_f-interpolation is the consistent choice.

  • at the \(\alpha_m\) point: the inertia term.

  • at the end of the step, i.e. the true current self.u: the cavity-volume constraint, \(J - 1\), and the circuit. These are algebraic constraints (Lagrange multipliers), not part of the differential dynamics an alpha-blend applies to – enforcing them against the alpha_f-interpolated configuration instead would leave self.u’s actual volume/incompressibility/circuit coupling unconstrained. The active stress joins this group precisely when evaluate_at_end_of_step is set on the active model: such a model’s stress depends on state advanced over the step (a stretch rate, or ODE states), and the alpha_f point would feed it a blended stretch and a rate scaled by \(1 - \alpha_f\) rather than the ones it actually advanced with.

Parameters:
class Function(*args, **kw)#

A finite element function.

A finite element function is represented by a function space (domain, element and dofmap) and a vector holding the degrees-of-freedom.

Parameters:

  • V (FunctionSpace)

  • x (la.Vector[Scalar] | None)

  • name (str | None)

  • dtype (npt.DTypeLike | None)

property T#

Transpose a rank-2 tensor expression.

For more general transpose operations of higher order tensor expressions, use indexing and Tensor.

collapse()#

Create a collapsed version of this Function.

Return type:

Function[Scalar]

copy()#

Create a copy of the Function.

The function space is shared and the degree-of-freedom vector is copied.

Returns:

A new Function with a copy of the degree-of-freedom vector.

Return type:

Function[Scalar]

count()#

Get count.

property dtype: DTypeLike#

Function value dtype.

dx(*ii)#

Return the partial derivative with respect to spatial variable number ii.

eval(x, cells, u=None, tol=1e-06, maxit=15)#

Evaluate Function at points x.

Parameters:

  • x (ArrayLike) – Points with shape (num_points, 3)

  • cells (ndarray[tuple[Any, ...], dtype[int32]]) – Array with cell indices, with shape (num_points,), where cell[i] is the index of the cell containing point x[i]. If the cell index is negative the point is ignored.

  • u (None | ndarray[tuple[Any, ...], dtype[Scalar]]) – Array to put evaluated data in.

  • tol (float) – Tolerance for convergence in Newton method for nonaffine pullbacks.

  • maxit (int) – Maximum number of Newton iterations for nonaffine pullbacks.

Return type:

ndarray[tuple[Any, …], dtype[Scalar]]

evaluate(x, mapping, component, index_values, derivatives=())#

Get self from mapping and return the component asked for.

property function_space: FunctionSpace#

FunctionSpace that the Function is defined on.

interpolate(u0, cells0=None, cells1=None)#

Interpolate an expression.

Parameters:

  • u0 (Callable | Expression[Scalar] | Function[Scalar]) – Callable function, Expression or Function to interpolate.

  • cells0 (ndarray[tuple[Any, ...], dtype[int32]] | None) – Cells in mesh associated with u0 to interpolate over. If None then all cells are interpolated over.

  • cells1 (ndarray[tuple[Any, ...], dtype[int32]] | None) – Cells in the mesh associated with self to interpolate over. If None, then taken to be the same cells as cells0. If cells1 is not None it must have the same length as cells0.

Return type:

None

interpolate_nonmatching(u0, cells, interpolation_data, tol=1e-06, maxit=15)#

Interpolate a Function on a non-matching mesh.

Parameters:

  • u0 (Function[Scalar]) – Function to interpolate.

  • cells (ndarray[tuple[Any, ...], dtype[int32]]) – The cells to interpolate over. If None then all cells are interpolated over.

  • interpolation_data (PointOwnershipData) – Data needed to interpolate functions defined on other meshes. Created by dolfinx.fem.create_interpolation_data().

  • tol (float) – Tolerance for convergence in Newton method for nonaffine pullbacks. Ignored if mesh geometry is affine.

  • maxit (int) – Maximum number of iterations for nonaffine pullback. Ignored if mesh geometry is affine.

Return type:

None

is_cellwise_constant()#

Return whether this expression is spatially constant over each cell.

property name: str#

Name of the Function.

split()#

Extract (any) sub-functions.

A sub-function can be extracted from a discrete function that is in a mixed, vector, or tensor FunctionSpace. The sub-function resides in the subspace of the mixed space.

Returns:

First level of subspaces of the function space.

Return type:

tuple[Function[Scalar], …]

sub(i)#

Return a sub-function (a view into the Function).

Sub-functions are indexed i = 0, ..., N-1, where N is the number of sub-spaces.

Parameters:

i (int) – Index of the sub-function to extract.

Returns:

A view into the parent Function.

Return type:

Function[Scalar]

Note

If the sub-Function is re-used, for performance reasons the returned Function should be stored by the caller to avoid repeated re-computation of the subspace.

static ufl_disable_profiling()#

Turn off the object counting mechanism. Return object init and del counts.

ufl_domain()#

Return the single unique domain this expression is defined on, or throw an error.

ufl_domains()#

Get the UFL domains.

ufl_element()#

Shortcut to get the finite element of the function space of this coefficient.

static ufl_enable_profiling()#

Turn on the object counting mechanism and reset counts to zero.

ufl_function_space()#

Get the function space of this coefficient.

property ufl_shape#

Return the associated UFL shape.

property x: Vector[Scalar]#

Vector holding the degrees-of-freedom.

class NonlinearProblem(F, u, *, petsc_options_prefix, bcs=None, J=None, P=None, kind=None, petsc_options=None, form_compiler_options=None, jit_options=None, entity_maps=None)#

High-level class for solving nonlinear problems with PETSc SNES.

Solves problems of the form \(F_i(u, v) = 0, i=0,\ldots,N\ \forall v \in V\) where \(u=(u_0,\ldots,u_N), v=(v_0,\ldots,v_N)\) using PETSc SNES as the non-linear solver.

Note

The deprecated version of this class for use with dolfinx.nls.petsc.NewtonSolver has been renamed dolfinx.fem.petsc.NewtonSolverNonlinearProblem.

Note

This high-level class automatically handles PETSc memory management. The user does not need to manually call .destroy() on returned PETSc objects.

Parameters:

  • F (ufl.form.Form | Sequence[ufl.form.Form])

  • u (_Function | Sequence[_Function])

  • petsc_options_prefix (str)

  • bcs (Sequence[DirichletBC] | None)

  • J (ufl.form.Form | Sequence[Sequence[ufl.form.Form]] | None)

  • P (ufl.form.Form | Sequence[Sequence[ufl.form.Form]] | None)

  • kind (str | Sequence[Sequence[str]] | None)

  • petsc_options (dict | None)

  • form_compiler_options (dict | None)

  • jit_options (dict | None)

  • entity_maps (Sequence[_EntityMap] | None)

property A: Mat#

Jacobian matrix.

property F: Form | Sequence[Form]#

The compiled residual.

property J: Form | Sequence[Sequence[Form]]#

The compiled Jacobian.

property P_mat: Mat | None#

Preconditioner matrix.

property b: Vec#

Residual vector.

property preconditioner: Form | Sequence[Sequence[Form]] | None#

The compiled preconditioner.

solve()#

Solve the problem.

This method updates the solution u function(s) stored in the problem instance.

Note

The user is responsible for asserting convergence of the SNES solver e.g. assert problem.solver.getConvergedReason() > 0. Alternatively, pass "snes_error_if_not_converged": True and "ksp_error_if_not_converged" : True in petsc_options to raise a PETScError on failure.

Returns:

The solution function(s).

Return type:

Function | Sequence[Function]

property solver: SNES#

The SNES solver.

property u: Function | Sequence[Function]#

Solution function(s).

Note

The function(s) do not share memory with the solution vector x.

property x: Vec#

Solution vector.

Note

The vector does not share memory with the solution function(s) u.

a(u, u_old, v_old, a_old, dt=None)[source]#

Acceleration computed using the generalized \(alpha\)-method

\[a_{i+1} = \frac{u_{i+1} - (u_i + \Delta t v_i + (0.5 - \beta) \Delta t^2 a_i)}{\beta \Delta t^2}\]
Parameters:

  • u (T) – Current displacement

  • u_old (T) – Previous displacement

  • v_old (T) – Previous velocity

  • a_old (T) – Previous acceleration

  • dt (optional) – Time step; see v().

Returns:

The current acceleration

Return type:

T

property incompressibility_index: int#

Index of the incompressibility row in the block system.

The order is (u, cavity pressures, rigid body, p, circulation states), so p is only the last row when there is no circulation model. Naming the row rather than counting back from the end keeps the constraint where it belongs when unknowns are added after it.

load_restart_metadata(data)#

Restore what restart_metadata() returned, after the Functions.

Sets the circuit’s step count and the stencil it selects, backward Euler included, so a snapshot taken before the step that switched to BDF2 switches it back. Raises ValueError if data and this problem disagree on whether there is a circulation.

Parameters:

data (Mapping[str, Any])

Return type:

None

reset_states()#

Reset states to old values

restart_functions()[source]#

StaticProblem.restart_functions, then the velocity and acceleration history.

Return type:

list[tuple[str, Function]]

restart_metadata()#

The restart state that is not a Function, JSON-able.

With a circulation, that is the number of converged steps the circuit has taken, which selects its stencil (BDF2 takes its first step as backward Euler). Without one, nothing.

Return type:

dict[str, Any]

solve(update_old_states=True, raise_on_failure=None)#

Solve nonlinear problem with Newton solver

Parameters:

  • update_old_states (bool, optional) – Whether to update old states before solving, by default True

  • raise_on_failure (bool | None, optional) – Whether a Newton solve that does not converge should raise rather than return False. None, the default, defers to parameters["raise_on_failure"], which is itself False.

Returns:

True if converged, False otherwise. When this is False the state is whatever the failed solve left behind and should not be used; reset_states() puts it back.

Return type:

bool

Notes

This overrides snes_error_if_not_converged (and the corresponding KSP option) on every call, so setting them through petsc_options has no effect – raise_on_failure is the single control.

Non-convergence is logged at warning level whichever way it is reported, because a caller that ignores the return value would otherwise carry on against a failed solve in silence.

update_fields()[source]#

Update old values of displacement, velocity and acceleration

Return type:

None

update_old_states()#

Update old states to current values

v(a, v_old, a_old, dt=None)[source]#

Velocity computed using the generalized \(alpha\)-method

\[v_{i+1} = v_i + (1-\gamma) \Delta t a_i + \gamma \Delta t a_{i+1}\]
Parameters:

  • a (T) – Current acceleration

  • v_old (T) – Previous velocity

  • a_old (T) – Previous acceleration

  • dt (optional) – Time step; defaults to parameters["dt"] in base units (a UFL expression when dt wraps a Constant). Pass a float for numpy input.

Returns:

The current velocity

Return type:

T

class pulse.problem.Geometry(*args, **kwargs)[source]#

Protocol for geometry objects used in mechanics problems.

class pulse.problem.StaticProblem(model: pulse.cardiac_model.CardiacModel, geometry: pulse.problem.Geometry, parameters: dict[str, typing.Any]=<factory>, bcs: pulse.boundary_conditions.BoundaryConditions = <factory>, cavities: list[pulse.problem.Cavity] = <factory>, circulation: pulse.circulation.CirculationModel | None = None, chambers: list[pulse.circulation.ChamberCoupling] = <factory>, circulation_missing: dict[str, typing.Any]=<factory>, NonlinearProblem: Type[dolfinx.fem.petsc.NonlinearProblem] = <class 'dolfinx.fem.petsc.NonlinearProblem'>, Function: Type[dolfinx.fem.function.Function] = <class 'dolfinx.fem.function.Function'>, monitor: pulse.telemetry.BaseMonitor = <factory>)[source]#
Parameters:
class Function(*args, **kw)#

A finite element function.

A finite element function is represented by a function space (domain, element and dofmap) and a vector holding the degrees-of-freedom.

Parameters:

  • V (FunctionSpace)

  • x (la.Vector[Scalar] | None)

  • name (str | None)

  • dtype (npt.DTypeLike | None)

property T#

Transpose a rank-2 tensor expression.

For more general transpose operations of higher order tensor expressions, use indexing and Tensor.

collapse()#

Create a collapsed version of this Function.

Return type:

Function[Scalar]

copy()#

Create a copy of the Function.

The function space is shared and the degree-of-freedom vector is copied.

Returns:

A new Function with a copy of the degree-of-freedom vector.

Return type:

Function[Scalar]

count()#

Get count.

property dtype: DTypeLike#

Function value dtype.

dx(*ii)#

Return the partial derivative with respect to spatial variable number ii.

eval(x, cells, u=None, tol=1e-06, maxit=15)#

Evaluate Function at points x.

Parameters:

  • x (ArrayLike) – Points with shape (num_points, 3)

  • cells (ndarray[tuple[Any, ...], dtype[int32]]) – Array with cell indices, with shape (num_points,), where cell[i] is the index of the cell containing point x[i]. If the cell index is negative the point is ignored.

  • u (None | ndarray[tuple[Any, ...], dtype[Scalar]]) – Array to put evaluated data in.

  • tol (float) – Tolerance for convergence in Newton method for nonaffine pullbacks.

  • maxit (int) – Maximum number of Newton iterations for nonaffine pullbacks.

Return type:

ndarray[tuple[Any, …], dtype[Scalar]]

evaluate(x, mapping, component, index_values, derivatives=())#

Get self from mapping and return the component asked for.

property function_space: FunctionSpace#

FunctionSpace that the Function is defined on.

interpolate(u0, cells0=None, cells1=None)#

Interpolate an expression.

Parameters:

  • u0 (Callable | Expression[Scalar] | Function[Scalar]) – Callable function, Expression or Function to interpolate.

  • cells0 (ndarray[tuple[Any, ...], dtype[int32]] | None) – Cells in mesh associated with u0 to interpolate over. If None then all cells are interpolated over.

  • cells1 (ndarray[tuple[Any, ...], dtype[int32]] | None) – Cells in the mesh associated with self to interpolate over. If None, then taken to be the same cells as cells0. If cells1 is not None it must have the same length as cells0.

Return type:

None

interpolate_nonmatching(u0, cells, interpolation_data, tol=1e-06, maxit=15)#

Interpolate a Function on a non-matching mesh.

Parameters:

  • u0 (Function[Scalar]) – Function to interpolate.

  • cells (ndarray[tuple[Any, ...], dtype[int32]]) – The cells to interpolate over. If None then all cells are interpolated over.

  • interpolation_data (PointOwnershipData) – Data needed to interpolate functions defined on other meshes. Created by dolfinx.fem.create_interpolation_data().

  • tol (float) – Tolerance for convergence in Newton method for nonaffine pullbacks. Ignored if mesh geometry is affine.

  • maxit (int) – Maximum number of iterations for nonaffine pullback. Ignored if mesh geometry is affine.

Return type:

None

is_cellwise_constant()#

Return whether this expression is spatially constant over each cell.

property name: str#

Name of the Function.

split()#

Extract (any) sub-functions.

A sub-function can be extracted from a discrete function that is in a mixed, vector, or tensor FunctionSpace. The sub-function resides in the subspace of the mixed space.

Returns:

First level of subspaces of the function space.

Return type:

tuple[Function[Scalar], …]

sub(i)#

Return a sub-function (a view into the Function).

Sub-functions are indexed i = 0, ..., N-1, where N is the number of sub-spaces.

Parameters:

i (int) – Index of the sub-function to extract.

Returns:

A view into the parent Function.

Return type:

Function[Scalar]

Note

If the sub-Function is re-used, for performance reasons the returned Function should be stored by the caller to avoid repeated re-computation of the subspace.

static ufl_disable_profiling()#

Turn off the object counting mechanism. Return object init and del counts.

ufl_domain()#

Return the single unique domain this expression is defined on, or throw an error.

ufl_domains()#

Get the UFL domains.

ufl_element()#

Shortcut to get the finite element of the function space of this coefficient.

static ufl_enable_profiling()#

Turn on the object counting mechanism and reset counts to zero.

ufl_function_space()#

Get the function space of this coefficient.

property ufl_shape#

Return the associated UFL shape.

property x: Vector[Scalar]#

Vector holding the degrees-of-freedom.

class NonlinearProblem(F, u, *, petsc_options_prefix, bcs=None, J=None, P=None, kind=None, petsc_options=None, form_compiler_options=None, jit_options=None, entity_maps=None)#

High-level class for solving nonlinear problems with PETSc SNES.

Solves problems of the form \(F_i(u, v) = 0, i=0,\ldots,N\ \forall v \in V\) where \(u=(u_0,\ldots,u_N), v=(v_0,\ldots,v_N)\) using PETSc SNES as the non-linear solver.

Note

The deprecated version of this class for use with dolfinx.nls.petsc.NewtonSolver has been renamed dolfinx.fem.petsc.NewtonSolverNonlinearProblem.

Note

This high-level class automatically handles PETSc memory management. The user does not need to manually call .destroy() on returned PETSc objects.

Parameters:

  • F (ufl.form.Form | Sequence[ufl.form.Form])

  • u (_Function | Sequence[_Function])

  • petsc_options_prefix (str)

  • bcs (Sequence[DirichletBC] | None)

  • J (ufl.form.Form | Sequence[Sequence[ufl.form.Form]] | None)

  • P (ufl.form.Form | Sequence[Sequence[ufl.form.Form]] | None)

  • kind (str | Sequence[Sequence[str]] | None)

  • petsc_options (dict | None)

  • form_compiler_options (dict | None)

  • jit_options (dict | None)

  • entity_maps (Sequence[_EntityMap] | None)

property A: Mat#

Jacobian matrix.

property F: Form | Sequence[Form]#

The compiled residual.

property J: Form | Sequence[Sequence[Form]]#

The compiled Jacobian.

property P_mat: Mat | None#

Preconditioner matrix.

property b: Vec#

Residual vector.

property preconditioner: Form | Sequence[Sequence[Form]] | None#

The compiled preconditioner.

solve()#

Solve the problem.

This method updates the solution u function(s) stored in the problem instance.

Note

The user is responsible for asserting convergence of the SNES solver e.g. assert problem.solver.getConvergedReason() > 0. Alternatively, pass "snes_error_if_not_converged": True and "ksp_error_if_not_converged" : True in petsc_options to raise a PETScError on failure.

Returns:

The solution function(s).

Return type:

Function | Sequence[Function]

property solver: SNES#

The SNES solver.

property u: Function | Sequence[Function]#

Solution function(s).

Note

The function(s) do not share memory with the solution vector x.

property x: Vec#

Solution vector.

Note

The vector does not share memory with the solution function(s) u.

property incompressibility_index: int#

Index of the incompressibility row in the block system.

The order is (u, cavity pressures, rigid body, p, circulation states), so p is only the last row when there is no circulation model. Naming the row rather than counting back from the end keeps the constraint where it belongs when unknowns are added after it.

load_restart_metadata(data)[source]#

Restore what restart_metadata() returned, after the Functions.

Sets the circuit’s step count and the stencil it selects, backward Euler included, so a snapshot taken before the step that switched to BDF2 switches it back. Raises ValueError if data and this problem disagree on whether there is a circulation.

Parameters:

data (Mapping[str, Any])

Return type:

None

reset_states()[source]#

Reset states to old values

restart_functions()[source]#

Every Function a restart needs, under mechanics_* names, in a fixed order.

These are the problem’s own Functions, not copies: write a checkpoint from them, or restore one by writing into them in place, then call load_restart_metadata(). The order is u; p if incompressible; each cavity’s pressure, in cavities order; the rigid-body multiplier r; then each circulation state. Each is followed by its _old copy, and a circulation state also by its _prev (BDF2’s second level). DynamicProblem appends v_old and a_old. With no cavity, rigid body or circulation, this is the list the CLI wrote before it called this method, so its earlier checkpoints still restore.

A solve overwrites the _old copies before it reads them, unless called with update_old_states=False.

Return type:

list[tuple[str, Function]]

restart_metadata()[source]#

The restart state that is not a Function, JSON-able.

With a circulation, that is the number of converged steps the circuit has taken, which selects its stencil (BDF2 takes its first step as backward Euler). Without one, nothing.

Return type:

dict[str, Any]

solve(update_old_states=True, raise_on_failure=None)[source]#

Solve nonlinear problem with Newton solver

Parameters:

  • update_old_states (bool, optional) – Whether to update old states before solving, by default True

  • raise_on_failure (bool | None, optional) – Whether a Newton solve that does not converge should raise rather than return False. None, the default, defers to parameters["raise_on_failure"], which is itself False.

Returns:

True if converged, False otherwise. When this is False the state is whatever the failed solve left behind and should not be used; reset_states() puts it back.

Return type:

bool

Notes

This overrides snes_error_if_not_converged (and the corresponding KSP option) on every call, so setting them through petsc_options has no effect – raise_on_failure is the single control.

Non-convergence is logged at warning level whichever way it is reported, because a caller that ignores the return value would otherwise carry on against a failed solve in silence.

update_fields()[source]#

Shift the circuit’s second history level, for BDF2.

Called only after a converged solve. update_old_states runs before Newton and again on every retry, so shifting there would consume a history level per attempt instead of per step.

circulation_states_old still holds the level update_old_states is about to overwrite, so copying it into _prev here lines up (y_{n+1}, y_n, y_{n-1}) for the next step.

update_old_states()[source]#

Update old states to current values

pulse.problem.interpolate(x0, x1, alpha)[source]#

Interpolate between \(x_0\) and \(x_1\) to find math:`x_{1-alpha}

Parameters:

  • x0 (T) – First point

  • x1 (T) – Second point

  • alpha (float) – Amount of interpolate

Returns:

math:`x_{1-alpha}

Return type:

T

pulse.problem.kPa = 1000.0#

1 kPa in pascals – the pressure analogue of mL (1 mL in cubic metres) above. There is no shared kPa constant to import for this, unlike mL, so it is defined right here, next to the one row that uses it.

pulse.problem.volume_scale(mesh_unit)[source]#

The factor turning a volume in ``mesh_unit``^3 into mL.

Every cavity’s volume row is multiplied by it, so the row’s residual is in mL whatever the mesh unit. In m^3 a heart’s volume changes by less than snes_atol (1e-6, i.e. 1 mL) in a typical step, so Newton could stop before the wall had moved at all.

Parameters:

mesh_unit (str)

Return type:

float

boundary conditions#

This module defines boundary conditions.

Boundary conditions are used to specify the behavior of the solution on the boundary of the domain. The boundary conditions can be Dirichlet, Neumann, or Robin boundary conditions.

Dirichlet boundary conditions are used to specify the solution on the boundary of the domain. Neumann boundary conditions are used to specify the traction on the boundary of the domain. Robin boundary conditions are used to specify a Robin type boundary condition on the boundary of the domain.

The boundary conditions are collected in a BoundaryConditions object.

class pulse.boundary_conditions.BoundaryConditions(neumann, dirichlet, robin, body_force)[source]#
Parameters:
body_force: Sequence[float | Constant | Function]#

Alias for field number 3

count(value, /)#

Return number of occurrences of value.

dirichlet: Sequence[Callable[[FunctionSpace], Sequence[DirichletBC]]]#

Alias for field number 1

index(value, start=0, stop=9223372036854775807, /)#

Return first index of value.

Raises ValueError if the value is not present.

neumann: Sequence[NeumannBC]#

Alias for field number 0

robin: Sequence[RobinBC]#

Alias for field number 2

class pulse.boundary_conditions.NeumannBC(traction: pulse.units.Variable, marker: int)[source]#
Parameters:
class pulse.boundary_conditions.RobinBC(value, marker, damping=False, perpendicular=False, normal=None)[source]#

A spring (damping=False) or dashpot (damping=True) on the facets marker.

It acts along the surface normal, or, with perpendicular=True, in the tangent plane. normal chooses the current or the reference normal; see RobinNormal. Left unset, it is current for a spring and reference for a dashpot, as before the option existed. The spring is at rest in the reference configuration.

Parameters:

  • value (Variable)

  • marker (int)

  • damping (bool)

  • perpendicular (bool)

  • normal (RobinNormal | None)

projection(N, F, mesh_is_reference=True)[source]#

Return the projector this BC acts with and the area ratio from the mesh to it.

N is the unit normal of the mesh. F is the deformation gradient from the mesh to the other configuration: the current one in a forward problem, and the reference one in the inverse (prestress) problem, where mesh_is_reference=False. Integrate the traction against the mesh’s ds times the returned ratio.

Parameters:

  • N (Expr)

  • F (Expr)

  • mesh_is_reference (bool)

class pulse.boundary_conditions.RobinNormal(*values)[source]#

Which surface normal a RobinBC spring or dashpot acts along.

Both measure the spring’s extension by projecting the displacement (or velocity) onto a normal, so neither is the true distance to the surface the spring is anchored to.

current

The normal \(\mathbf{n}\) and area \(da\) of the current configuration, pushed forward with Nanson’s formula, \(\mathbf{t} = k (\mathbf{u} \cdot \mathbf{n}) \mathbf{n}\). It derives from no energy. It follows the surface as it deforms, so it also resists displacement that turns normal as the wall rotates.

reference

The normal \(\mathbf{N}\) and area \(dA\) of the reference configuration, \(\mathbf{t} = k (\mathbf{u} \cdot \mathbf{N}) \mathbf{N}\) (Pfaller et al. 2019, eqs. 4-5). It is the derivative of the energy \(\frac{1}{2} \int k (\mathbf{u} \cdot \mathbf{N})^2 \, dA\), so its Jacobian is symmetric. It assumes small rotations of the surface. Under large deformation it lets a free base flare: in the fixed-point unloader demo, the same springs nearly double the inflation, and the unloaded base comes out wider than the loaded one.

Leaving RobinBC.normal unset keeps what pulse did before the option existed, which the demos and templates were tuned with: springs act along current, dashpots along reference.

On a curved surface both forms register pure tangential sliding as a change in the normal gap, of second order in the rotation. So a stiff epicardial spring holds back ventricular twist. reference does so more than current.

capitalize()#

Return a capitalized version of the string.

More specifically, make the first character have upper case and the rest lower case.

casefold()#

Return a version of the string suitable for caseless comparisons.

center(width, fillchar=' ', /)#

Return a centered string of length width.

Padding is done using the specified fill character (default is a space).

count(sub[, start[, end]]) → int#

Return the number of non-overlapping occurrences of substring sub in string S[start:end]. Optional arguments start and end are interpreted as in slice notation.

encode(encoding='utf-8', errors='strict')#

Encode the string using the codec registered for encoding.

encoding

The encoding in which to encode the string.

errors

The error handling scheme to use for encoding errors. The default is ‘strict’ meaning that encoding errors raise a UnicodeEncodeError. Other possible values are ‘ignore’, ‘replace’ and ‘xmlcharrefreplace’ as well as any other name registered with codecs.register_error that can handle UnicodeEncodeErrors.

endswith(suffix[, start[, end]]) → bool#

Return True if S ends with the specified suffix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. suffix can also be a tuple of strings to try.

expandtabs(tabsize=8)#

Return a copy where all tab characters are expanded using spaces.

If tabsize is not given, a tab size of 8 characters is assumed.

find(sub[, start[, end]]) → int#

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

format(*args, **kwargs) → str#

Return a formatted version of S, using substitutions from args and kwargs. The substitutions are identified by braces (‘{’ and ‘}’).

format_map(mapping) → str#

Return a formatted version of S, using substitutions from mapping. The substitutions are identified by braces (‘{’ and ‘}’).

index(sub[, start[, end]]) → int#

Return the lowest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Raises ValueError when the substring is not found.

isalnum()#

Return True if the string is an alpha-numeric string, False otherwise.

A string is alpha-numeric if all characters in the string are alpha-numeric and there is at least one character in the string.

isalpha()#

Return True if the string is an alphabetic string, False otherwise.

A string is alphabetic if all characters in the string are alphabetic and there is at least one character in the string.

isascii()#

Return True if all characters in the string are ASCII, False otherwise.

ASCII characters have code points in the range U+0000-U+007F. Empty string is ASCII too.

isdecimal()#

Return True if the string is a decimal string, False otherwise.

A string is a decimal string if all characters in the string are decimal and there is at least one character in the string.

isdigit()#

Return True if the string is a digit string, False otherwise.

A string is a digit string if all characters in the string are digits and there is at least one character in the string.

isidentifier()#

Return True if the string is a valid Python identifier, False otherwise.

Call keyword.iskeyword(s) to test whether string s is a reserved identifier, such as “def” or “class”.

islower()#

Return True if the string is a lowercase string, False otherwise.

A string is lowercase if all cased characters in the string are lowercase and there is at least one cased character in the string.

isnumeric()#

Return True if the string is a numeric string, False otherwise.

A string is numeric if all characters in the string are numeric and there is at least one character in the string.

isprintable()#

Return True if the string is printable, False otherwise.

A string is printable if all of its characters are considered printable in repr() or if it is empty.

isspace()#

Return True if the string is a whitespace string, False otherwise.

A string is whitespace if all characters in the string are whitespace and there is at least one character in the string.

istitle()#

Return True if the string is a title-cased string, False otherwise.

In a title-cased string, upper- and title-case characters may only follow uncased characters and lowercase characters only cased ones.

isupper()#

Return True if the string is an uppercase string, False otherwise.

A string is uppercase if all cased characters in the string are uppercase and there is at least one cased character in the string.

join(iterable, /)#

Concatenate any number of strings.

The string whose method is called is inserted in between each given string. The result is returned as a new string.

Example: ‘.’.join([‘ab’, ‘pq’, ‘rs’]) -> ‘ab.pq.rs’

ljust(width, fillchar=' ', /)#

Return a left-justified string of length width.

Padding is done using the specified fill character (default is a space).

lower()#

Return a copy of the string converted to lowercase.

lstrip(chars=None, /)#

Return a copy of the string with leading whitespace removed.

If chars is given and not None, remove characters in chars instead.

static maketrans()#

Return a translation table usable for str.translate().

If there is only one argument, it must be a dictionary mapping Unicode ordinals (integers) or characters to Unicode ordinals, strings or None. Character keys will be then converted to ordinals. If there are two arguments, they must be strings of equal length, and in the resulting dictionary, each character in x will be mapped to the character at the same position in y. If there is a third argument, it must be a string, whose characters will be mapped to None in the result.

partition(sep, /)#

Partition the string into three parts using the given separator.

This will search for the separator in the string. If the separator is found, returns a 3-tuple containing the part before the separator, the separator itself, and the part after it.

If the separator is not found, returns a 3-tuple containing the original string and two empty strings.

removeprefix(prefix, /)#

Return a str with the given prefix string removed if present.

If the string starts with the prefix string, return string[len(prefix):]. Otherwise, return a copy of the original string.

removesuffix(suffix, /)#

Return a str with the given suffix string removed if present.

If the string ends with the suffix string and that suffix is not empty, return string[:-len(suffix)]. Otherwise, return a copy of the original string.

replace(old, new, count=-1, /)#

Return a copy with all occurrences of substring old replaced by new.

count

Maximum number of occurrences to replace. -1 (the default value) means replace all occurrences.

If the optional argument count is given, only the first count occurrences are replaced.

rfind(sub[, start[, end]]) → int#

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Return -1 on failure.

rindex(sub[, start[, end]]) → int#

Return the highest index in S where substring sub is found, such that sub is contained within S[start:end]. Optional arguments start and end are interpreted as in slice notation.

Raises ValueError when the substring is not found.

rjust(width, fillchar=' ', /)#

Return a right-justified string of length width.

Padding is done using the specified fill character (default is a space).

rpartition(sep, /)#

Partition the string into three parts using the given separator.

This will search for the separator in the string, starting at the end. If the separator is found, returns a 3-tuple containing the part before the separator, the separator itself, and the part after it.

If the separator is not found, returns a 3-tuple containing two empty strings and the original string.

rsplit(sep=None, maxsplit=-1)#

Return a list of the substrings in the string, using sep as the separator string.

sep

The separator used to split the string.

When set to None (the default value), will split on any whitespace character (including n r t f and spaces) and will discard empty strings from the result.

maxsplit

Maximum number of splits. -1 (the default value) means no limit.

Splitting starts at the end of the string and works to the front.

rstrip(chars=None, /)#

Return a copy of the string with trailing whitespace removed.

If chars is given and not None, remove characters in chars instead.

split(sep=None, maxsplit=-1)#

Return a list of the substrings in the string, using sep as the separator string.

sep

The separator used to split the string.

When set to None (the default value), will split on any whitespace character (including n r t f and spaces) and will discard empty strings from the result.

maxsplit

Maximum number of splits. -1 (the default value) means no limit.

Splitting starts at the front of the string and works to the end.

Note, str.split() is mainly useful for data that has been intentionally delimited. With natural text that includes punctuation, consider using the regular expression module.

splitlines(keepends=False)#

Return a list of the lines in the string, breaking at line boundaries.

Line breaks are not included in the resulting list unless keepends is given and true.

startswith(prefix[, start[, end]]) → bool#

Return True if S starts with the specified prefix, False otherwise. With optional start, test S beginning at that position. With optional end, stop comparing S at that position. prefix can also be a tuple of strings to try.

strip(chars=None, /)#

Return a copy of the string with leading and trailing whitespace removed.

If chars is given and not None, remove characters in chars instead.

swapcase()#

Convert uppercase characters to lowercase and lowercase characters to uppercase.

title()#

Return a version of the string where each word is titlecased.

More specifically, words start with uppercased characters and all remaining cased characters have lower case.

translate(table, /)#

Replace each character in the string using the given translation table.

table

Translation table, which must be a mapping of Unicode ordinals to Unicode ordinals, strings, or None.

The table must implement lookup/indexing via __getitem__, for instance a dictionary or list. If this operation raises LookupError, the character is left untouched. Characters mapped to None are deleted.

upper()#

Return a copy of the string converted to uppercase.

zfill(width, /)#

Pad a numeric string with zeros on the left, to fill a field of the given width.

The string is never truncated.

pulse.boundary_conditions.nanson(F, N)[source]#

Push the unit normal N through F: \(\mathbf{n}\, da = J \mathbf{F}^{-T} \mathbf{N}\, dA\). Returns the unit normal \(\mathbf{n}\) and the area ratio \(da / dA\).

Parameters:

units#

class pulse.units.Variable(value: T, unit: pint.registry.Unit | str | None)[source]#
Parameters:

  • value (T)

  • unit (Unit | str | None)

utils#

pulse.utils.evaluate_at_vertex_tag(u, vt, tag)[source]#

Given a function u and a vertex tag vt return the values of u at the vertices tagged with tag

Parameters:
Returns:

The values of u at the vertices tagged with tag

Return type:

npt.NDArray[np.int32]

pulse.utils.gather_broadcast_array(comm, local_array)[source]#

Collects local arrays from all processes on the root process and distributes the global array to all processes. Assumes that the local arrays are either the same of empty / None on all processes.

Parameters:

  • comm (MPI.Comm) – The MPI communicator

  • local_array (np.ndarray) – The local array

Returns:

The global array on the root process

Return type:

np.ndarray

pulse.utils.map_vector_field(f, u=None, normalize=True, name='mapped_f0')[source]#

Map a vector field defined on a mesh to a new mesh.

Parameters:

  • f (dolfinx.fem.Function) – The vector field to map.

  • new_mesh (dolfinx.mesh.Mesh) – The new mesh to map the vector field to.

  • u (dolfinx.fem.Function | None, optional) – The displacement field used to map the points, by default None.

  • normalize (bool, optional) – Whether to normalize the vector field after mapping, by default True.

  • name (str)

Returns:

The mapped vector field on the new mesh.

Return type:

dolfinx.fem.Function

References#

[GMW91]

Julius M Guccione, Andrew D McCulloch, and Lewis K Waldman. Passive material properties of intact ventricular myocardium determined from a cylindrical model. Journal of Biomechanical Engineering, 113(1):42–55, 1991. doi:10.1115/1.2894084.

[HO09]

Gerhard A Holzapfel and Ray W Ogden. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367(1902):3445–3475, 2009. doi:10.1098/rsta.2009.0091.

[RQ21] (1,2,3)

Francesco Regazzoni and Alfio Quarteroni. An oscillation-free fully staggered algorithm for velocity-dependent active models of cardiac mechanics. Computer Methods in Applied Mechanics and Engineering, 373:113506, 2021. doi:10.1016/j.cma.2020.113506.

[ULM02]

Taras P Usyk, Ian J LeGrice, and Andrew D McCulloch. Computational model of three-dimensional cardiac electromechanics. Computing and visualization in science, 4(4):249–257, 2002. doi:10.1007/s00791-002-0081-9.

Contents