cbcbeat package

1. cbcbeat package#

1.1. Subpackages#

1.1.1. cbcbeat.cellmodels package

1.2. Submodules#

1.3. cbcbeat.bidomainsolver module#

These solvers solve the (pure) bidomain equations on the form: find the transmembrane potential \(v = v(x, t)\) and the extracellular potential \(u = u(x, t)\) such that

\[ \begin{align}\begin{aligned}v_t - \mathrm{div} ( G_i v + G_i u) = I_s\\\mathrm{div} (G_i v + (G_i + G_e) u) = I_a\end{aligned}\end{align} \]

where the subscript \(t\) denotes the time derivative; \(G_x\) denotes a weighted gradient: \(G_x = M_x \mathrm{grad}(v)\) for \(x \in \{i, e\}\), where \(M_i\) and \(M_e\) are the intracellular and extracellular cardiac conductivity tensors, respectively; \(I_s\) and \(I_a\) are prescribed input. In addition, initial conditions are given for \(v\):

\[v(x, 0) = v_0\]

Finally, boundary conditions must be prescribed. For now, this solver assumes pure homogeneous Neumann boundary conditions for \(v\) and \(u\) and enforces the additional average value zero constraint for u.

class cbcbeat.bidomainsolver.BasicBidomainSolver(mesh, time, M_i, M_e, I_s=None, I_a=None, v_=None, params=None)[source]#

Bases: object

This solver is based on a theta-scheme discretization in time and CG_1 x CG_1 (x R) elements in space.

Note

For the sake of simplicity and consistency with other solver objects, this solver operates on its solution fields (as state variables) directly internally. More precisely, solve (and step) calls will act by updating the internal solution fields. It implies that initial conditions can be set (and are intended to be set) by modifying the solution fields prior to simulation.

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
M_i (ufl.Expr): The intracellular conductivity tensor (as an UFL expression)
M_e (ufl.Expr): The extracellular conductivity tensor (as an UFL expression)
I_s (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
I_a (dolfin.Expression, optional): A (typically time-dependent) external applied current
v_ (ufl.Expr, optional): Initial condition for v. A new dolfin.Function will be created if none is given.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(BasicBidomainSolver.default_parameters(), True)

solution_fields()[source]#

Return tuple of previous and current solution objects.

Modifying these will modify the solution objects of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous v, current vur) (tuple of dolfin.Function)

solve(interval, dt=None)[source]#

Solve the discretization on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the interval and the current solution.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, optional): The timestep for the solve. Defaults to length of interval

Returns

(timestep, solution_fields) via (genexpr)

Example of usage:

# Create generator
solutions = solver.solve((0.0, 1.0), 0.1)

# Iterate over generator (computes solutions as you go)
for (interval, solution_fields) in solutions:
  (t0, t1) = interval
  v_, vur = solution_fields
  # do something with the solutions

step(interval)[source]#

Solve on the given time interval (t0, t1).

Arguments

interval (tuple): The time interval (t0, t1) for the step

Invariants

Assuming that v_ is in the correct state for t0, gives self.vur in correct state at t1.

property time#: The internal time of the solver.

class cbcbeat.bidomainsolver.BidomainSolver(mesh, time, M_i, M_e, I_s=None, I_a=None, v_=None, params=None)[source]#

Bases: BasicBidomainSolver

This solver is based on a theta-scheme discretization in time and CG_1 x CG_1 (x R) elements in space.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
M_i (ufl.Expr): The intracellular conductivity tensor (as an UFL expression)
M_e (ufl.Expr): The extracellular conductivity tensor (as an UFL expression)
I_s (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
I_a (dolfin.Expression, optional): A (typically time-dependent) external applied current
v_ (ufl.Expr, optional): Initial condition for v. A new dolfin.Function will be created if none is given.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(BidomainSolver.default_parameters(), True)

property linear_solver#: The linear solver (dolfin.LUSolver or dolfin.PETScKrylovSolver).

property nullspace#

step(interval)[source]#

Solve on the given time step (t0, t1).

Arguments

interval (tuple): The time interval (t0, t1) for the step

Invariants

Assuming that v_ is in the correct state for t0, gives self.vur in correct state at t1.

variational_forms(k_n)[source]#

Create the variational forms corresponding to the given discretization of the given system of equations.

Arguments

k_n (ufl.Expr or float): The time step

Returns

(lhs, rhs) (tuple of ufl.Form)

1.4. cbcbeat.cardiacmodels module#

This module contains a container class for cardiac models: CardiacModel. This class should be instantiated for setting up specific cardiac simulation scenarios.

class cbcbeat.cardiacmodels.CardiacModel(domain, time, M_i, M_e, cell_models, stimulus=None, applied_current=None)[source]#

Bases: object

A container class for cardiac models. Objects of this class represent a specific cardiac simulation set-up and should provide

A computational domain
A cardiac cell model
Intra-cellular and extra-cellular conductivities
Various forms of stimulus (optional).

This container class is designed for use with the splitting solvers (cbcbeat.splittingsolver), see their documentation for more information on how the attributes are interpreted in that context.

Arguments

domain (dolfin.Mesh): the computational domain in space
time (dolfin.Constant or None ): A constant holding the current time.
M_i (ufl.Expr): the intra-cellular conductivity as an ufl Expression
M_e (ufl.Expr): the extra-cellular conductivity as an ufl Expression
cell_models (CardiacCellModel): a cell model or a dict with cell models associated with a cell model domain
stimulus (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
applied_current (ufl.Expr, optional): an applied current as an ufl Expression

applied_current()[source]#: An applied current: used as a source in the elliptic bidomain equation

cell_models()[source]#: Return the cell models

conductivities()[source]#

Return the intracellular and extracellular conductivities as a tuple of UFL Expressions.

Returns (M_i, M_e) (tuple of ufl.Expr)

domain()[source]#: The spatial domain (dolfin.Mesh).

extracellular_conductivity()[source]#: The intracellular conductivity (ufl.Expr).

intracellular_conductivity()[source]#: The intracellular conductivity (ufl.Expr).

stimulus()[source]#: A stimulus: used as a source in the parabolic bidomain equation

time()[source]#: The current time (dolfin.Constant or None).

1.5. cbcbeat.cellsolver module#

This module contains solvers for (subclasses of) CardiacCellModel.

class cbcbeat.cellsolver.BasicCardiacODESolver(mesh, time, model, I_s=None, params=None)[source]#

Bases: object

A basic, non-optimised solver for systems of ODEs typically encountered in cardiac applications of the form: find a scalar field \(v = v(x, t)\) and a vector field \(s = s(x, t)\)

\[ \begin{align}\begin{aligned}v_t = - I_{ion}(v, s) + I_s\\s_t = F(v, s)\end{aligned}\end{align} \]

where \(I_{ion}\) and \(F\) are given non-linear functions, and \(I_s\) is some prescribed stimulus.

Here, this nonlinear ODE system is solved via a theta-scheme. By default theta=0.5, which corresponds to a Crank-Nicolson scheme. This can be changed by modifying the solver parameters.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
model (cbcbeat.CardiacCellModel): A representation of the cardiac cell model(s)
I_s (optional) A typically time-dependent external stimulus: given as a dolfin.cpp.function.GenericFunction or a Markerwise. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

solution_fields()[source]#

Return tuple of previous and current solution objects.

Modifying these will modify the solution objects of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous vs, current vs) (tuple of dolfin.Function)

solve(interval, dt=None)[source]#

Solve the problem given by the model on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the interval and the current vs solution.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, optional): The timestep for the solve. Defaults to length of interval

Returns

(timestep, current vs) via (genexpr)

Example of usage:

# Create generator
solutions = solver.solve((0.0, 1.0), 0.1)

# Iterate over generator (computes solutions as you go)
for (interval, vs) in solutions:
  # do something with the solutions

step(interval)[source]#

Solve on the given time step (t0, t1).

End users are recommended to use solve instead.

Arguments

interval (tuple): The time interval (t0, t1) for the step

property time#: The internal time of the solver.

class cbcbeat.cellsolver.BasicSingleCellSolver(model, time, params=None)[source]#

Bases: BasicCardiacODESolver

A basic, non-optimised solver for systems of ODEs typically encountered in cardiac applications of the form: find a scalar field \(v = v(t)\) and a vector field \(s = s(t)\)

\[ \begin{align}\begin{aligned}v_t = - I_{ion}(v, s) + I_s\\s_t = F(v, s)\end{aligned}\end{align} \]

where \(I_{ion}\) and \(F\) are given non-linear functions, \(I_s\) is some prescribed stimulus. If \(I_s\) depends on time, it is assumed that \(I_s\) is a dolfin.Expression with parameter ‘t’.

Use this solver if you just want to test the results from a cardiac cell model without any spatial mesh dependence.

Here, this nonlinear ODE system is solved via a theta-scheme. By default theta=0.5, which corresponds to a Crank-Nicolson scheme. This can be changed by modifying the solver parameters.

Note

Arguments

model (CardiacCellModel): A cardiac cell model
time (Constant or None): A constant holding the current time.
params (dolfin.Parameters, optional): Solver parameters

class cbcbeat.cellsolver.CardiacODESolver(mesh, time, model, I_s=None, params=None)[source]#

Bases: object

An optimised solver for systems of ODEs typically encountered in cardiac applications of the form: find a scalar field \(v = v(x, t)\) and a vector field \(s = s(x, t)\)

\[ \begin{align}\begin{aligned}v_t = - I_{ion}(v, s) + I_s\\s_t = F(v, s)\end{aligned}\end{align} \]

where \(I_{ion}\) and \(F\) are given non-linear functions, and \(I_s\) is some prescribed stimulus.

Note

Arguments

mesh (dolfin.Mesh): The spatial mesh (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
model (cbcbeat.CardiacCellModel): A representation of the cardiac cell model(s)
I_s (dolfin.Expression, optional): A typically time-dependent external stimulus. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

solution_fields()[source]#

Return current solution object.

Modifying this will modify the solution object of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous vs_, current vs) (dolfin.Function)

solve(interval, dt=None)[source]#

Solve the problem given by the model on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the interval and the current vs solution.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, optional): The timestep for the solve. Defaults to length of interval

Returns

(timestep, current vs) via (genexpr)

Example of usage:

# Create generator
solutions = solver.solve((0.0, 1.0), 0.1)

# Iterate over generator (computes solutions as you go)
for (interval, vs) in solutions:
  # do something with the solutions

step(interval)[source]#

Solve on the given time step (t0, t1).

End users are recommended to use solve instead.

Arguments

interval (tuple): The time interval (t0, t1) for the step

class cbcbeat.cellsolver.SingleCellSolver(model, time, params=None)[source]#: Bases: CardiacODESolver

1.6. cbcbeat.dolfinimport module#

This module handles all dolfin import in cbcbeat. Here dolfin and dolfin_adjoint gets imported. If dolfin_adjoint is not present it will not be imported.

1.7. cbcbeat.gossplittingsolver module#

1.8. cbcbeat.gotran2cellmodel module#

1.9. cbcbeat.markerwisefield module#

class cbcbeat.markerwisefield.Markerwise(objects, keys, markers)[source]#

Bases: object

A container class representing an object defined by a number of objects combined with a mesh function defining mesh domains and a map between the these.

Arguments

objects (tuple): the different objects
keys (tuple of ints): a map from the objects to the domains marked in markers
markers (dolfin.MeshFunction): a mesh function mapping which domains the mesh consist of

Example of usage

Given (g0, g1), (2, 5) and markers, let

g = g0 on domains marked by 2 in markers g = g1 on domains marked by 5 in markers

letting:

g = Markerwise((g0, g1), (2, 5), markers)

keys()[source]#: The keys or domain numbers

markers()[source]#: The markers

values()[source]#: The objects

cbcbeat.markerwisefield.handle_markerwise(g, classtype)[source]#

cbcbeat.markerwisefield.rhs_with_markerwise_field(g, mesh, v)[source]#

1.10. cbcbeat.monodomainsolver module#

These solvers solve the (pure) monodomain equations on the form: find the transmembrane potential \(v = v(x, t)\) such that

\[v_t - \mathrm{div} ( G_i v) = I_s\]

where the subscript \(t\) denotes the time derivative; \(G_i\) denotes a weighted gradient: \(G_i = M_i \mathrm{grad}(v)\) for, where \(M_i\) is the intracellular cardiac conductivity tensor; \(I_s\) ise prescribed input. In addition, initial conditions are given for \(v\):

\[v(x, 0) = v_0\]

Finally, boundary conditions must be prescribed. For now, this solver assumes pure homogeneous Neumann boundary conditions for \(v\).

class cbcbeat.monodomainsolver.BasicMonodomainSolver(mesh, time, M_i, I_s=None, v_=None, params=None)[source]#

Bases: object

This solver is based on a theta-scheme discretization in time and CG_1 elements in space.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
M_i (ufl.Expr): The intracellular conductivity tensor (as an UFL expression)
I_s (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
v_ (ufl.Expr, optional): Initial condition for v. A new dolfin.Function will be created if none is given.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(BasicMonodomainSolver.default_parameters(), True)

solution_fields()[source]#

Return tuple of previous and current solution objects.

Modifying these will modify the solution objects of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous v, current v) (tuple of dolfin.Function)

solve(interval, dt=None)[source]#

Solve the discretization on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the interval and the current solution.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, optional): The timestep for the solve. Defaults to length of interval

Returns

(timestep, solution_field) via (genexpr)

Example of usage:

# Create generator
solutions = solver.solve((0.0, 1.0), 0.1)

# Iterate over generator (computes solutions as you go)
for (interval, solution_fields) in solutions:
  (t0, t1) = interval
  v_, v = solution_fields
  # do something with the solutions

step(interval)[source]#

Solve on the given time interval (t0, t1).

Arguments

interval (tuple): The time interval (t0, t1) for the step

Invariants

Assuming that v_ is in the correct state for t0, gives self.v in correct state at t1.

property time#: The internal time of the solver.

class cbcbeat.monodomainsolver.MonodomainSolver(mesh, time, M_i, I_s=None, v_=None, params=None)[source]#

Bases: BasicMonodomainSolver

This solver is based on a theta-scheme discretization in time and CG_1 elements in space.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
M_i (ufl.Expr): The intracellular conductivity tensor (as an UFL expression)
I_s (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
v_ (ufl.Expr, optional): Initial condition for v. A new dolfin.Function will be created if none is given.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(MonodomainSolver.default_parameters(), True)

property linear_solver#: The linear solver (dolfin.LUSolver or dolfin.KrylovSolver).

step(interval)[source]#

Solve on the given time step (t0, t1).

Arguments

interval (tuple): The time interval (t0, t1) for the step

Invariants

Assuming that v_ is in the correct state for t0, gives self.v in correct state at t1.

variational_forms(k_n)[source]#

Create the variational forms corresponding to the given discretization of the given system of equations.

Arguments

k_n (ufl.Expr or float): The time step

Returns

(lhs, rhs, prec) (tuple of ufl.Form)

1.11. cbcbeat.splittingsolver module#

This module contains splitting solvers for CardiacModel objects. In particular, the classes

SplittingSolver

BasicSplittingSolver

These solvers solve the bidomain (or monodomain) equations on the form: find the transmembrane potential \(v = v(x, t)\) in mV, the extracellular potential \(u = u(x, t)\) in mV, and any additional state variables \(s = s(x, t)\) such that

\[ \begin{align}\begin{aligned}v_t - \mathrm{div} (M_i \mathrm{grad} v + M_i \mathrm{grad} u) = - I_{ion}(v, s) + I_s\\ \mathrm{div} (M_i \mathrm{grad} v + (M_i + M_e) \mathrm{grad} u) = I_a\\s_t = F(v, s)\end{aligned}\end{align} \]

where

the subscript \(t\) denotes the time derivative,

\(M_i\) and \(M_e\) are conductivity tensors (in mm^2/ms)

\(I_s\) is prescribed input current (in mV/ms)

\(I_a\) is prescribed input current (in mV/ms)

\(I_{ion}\) and \(F\) are typically specified by a cell model

Note that M_i and M_e can be viewed as scaled by \(\chi*C_m\) where

\(\chi\) is the surface-to volume ratio of cells (in 1/mm) ,
\(C_m\) is the specific membrane capacitance (in mu F/(mm^2) ),

In addition, initial conditions are given for \(v\) and \(s\):

\[ \begin{align}\begin{aligned}v(x, 0) = v_0\\s(x, 0) = s_0\end{aligned}\end{align} \]

Finally, boundary conditions must be prescribed. These solvers assume pure Neumann boundary conditions for \(v\) and \(u\) and enforce the additional average value zero constraint for u.

The solvers take as input a CardiacModel providing the required input specification of the problem. In particular, the applied current \(I_a\) is extracted from the applied_current attribute, while the stimulus \(I_s\) is extracted from the stimulus attribute.

It should be possible to use the solvers interchangably. However, note that the BasicSplittingSolver is not optimised and should be used for testing or debugging purposes primarily.

class cbcbeat.splittingsolver.BasicSplittingSolver(model, params=None, V_index=0)[source]#

Bases: object

A non-optimised solver for the bidomain equations based on the operator splitting scheme described in Sundnes et al 2006, p. 78 ff.

The solver computes as solutions:

“vs” (dolfin.Function) representing the solution for the transmembrane potential and any additional state variables, and

“vur” (dolfin.Function) representing the transmembrane potential in combination with the extracellular potential and an additional Lagrange multiplier.

The algorithm can be controlled by a number of parameters. In particular, the splitting algorithm can be controlled by the parameter “theta”: “theta” set to 1.0 corresponds to a (1st order) Godunov splitting while “theta” set to 0.5 to a (2nd order) Strang splitting.

This solver has not been optimised for computational efficiency and should therefore primarily be used for debugging purposes. For an equivalent, but more efficient, solver, see cbcbeat.splittingsolver.SplittingSolver.

Arguments

model (cbcbeat.cardiacmodels.CardiacModel): a CardiacModel object describing the simulation set-up
params (dolfin.Parameters, optional): a Parameters object controlling solver parameters
V_index (int): Index for the membrane potential among the state variables, by default 0

Assumptions

The cardiac conductivities do not vary in time

static default_parameters()[source]#

Initialize and return a set of default parameters for the splitting solver

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(BasicSplittingSolver.default_parameters(), True)

merge(solution)[source]#

Combine solutions from the PDE solve and the ODE solve to form a single mixed function.

Arguments

solution (dolfin.Function): Function holding the combined result

solution_fields()[source]#

Return tuple of previous and current solution objects.

Modifying these will modify the solution objects of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous vs, current vs, current vur) (tuple of dolfin.Function)

solve(interval, dt)[source]#

Solve the problem given by the model on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the time step and the solution fields.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, list of tuples of floats): The timestep for the solve. A list of tuples of floats can also be passed. Each tuple should contain two floats where the first includes the start time and the second the dt.

Returns

(timestep, solution_fields) via (genexpr)

Example of usage:

# Create generator
dts = [(0., 0.1), (1.0, 0.05), (2.0, 0.1)]
solutions = solver.solve((0.0, 1.0), dts)

# Iterate over generator (computes solutions as you go)
for ((t0, t1), (vs_, vs, vur)) in solutions:
  # do something with the solutions

step(interval)[source]#

Solve the problem given by the model on a given time interval (t0, t1) with timestep given by the interval length.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)

Invariants

Given self._vs in a correct state at t0, provide v and s (in self.vs) and u (in self.vur) in a correct state at t1. (Note that self.vur[0] == self.vs[0] only if theta = 1.0.)

class cbcbeat.splittingsolver.SplittingSolver(model, params=None, V_index=0)[source]#

Bases: BasicSplittingSolver

An optimised solver for the bidomain equations based on the operator splitting scheme described in Sundnes et al 2006, p. 78 ff.

The solver computes as solutions:

“vs” (dolfin.Function) representing the solution for the transmembrane potential and any additional state variables, and

“vur” (dolfin.Function) representing the transmembrane potential in combination with the extracellular potential and an additional Lagrange multiplier.

Arguments

model (cbcbeat.cardiacmodels.CardiacModel): a CardiacModel object describing the simulation set-up
params (dolfin.Parameters, optional): a Parameters object controlling solver parameters

Example of usage:

from cbcbeat import *

# Describe the cardiac model
mesh = UnitSquareMesh(100, 100)
time = Constant(0.0)
cell_model = FitzHughNagumoManual()
stimulus = Expression("10*t*x[0]", t=time, degree=1)
cm = CardiacModel(mesh, time, 1.0, 1.0, cell_model, stimulus)

# Extract default solver parameters
ps = SplittingSolver.default_parameters()

# Examine the default parameters
info(ps, True)

# Update parameter dictionary
ps["theta"] = 1.0 # Use first order splitting
ps["CardiacODESolver"]["scheme"] = "GRL1" # Use Generalized Rush-Larsen scheme

# Use monodomain equations as the PDE model
ps["pde_solver"] = "monodomain"
# Use direct linear solver of the PDEs
ps["MonodomainSolver"]["linear_solver_type"] = "direct"
# Use backward Euler for temporal discretization for the PDEs
ps["MonodomainSolver"]["theta"] = 1.0

solver = SplittingSolver(cm, params=ps)

# Extract the solution fields and set the initial conditions
(vs_, vs, vur) = solver.solution_fields()
vs_.assign(cell_model.initial_conditions())

# Solve
dt = 0.1
T = 1.0
interval = (0.0, T)
for (timestep, fields) in solver.solve(interval, dt):
    (vs_, vs, vur) = fields
    # Do something with the solutions

Assumptions

The cardiac conductivities do not vary in time

static default_parameters()[source]#

Initialize and return a set of default parameters for the splitting solver

Returns: The set of default parameters (dolfin.Parameters)

Example of usage:

info(SplittingSolver.default_parameters(), True)

1.12. cbcbeat.utils module#

This module provides various utilities for internal use.

class cbcbeat.utils.Projecter(V, params=None)[source]#

Bases: object

Customized class for repeated projection.

Arguments

V (dolfin.FunctionSpace): The function space to project into
solver_type (string, optional): “iterative” (default) or “direct”

Example of usage::

my_project = Projecter(V, solver_type=”direct”) u = Function(V) f = Function(W) my_project(f, u)

static default_parameters()[source]#

cbcbeat.utils.convergence_rate(hs, errors)[source]#: Compute and return rates of convergence \(r_i\) such that

\[errors = C hs^r\]

cbcbeat.utils.end_of_time(T, t0, t1, dt)[source]#: Return True if the interval (t0, t1) is the last before the end time T, otherwise False.

cbcbeat.utils.state_space(domain, d, family=None, k=1)[source]#

Return function space for the state variables.

Arguments

domain (dolfin.Mesh): The computational domain
d (int): The number of states
family (string, optional): The finite element family, defaults to “CG” if None is given.
k (int, optional): The finite element degree, defaults to 1

Returns

a function space (dolfin.FunctionSpace)

1.13. Module contents#

The cbcbeat Python module is a problem and solver collection for cardiac electrophysiology models.

To import the module, type:

from cbcbeat import *

class cbcbeat.AbstractCell(topological_dimension, geometric_dimension)[source]#

Bases: object

Representation of an abstract finite element cell with only the dimensions known.

geometric_dimension()[source]#: Return the dimension of the space this cell is embedded in.

has_simplex_facets()[source]#: Return True if all the facets of this cell are simplex cells.

is_simplex()[source]#: Return True if this is a simplex cell.

topological_dimension()[source]#: Return the dimension of the topology of this cell.

class cbcbeat.AbstractDomain(topological_dimension, geometric_dimension)[source]#

Bases: object

Symbolic representation of a geometric domain with only a geometric and topological dimension.

geometric_dimension()[source]#: Return the dimension of the space this domain is embedded in.

topological_dimension()[source]#: Return the dimension of the topology of this domain.

cbcbeat.And(left, right)[source]#: UFL operator: A boolean expression (left and right) for use with conditional.

class cbcbeat.Argument(function_space, number, part=None)[source]#

Bases: FormArgument

UFL value: Representation of an argument to a form.

is_cellwise_constant()[source]#: Return whether this expression is spatially constant over each cell.

number()[source]#: Return the Argument number.

part()[source]#

ufl_domain()[source]#: Deprecated, please use .ufl_function_space().ufl_domain() instead.

ufl_domains()[source]#: Deprecated, please use .ufl_function_space().ufl_domains() instead.

ufl_element()[source]#: Deprecated, please use .ufl_function_space().ufl_element() instead.

ufl_function_space()[source]#: Get the function space of this Argument.

property ufl_shape#: Return the associated UFL shape.

cbcbeat.Arguments(function_space, number)[source]#: UFL value: Create an Argument in a mixed space, and return a tuple with the function components corresponding to the subelements.

class cbcbeat.BasicBidomainSolver(mesh, time, M_i, M_e, I_s=None, I_a=None, v_=None, params=None)[source]#

Bases: object

This solver is based on a theta-scheme discretization in time and CG_1 x CG_1 (x R) elements in space.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
M_i (ufl.Expr): The intracellular conductivity tensor (as an UFL expression)
M_e (ufl.Expr): The extracellular conductivity tensor (as an UFL expression)
I_s (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
I_a (dolfin.Expression, optional): A (typically time-dependent) external applied current
v_ (ufl.Expr, optional): Initial condition for v. A new dolfin.Function will be created if none is given.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(BasicBidomainSolver.default_parameters(), True)

solution_fields()[source]#

Return tuple of previous and current solution objects.

Modifying these will modify the solution objects of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous v, current vur) (tuple of dolfin.Function)

solve(interval, dt=None)[source]#

Solve the discretization on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the interval and the current solution.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, optional): The timestep for the solve. Defaults to length of interval

Returns

(timestep, solution_fields) via (genexpr)

Example of usage:

# Create generator
solutions = solver.solve((0.0, 1.0), 0.1)

# Iterate over generator (computes solutions as you go)
for (interval, solution_fields) in solutions:
  (t0, t1) = interval
  v_, vur = solution_fields
  # do something with the solutions

step(interval)[source]#

Solve on the given time interval (t0, t1).

Arguments

interval (tuple): The time interval (t0, t1) for the step

Invariants

Assuming that v_ is in the correct state for t0, gives self.vur in correct state at t1.

property time#: The internal time of the solver.

class cbcbeat.BasicCardiacODESolver(mesh, time, model, I_s=None, params=None)[source]#

Bases: object

A basic, non-optimised solver for systems of ODEs typically encountered in cardiac applications of the form: find a scalar field \(v = v(x, t)\) and a vector field \(s = s(x, t)\)

\[ \begin{align}\begin{aligned}v_t = - I_{ion}(v, s) + I_s\\s_t = F(v, s)\end{aligned}\end{align} \]

where \(I_{ion}\) and \(F\) are given non-linear functions, and \(I_s\) is some prescribed stimulus.

Here, this nonlinear ODE system is solved via a theta-scheme. By default theta=0.5, which corresponds to a Crank-Nicolson scheme. This can be changed by modifying the solver parameters.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
model (cbcbeat.CardiacCellModel): A representation of the cardiac cell model(s)
I_s (optional) A typically time-dependent external stimulus: given as a dolfin.cpp.function.GenericFunction or a Markerwise. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

solution_fields()[source]#

Return tuple of previous and current solution objects.

Modifying these will modify the solution objects of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous vs, current vs) (tuple of dolfin.Function)

solve(interval, dt=None)[source]#

Solve the problem given by the model on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the interval and the current vs solution.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, optional): The timestep for the solve. Defaults to length of interval

Returns

(timestep, current vs) via (genexpr)

Example of usage:

# Create generator
solutions = solver.solve((0.0, 1.0), 0.1)

# Iterate over generator (computes solutions as you go)
for (interval, vs) in solutions:
  # do something with the solutions

step(interval)[source]#

Solve on the given time step (t0, t1).

End users are recommended to use solve instead.

Arguments

interval (tuple): The time interval (t0, t1) for the step

property time#: The internal time of the solver.

class cbcbeat.BasicMonodomainSolver(mesh, time, M_i, I_s=None, v_=None, params=None)[source]#

Bases: object

This solver is based on a theta-scheme discretization in time and CG_1 elements in space.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
M_i (ufl.Expr): The intracellular conductivity tensor (as an UFL expression)
I_s (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
v_ (ufl.Expr, optional): Initial condition for v. A new dolfin.Function will be created if none is given.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(BasicMonodomainSolver.default_parameters(), True)

solution_fields()[source]#

Return tuple of previous and current solution objects.

Modifying these will modify the solution objects of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous v, current v) (tuple of dolfin.Function)

solve(interval, dt=None)[source]#

Solve the discretization on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the interval and the current solution.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, optional): The timestep for the solve. Defaults to length of interval

Returns

(timestep, solution_field) via (genexpr)

Example of usage:

# Create generator
solutions = solver.solve((0.0, 1.0), 0.1)

# Iterate over generator (computes solutions as you go)
for (interval, solution_fields) in solutions:
  (t0, t1) = interval
  v_, v = solution_fields
  # do something with the solutions

step(interval)[source]#

Solve on the given time interval (t0, t1).

Arguments

interval (tuple): The time interval (t0, t1) for the step

Invariants

Assuming that v_ is in the correct state for t0, gives self.v in correct state at t1.

property time#: The internal time of the solver.

class cbcbeat.BasicSingleCellSolver(model, time, params=None)[source]#

Bases: BasicCardiacODESolver

A basic, non-optimised solver for systems of ODEs typically encountered in cardiac applications of the form: find a scalar field \(v = v(t)\) and a vector field \(s = s(t)\)

\[ \begin{align}\begin{aligned}v_t = - I_{ion}(v, s) + I_s\\s_t = F(v, s)\end{aligned}\end{align} \]

Use this solver if you just want to test the results from a cardiac cell model without any spatial mesh dependence.

Here, this nonlinear ODE system is solved via a theta-scheme. By default theta=0.5, which corresponds to a Crank-Nicolson scheme. This can be changed by modifying the solver parameters.

Note

Arguments

model (CardiacCellModel): A cardiac cell model
time (Constant or None): A constant holding the current time.
params (dolfin.Parameters, optional): Solver parameters

class cbcbeat.BasicSplittingSolver(model, params=None, V_index=0)[source]#

Bases: object

A non-optimised solver for the bidomain equations based on the operator splitting scheme described in Sundnes et al 2006, p. 78 ff.

The solver computes as solutions:

“vs” (dolfin.Function) representing the solution for the transmembrane potential and any additional state variables, and

“vur” (dolfin.Function) representing the transmembrane potential in combination with the extracellular potential and an additional Lagrange multiplier.

Arguments

model (cbcbeat.cardiacmodels.CardiacModel): a CardiacModel object describing the simulation set-up
params (dolfin.Parameters, optional): a Parameters object controlling solver parameters
V_index (int): Index for the membrane potential among the state variables, by default 0

Assumptions

The cardiac conductivities do not vary in time

static default_parameters()[source]#

Initialize and return a set of default parameters for the splitting solver

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(BasicSplittingSolver.default_parameters(), True)

merge(solution)[source]#

Combine solutions from the PDE solve and the ODE solve to form a single mixed function.

Arguments

solution (dolfin.Function): Function holding the combined result

solution_fields()[source]#

Return tuple of previous and current solution objects.

Modifying these will modify the solution objects of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous vs, current vs, current vur) (tuple of dolfin.Function)

solve(interval, dt)[source]#

Solve the problem given by the model on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the time step and the solution fields.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, list of tuples of floats): The timestep for the solve. A list of tuples of floats can also be passed. Each tuple should contain two floats where the first includes the start time and the second the dt.

Returns

(timestep, solution_fields) via (genexpr)

Example of usage:

# Create generator
dts = [(0., 0.1), (1.0, 0.05), (2.0, 0.1)]
solutions = solver.solve((0.0, 1.0), dts)

# Iterate over generator (computes solutions as you go)
for ((t0, t1), (vs_, vs, vur)) in solutions:
  # do something with the solutions

step(interval)[source]#

Solve the problem given by the model on a given time interval (t0, t1) with timestep given by the interval length.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)

Invariants

Given self._vs in a correct state at t0, provide v and s (in self.vs) and u (in self.vur) in a correct state at t1. (Note that self.vur[0] == self.vs[0] only if theta = 1.0.)

class cbcbeat.Beeler_reuter_1977(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

F(v, s, time=None)[source]#: Right hand side for ODE system

I(v, s, time=None)[source]#: Transmembrane current

I = -dV/dt

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

class cbcbeat.BidomainSolver(mesh, time, M_i, M_e, I_s=None, I_a=None, v_=None, params=None)[source]#

Bases: BasicBidomainSolver

This solver is based on a theta-scheme discretization in time and CG_1 x CG_1 (x R) elements in space.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
M_i (ufl.Expr): The intracellular conductivity tensor (as an UFL expression)
M_e (ufl.Expr): The extracellular conductivity tensor (as an UFL expression)
I_s (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
I_a (dolfin.Expression, optional): A (typically time-dependent) external applied current
v_ (ufl.Expr, optional): Initial condition for v. A new dolfin.Function will be created if none is given.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(BidomainSolver.default_parameters(), True)

property linear_solver#: The linear solver (dolfin.LUSolver or dolfin.PETScKrylovSolver).

property nullspace#

step(interval)[source]#

Solve on the given time step (t0, t1).

Arguments

interval (tuple): The time interval (t0, t1) for the step

Invariants

Assuming that v_ is in the correct state for t0, gives self.vur in correct state at t1.

variational_forms(k_n)[source]#

Create the variational forms corresponding to the given discretization of the given system of equations.

Arguments

k_n (ufl.Expr or float): The time step

Returns

(lhs, rhs) (tuple of ufl.Form)

class cbcbeat.BrokenElement(element)[source]#

Bases: FiniteElementBase

The discontinuous version of an existing Finite Element space.

mapping()[source]#: Not implemented.

reconstruct(**kwargs)[source]#

shortstr()[source]#: Format as string for pretty printing.

class cbcbeat.CardiacCellModel(params=None, init_conditions=None)[source]#

Bases: object

Base class for cardiac cell models. Specialized cell models should subclass this class.

Essentially, a cell model represents a system of ordinary differential equations. A cell model is here described by two (Python) functions, named F and I. The model describes the behavior of the transmembrane potential ‘v’ and a number of state variables ‘s’

The function F gives the right-hand side for the evolution of the state variables:

d/dt s = F(v, s)

The function I gives the ionic current. If a single cell is considered, I gives the (negative) right-hand side for the evolution of the transmembrane potential

(*) d/dt v = - I(v, s)

If used in a bidomain setting, the ionic current I enters into the parabolic partial differential equation of the bidomain equations.

If a stimulus is provided via

cell = CardiacCellModel() cell.stimulus = Expression(“I_s(t)”, degree=1)

then I_s is added to the right-hand side of (*), which thus reads

d/dt v = - I(v, s) + I_s

Note that the cardiac cell model stimulus is ignored when the cell model is used a spatially-varying setting (for instance in the bidomain setting). In this case, the user is expected to specify a stimulus for the cardiac model instead.

F(v, s, time=None)[source]#: Return right-hand side for state variable evolution.

I(v, s, time=None)[source]#: Return the ionic current.

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

initial_conditions()[source]#: Return initial conditions for v and s as an Expression.

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

parameters()[source]#: Return the current parameters.

set_initial_conditions(**init)[source]#: Update initial_conditions in model

set_parameters(**params)[source]#: Update parameters in model

class cbcbeat.CardiacModel(domain, time, M_i, M_e, cell_models, stimulus=None, applied_current=None)[source]#

Bases: object

A container class for cardiac models. Objects of this class represent a specific cardiac simulation set-up and should provide

A computational domain
A cardiac cell model
Intra-cellular and extra-cellular conductivities
Various forms of stimulus (optional).

This container class is designed for use with the splitting solvers (cbcbeat.splittingsolver), see their documentation for more information on how the attributes are interpreted in that context.

Arguments

domain (dolfin.Mesh): the computational domain in space
time (dolfin.Constant or None ): A constant holding the current time.
M_i (ufl.Expr): the intra-cellular conductivity as an ufl Expression
M_e (ufl.Expr): the extra-cellular conductivity as an ufl Expression
cell_models (CardiacCellModel): a cell model or a dict with cell models associated with a cell model domain
stimulus (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
applied_current (ufl.Expr, optional): an applied current as an ufl Expression

applied_current()[source]#: An applied current: used as a source in the elliptic bidomain equation

cell_models()[source]#: Return the cell models

conductivities()[source]#

Return the intracellular and extracellular conductivities as a tuple of UFL Expressions.

Returns (M_i, M_e) (tuple of ufl.Expr)

domain()[source]#: The spatial domain (dolfin.Mesh).

extracellular_conductivity()[source]#: The intracellular conductivity (ufl.Expr).

intracellular_conductivity()[source]#: The intracellular conductivity (ufl.Expr).

stimulus()[source]#: A stimulus: used as a source in the parabolic bidomain equation

time()[source]#: The current time (dolfin.Constant or None).

class cbcbeat.CardiacODESolver(mesh, time, model, I_s=None, params=None)[source]#

Bases: object

An optimised solver for systems of ODEs typically encountered in cardiac applications of the form: find a scalar field \(v = v(x, t)\) and a vector field \(s = s(x, t)\)

\[ \begin{align}\begin{aligned}v_t = - I_{ion}(v, s) + I_s\\s_t = F(v, s)\end{aligned}\end{align} \]

where \(I_{ion}\) and \(F\) are given non-linear functions, and \(I_s\) is some prescribed stimulus.

Note

Arguments

mesh (dolfin.Mesh): The spatial mesh (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
model (cbcbeat.CardiacCellModel): A representation of the cardiac cell model(s)
I_s (dolfin.Expression, optional): A typically time-dependent external stimulus. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

solution_fields()[source]#

Return current solution object.

Modifying this will modify the solution object of the solver and thus provides a way for setting initial conditions for instance.

Returns: (previous vs_, current vs) (dolfin.Function)

solve(interval, dt=None)[source]#

Solve the problem given by the model on a given time interval (t0, t1) with a given timestep dt and return generator for a tuple of the interval and the current vs solution.

Arguments

interval (tuple): The time interval for the solve given by (t0, t1)
dt (int, optional): The timestep for the solve. Defaults to length of interval

Returns

(timestep, current vs) via (genexpr)

Example of usage:

# Create generator
solutions = solver.solve((0.0, 1.0), 0.1)

# Iterate over generator (computes solutions as you go)
for (interval, vs) in solutions:
  # do something with the solutions

step(interval)[source]#

Solve on the given time step (t0, t1).

End users are recommended to use solve instead.

Arguments

interval (tuple): The time interval (t0, t1) for the step

class cbcbeat.Cell(cellname, geometric_dimension=None)[source]#

Bases: AbstractCell

Representation of a named finite element cell with known structure.

cellname()[source]#: Return the cellname of the cell.

facet_types()[source]#: A tuple of ufl_legacy.Cell representing the facets of self.

has_simplex_facets()[source]#: Return True if all the facets of this cell are simplex cells.

is_simplex()[source]#: Return True if this is a simplex cell.

num_edges()[source]#: The number of cell edges.

num_facets()[source]#: The number of cell facets.

num_vertices()[source]#: The number of cell vertices.

reconstruct(geometric_dimension=None)[source]#

class cbcbeat.CellDiameter(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The diameter of the cell, i.e., maximal distance of two points in the cell.

name = 'diameter'#

class cbcbeat.CellNormal(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The upwards pointing normal vector of the current manifold cell.

name = 'cell_normal'#

property ufl_shape#: Return the number of coordinates defined (i.e. the geometric dimension of the domain).

class cbcbeat.CellVolume(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The volume of the cell.

name = 'volume'#

class cbcbeat.Circumradius(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The circumradius of the cell.

name = 'circumradius'#

class cbcbeat.Coefficient(function_space, count=None)[source]#

Bases: FormArgument

UFL form argument type: Representation of a form coefficient.

count()[source]#

is_cellwise_constant()[source]#: Return whether this expression is spatially constant over each cell.

ufl_domain()[source]#: Shortcut to get the domain of the function space of this coefficient.

ufl_domains()[source]#: Return tuple of domains related to this terminal object.

ufl_element()[source]#: Shortcut to get the finite element of the function space of this coefficient.

ufl_function_space()[source]#: Get the function space of this coefficient.

property ufl_shape#: Return the associated UFL shape.

cbcbeat.Coefficients(function_space)[source]#: UFL value: Create a Coefficient in a mixed space, and return a tuple with the function components corresponding to the subelements.

class cbcbeat.Constant(domain, shape=(), count=None)[source]#

Bases: Terminal

count()[source]#

is_cellwise_constant()[source]#

ufl_domain()[source]#: Return the single unique domain this expression is defined on, or throw an error.

ufl_domains()[source]#: Return tuple of domains related to this terminal object.

property ufl_shape#

cbcbeat.Dn(f)[source]#: UFL operator: Take the directional derivative of f in the facet normal direction, Dn(f) := dot(grad(f), n).

cbcbeat.Dx(f, *i)[source]#: UFL operator: Take the partial derivative of f with respect to spatial variable number i. Equivalent to f.dx(*i).

class cbcbeat.EnrichedElement(*elements)[source]#

Bases: EnrichedElementBase

The vector sum of several finite element spaces:

\[\textrm{EnrichedElement}(V, Q) = \{v + q | v \in V, q \in Q\}.\]

Dual basis is a concatenation of subelements dual bases; primal basis is a concatenation of subelements primal bases; resulting element is not nodal even when subelements are. Structured basis may be exploited in form compilers.

is_cellwise_constant()[source]#: Return whether the basis functions of this element is spatially constant over each cell.

shortstr()[source]#: Format as string for pretty printing.

class cbcbeat.FacetArea(domain)[source]#

Bases: GeometricFacetQuantity

UFL geometry representation: The area of the facet.

name = 'facetarea'#

cbcbeat.FacetElement(element)[source]#: Constructs the restriction of a finite element to the facets of the cell.

class cbcbeat.FacetNormal(domain)[source]#

Bases: GeometricFacetQuantity

UFL geometry representation: The outwards pointing normal vector of the current facet.

is_cellwise_constant()[source]#: Return whether this expression is spatially constant over each cell.

name = 'n'#

property ufl_shape#: Return the number of coordinates defined (i.e. the geometric dimension of the domain).

class cbcbeat.Fenton_karma_1998_BR_altered(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

F(v, s, time=None)[source]#: Right hand side for ODE system

I(v, s, time=None)[source]#: Transmembrane current

I = -dV/dt

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

class cbcbeat.Fenton_karma_1998_MLR1_altered(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

F(v, s, time=None)[source]#: Right hand side for ODE system

I(v, s, time=None)[source]#: Transmembrane current

I = -dV/dt

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

class cbcbeat.FiniteElement(family, cell=None, degree=None, form_degree=None, quad_scheme=None, variant=None)[source]#

Bases: FiniteElementBase

The basic finite element class for all simple finite elements.

mapping()[source]#: Not implemented.

reconstruct(family=None, cell=None, degree=None, quad_scheme=None, variant=None)[source]#: Construct a new FiniteElement object with some properties replaced with new values.

shortstr()[source]#: Format as string for pretty printing.

sobolev_space()[source]#: Return the underlying Sobolev space.

variant()[source]#: Return the variant used to initialise the element.

class cbcbeat.FiniteElementBase(family, cell, degree, quad_scheme, value_shape, reference_value_shape)[source]#

Bases: object

Base class for all finite elements.

cell()[source]#: Return cell of finite element.

degree(component=None)[source]#: Return polynomial degree of finite element.

extract_component(i)[source]#: Recursively extract component index relative to a (simple) element and that element for given value component index.

extract_reference_component(i)[source]#: Recursively extract reference component index relative to a (simple) element and that element for given reference value component index.

extract_subelement_component(i)[source]#: Extract direct subelement index and subelement relative component index for a given component index.

extract_subelement_reference_component(i)[source]#: Extract direct subelement index and subelement relative reference component index for a given reference component index.

family()[source]#: Return finite element family.

is_cellwise_constant(component=None)[source]#: Return whether the basis functions of this element is spatially constant over each cell.

mapping()[source]#: Not implemented.

num_sub_elements()[source]#: Return number of sub-elements.

quadrature_scheme()[source]#: Return quadrature scheme of finite element.

reference_value_shape()[source]#: Return the shape of the value space on the reference cell.

reference_value_size()[source]#: Return the integer product of the reference value shape.

sub_elements()[source]#: Return list of sub-elements.

symmetry()[source]#: Return the symmetry dict, which is a mapping \(c_0 \to c_1\) meaning that component \(c_0\) is represented by component \(c_1\). A component is a tuple of one or more ints.

value_shape()[source]#: Return the shape of the value space on the global domain.

value_size()[source]#: Return the integer product of the value shape.

variant()[source]#: Return the variant used to initialise the element.

class cbcbeat.FitzHughNagumoManual(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

A reparametrized FitzHughNagumo model, based on Section 2.4.1 in “Computing the electrical activity in the heart” by Sundnes et al, 2006.

This is a model containing two nonlinear, ODEs for the evolution of the transmembrane potential v and one additional state variable s.

F(v, s, time=None)[source]#: Return right-hand side for state variable evolution.

I(v, s, time=None)[source]#: Return the ionic current.

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

num_states()[source]#: Return number of state variables.

class cbcbeat.Form(integrals)[source]#

Bases: object

Description of a weak form consisting of a sum of integrals over subdomains.

arguments()[source]#: Return all Argument objects found in form.

coefficient_numbering()[source]#: Return a contiguous numbering of coefficients in a mapping {coefficient:number}.

coefficients()[source]#: Return all Coefficient objects found in form.

constant_numbering()[source]#: Return a contiguous numbering of constants in a mapping {constant:number}.

constants()[source]#

domain_numbering()[source]#: Return a contiguous numbering of domains in a mapping {domain:number}.

empty()[source]#: Returns whether the form has no integrals.

equals(other)[source]#: Evaluate bool(lhs_form == rhs_form).

geometric_dimension()[source]#: Return the geometric dimension shared by all domains and functions in this form.

integrals()[source]#: Return a sequence of all integrals in form.

integrals_by_domain(domain)[source]#: Return a sequence of all integrals with a particular integration domain.

integrals_by_type(integral_type)[source]#: Return a sequence of all integrals with a particular domain type.

max_subdomain_ids()[source]#: Returns a mapping on the form {domain:{integral_type:max_subdomain_id}}.

signature()[source]#: Signature for use with jit cache (independent of incidental numbering of indices etc.)

subdomain_data()[source]#: Returns a mapping on the form {domain:{integral_type: subdomain_data}}.

ufl_cell()[source]#: Return the single cell this form is defined on, fails if multiple cells are found.

ufl_domain()[source]#

Return the single geometric integration domain occuring in the form.

Fails if multiple domains are found.

NB! This does not include domains of coefficients defined on other meshes, look at form data for that additional information.

ufl_domains()[source]#

Return the geometric integration domains occuring in the form.

NB! This does not include domains of coefficients defined on other meshes.

The return type is a tuple even if only a single domain exists.

class cbcbeat.FunctionSpace(domain, element)[source]#

Bases: AbstractFunctionSpace

ufl_domain()[source]#: Return ufl domain.

ufl_domains()[source]#: Return ufl domains.

ufl_element()[source]#: Return ufl element.

ufl_sub_spaces()[source]#: Return ufl sub spaces.

class cbcbeat.HCurlElement(element)[source]#

Bases: FiniteElementBase

A curl-conforming version of an outer product element, assuming this makes mathematical sense.

mapping()[source]#: Not implemented.

reconstruct(**kwargs)[source]#

shortstr()[source]#: Format as string for pretty printing.

sobolev_space()[source]#: Return the underlying Sobolev space.

class cbcbeat.HDivElement(element)[source]#

Bases: FiniteElementBase

A div-conforming version of an outer product element, assuming this makes mathematical sense.

mapping()[source]#: Not implemented.

reconstruct(**kwargs)[source]#

shortstr()[source]#: Format as string for pretty printing.

sobolev_space()[source]#: Return the underlying Sobolev space.

class cbcbeat.Identity(dim)[source]#

Bases: ConstantValue

UFL literal type: Representation of an identity matrix.

evaluate(x, mapping, component, index_values)[source]#: Evaluates the identity matrix on the given components.

ufl_shape#

class cbcbeat.Index(count=None)[source]#

Bases: IndexBase

UFL value: An index with no value assigned.

Used to represent free indices in Einstein indexing notation.

count()[source]#

class cbcbeat.Integral(integrand, integral_type, domain, subdomain_id, metadata, subdomain_data)[source]#

Bases: object

An integral over a single domain.

integral_type()[source]#: Return the domain type of this integral.

integrand()[source]#: Return the integrand expression, which is an Expr instance.

metadata()[source]#: Return the compiler metadata this integral has been annotated with.

reconstruct(integrand=None, integral_type=None, domain=None, subdomain_id=None, metadata=None, subdomain_data=None)[source]#: Construct a new Integral object with some properties replaced with new values.

1. Example:#

<a = Integral instance> b = a.reconstruct(expand_compounds(a.integrand())) c = a.reconstruct(metadata={‘quadrature_degree’:2})

subdomain_data()[source]#: Return the domain data of this integral.

subdomain_id()[source]#: Return the subdomain id of this integral.

ufl_domain()[source]#: Return the integration domain of this integral.

cbcbeat.InteriorElement(element)[source]#: Constructs the restriction of a finite element to the interior of the cell.

class cbcbeat.Jacobian(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The Jacobian of the mapping from reference cell to spatial coordinates.

\[J_{ij} = \frac{dx_i}{dX_j}\]

is_cellwise_constant()[source]#: Return whether this expression is spatially constant over each cell.

name = 'J'#

property ufl_shape#: Return the number of coordinates defined (i.e. the geometric dimension of the domain).

class cbcbeat.JacobianDeterminant(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The determinant of the Jacobian.

Represents the signed determinant of a square Jacobian or the pseudo-determinant of a non-square Jacobian.

is_cellwise_constant()[source]#: Return whether this expression is spatially constant over each cell.

name = 'detJ'#

class cbcbeat.JacobianInverse(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The inverse of the Jacobian.

Represents the inverse of a square Jacobian or the pseudo-inverse of a non-square Jacobian.

is_cellwise_constant()[source]#: Return whether this expression is spatially constant over each cell.

name = 'K'#

property ufl_shape#: Return the number of coordinates defined (i.e. the geometric dimension of the domain).

class cbcbeat.Logger(name, exception_type=<class 'Exception'>)[source]#

Bases: object

add_indent(increment=1)[source]#: Add to indentation level.

add_logfile(filename=None, mode='a', level=10)[source]#: Add a log file.

begin(*message)[source]#: Begin task: write message and increase indentation level.

debug(*message)[source]#: Write debug message.

deprecate(*message)[source]#: Write deprecation message.

end()[source]#: End task: write a newline and decrease indentation level.

error(*message)[source]#: Write error message and raise an exception.

get_handler()[source]#: Get handler for logging.

get_logfile_handler(filename)[source]#: Gets the handler to the file identified by the given file name.

get_logger()[source]#: Return message logger.

info(*message)[source]#: Write info message.

info_blue(*message)[source]#: Write info message in blue.

info_green(*message)[source]#: Write info message in green.

info_red(*message)[source]#: Write info message in red.

log(level, *message)[source]#: Write a log message on given log level.

pop_level()[source]#: Pop log level from the level stack, reverting to before the last push_level.

push_level(level)[source]#: Push a log level on the level stack.

set_handler(handler)[source]#: Replace handler for logging. To add additional handlers instead of replacing the existing one, use log.get_logger().addHandler(myhandler). See the logging module for more details.

set_indent(level)[source]#: Set indentation level.

set_level(level)[source]#: Set log level.

set_prefix(prefix)[source]#: Set prefix for log messages.

warning(*message)[source]#: Write warning message.

warning_blue(*message)[source]#: Write warning message in blue.

warning_green(*message)[source]#: Write warning message in green.

warning_red(*message)[source]#: Write warning message in red.

class cbcbeat.Markerwise(objects, keys, markers)[source]#

Bases: object

A container class representing an object defined by a number of objects combined with a mesh function defining mesh domains and a map between the these.

Arguments

objects (tuple): the different objects
keys (tuple of ints): a map from the objects to the domains marked in markers
markers (dolfin.MeshFunction): a mesh function mapping which domains the mesh consist of

Example of usage

Given (g0, g1), (2, 5) and markers, let

g = g0 on domains marked by 2 in markers g = g1 on domains marked by 5 in markers

letting:

g = Markerwise((g0, g1), (2, 5), markers)

keys()[source]#: The keys or domain numbers

markers()[source]#: The markers

values()[source]#: The objects

cbcbeat.Max(x, y)[source]#: UFL operator: Take the maximum of x and y.

class cbcbeat.MaxCellEdgeLength(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The maximum edge length of the cell.

name = 'maxcelledgelength'#

class cbcbeat.MaxFacetEdgeLength(domain)[source]#

Bases: GeometricFacetQuantity

UFL geometry representation: The maximum edge length of the facet.

name = 'maxfacetedgelength'#

class cbcbeat.Measure(integral_type, domain=None, subdomain_id='everywhere', metadata=None, subdomain_data=None)[source]#

Bases: object

integral_type()[source]#

Return the domain type.

Valid domain types are “cell”, “exterior_facet”, “interior_facet”, etc.

metadata()[source]#: Return the integral metadata. This data is not interpreted by UFL. It is passed to the form compiler which can ignore it or use it to compile each integral of a form in a different way.

reconstruct(integral_type=None, subdomain_id=None, domain=None, metadata=None, subdomain_data=None)[source]#

Construct a new Measure object with some properties replaced with new values.

1. Example:#

<dm = Measure instance> b = dm.reconstruct(subdomain_id=2) c = dm.reconstruct(metadata={ “quadrature_degree”: 3 })

Used by the call operator, so this is equivalent:: b = dm(2) c = dm(0, { “quadrature_degree”: 3 })

subdomain_data()[source]#: Return the integral subdomain_data. This data is not interpreted by UFL. Its intension is to give a context in which the domain id is interpreted.

subdomain_id()[source]#: Return the domain id of this measure (integer).

ufl_domain()[source]#

Return the domain associated with this measure.

This may be None or a Domain object.

class cbcbeat.Mesh(coordinate_element, ufl_id=None, cargo=None)[source]#

Bases: AbstractDomain

Symbolic representation of a mesh.

is_piecewise_linear_simplex_domain()[source]#

ufl_cargo()[source]#: Return carried object that will not be used by UFL.

ufl_cell()[source]#

ufl_coordinate_element()[source]#

ufl_id()#: Return the ufl_id of this object.

class cbcbeat.MeshView(mesh, topological_dimension, ufl_id=None)[source]#

Bases: AbstractDomain

Symbolic representation of a mesh.

is_piecewise_linear_simplex_domain()[source]#

ufl_cell()[source]#

ufl_id()#: Return the ufl_id of this object.

ufl_mesh()[source]#

cbcbeat.Min(x, y)[source]#: UFL operator: Take the minimum of x and y.

class cbcbeat.MinCellEdgeLength(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The minimum edge length of the cell.

name = 'mincelledgelength'#

class cbcbeat.MinFacetEdgeLength(domain)[source]#

Bases: GeometricFacetQuantity

UFL geometry representation: The minimum edge length of the facet.

name = 'minfacetedgelength'#

class cbcbeat.MixedElement(*elements, **kwargs)[source]#

Bases: FiniteElementBase

A finite element composed of a nested hierarchy of mixed or simple elements.

degree(component=None)[source]#: Return polynomial degree of finite element.

extract_component(i)[source]#: Recursively extract component index relative to a (simple) element and that element for given value component index.

extract_reference_component(i)[source]#: Recursively extract reference_component index relative to a (simple) element and that element for given value reference_component index.

extract_subelement_component(i)[source]#: Extract direct subelement index and subelement relative component index for a given component index.

extract_subelement_reference_component(i)[source]#: Extract direct subelement index and subelement relative reference_component index for a given reference_component index.

is_cellwise_constant(component=None)[source]#: Return whether the basis functions of this element is spatially constant over each cell.

mapping()[source]#: Not implemented.

num_sub_elements()[source]#: Return number of sub elements.

reconstruct(**kwargs)[source]#

reconstruct_from_elements(*elements)[source]#: Reconstruct a mixed element from new subelements.

shortstr()[source]#: Format as string for pretty printing.

sub_elements()[source]#: Return list of sub elements.

symmetry()[source]#: Return the symmetry dict, which is a mapping \(c_0 \to c_1\) meaning that component \(c_0\) is represented by component \(c_1\). A component is a tuple of one or more ints.

class cbcbeat.MixedFunctionSpace(*args)[source]#

Bases: AbstractFunctionSpace

num_sub_spaces()[source]#

ufl_domain()[source]#: Return ufl domain.

ufl_domains()[source]#: Return ufl domains.

ufl_element()[source]#

ufl_elements()[source]#: Return ufl elements.

ufl_sub_space(i)[source]#: Return i-th ufl sub space.

ufl_sub_spaces()[source]#: Return ufl sub spaces.

class cbcbeat.MonodomainSolver(mesh, time, M_i, I_s=None, v_=None, params=None)[source]#

Bases: BasicMonodomainSolver

This solver is based on a theta-scheme discretization in time and CG_1 elements in space.

Note

Arguments

mesh (dolfin.Mesh): The spatial domain (mesh)
time (dolfin.Constant or None): A constant holding the current time. If None is given, time is created for you, initialized to zero.
M_i (ufl.Expr): The intracellular conductivity tensor (as an UFL expression)
I_s (dict, optional): A typically time-dependent external stimulus given as a dict, with domain markers as the key and a dolfin.Expression as values. NB: it is assumed that the time dependence of I_s is encoded via the ‘time’ Constant.
v_ (ufl.Expr, optional): Initial condition for v. A new dolfin.Function will be created if none is given.
params (dolfin.Parameters, optional): Solver parameters

static default_parameters()[source]#

Initialize and return a set of default parameters

Returns: A set of parameters (dolfin.Parameters)

To inspect all the default parameters, do:

info(MonodomainSolver.default_parameters(), True)

property linear_solver#: The linear solver (dolfin.LUSolver or dolfin.KrylovSolver).

step(interval)[source]#

Solve on the given time step (t0, t1).

Arguments

interval (tuple): The time interval (t0, t1) for the step

Invariants

Assuming that v_ is in the correct state for t0, gives self.v in correct state at t1.

variational_forms(k_n)[source]#

Create the variational forms corresponding to the given discretization of the given system of equations.

Arguments

k_n (ufl.Expr or float): The time step

Returns

(lhs, rhs, prec) (tuple of ufl.Form)

class cbcbeat.MultiCellModel(models, keys, markers)[source]#

Bases: CardiacCellModel

F(v, s, time=None, index=None)[source]#: Return right-hand side for state variable evolution.

I(v, s, time=None, index=None)[source]#: Return the ionic current.

initial_conditions()[source]#: Return initial conditions for v and s as a dolfin.dolfin.cpp.function.GenericFunction.

keys()[source]#

markers()[source]#

mesh()[source]#

models()[source]#

num_models()[source]#

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

class cbcbeat.NoCellModel(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

Class representing no cell model (only bidomain equations). It actually just represents a single completely decoupled ODE.

F(v, s, time=None)[source]#: Return right-hand side for state variable evolution.

I(v, s, time=None)[source]#: Return the ionic current.

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

class cbcbeat.NodalEnrichedElement(*elements)[source]#

Bases: EnrichedElementBase

The vector sum of several finite element spaces:

\[\textrm{EnrichedElement}(V, Q) = \{v + q | v \in V, q \in Q\}.\]

Primal basis is reorthogonalized to dual basis which is a concatenation of subelements dual bases; resulting element is nodal.

is_cellwise_constant()[source]#: Return whether the basis functions of this element is spatially constant over each cell.

shortstr()[source]#: Format as string for pretty printing.

cbcbeat.Not(condition)[source]#: UFL operator: A boolean expression (not condition) for use with conditional.

cbcbeat.Or(left, right)[source]#: UFL operator: A boolean expression (left or right) for use with conditional.

class cbcbeat.PermutationSymbol(dim)[source]#

Bases: ConstantValue

UFL literal type: Representation of a permutation symbol.

This is also known as the Levi-Civita symbol, antisymmetric symbol, or alternating symbol.

evaluate(x, mapping, component, index_values)[source]#: Evaluates the permutation symbol.

ufl_shape#

class cbcbeat.RestrictedElement(element, restriction_domain)[source]#

Bases: FiniteElementBase

Represents the restriction of a finite element to a type of cell entity.

is_cellwise_constant()[source]#: Return whether the basis functions of this element is spatially constant over each cell.

mapping()[source]#: Not implemented.

num_restricted_sub_elements()[source]#: Return number of restricted sub elements.

num_sub_elements()[source]#: Return number of sub elements.

reconstruct(**kwargs)[source]#

restricted_sub_elements()[source]#: Return list of restricted sub elements.

restriction_domain()[source]#: Return the domain onto which the element is restricted.

shortstr()[source]#: Format as string for pretty printing.

sub_element()[source]#: Return the element which is restricted.

sub_elements()[source]#: Return list of sub elements.

symmetry()[source]#: Return the symmetry dict, which is a mapping \(c_0 \to c_1\) meaning that component \(c_0\) is represented by component \(c_1\). A component is a tuple of one or more ints.

class cbcbeat.RogersMcCulloch(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

The Rogers-McCulloch model is a modified FitzHughNagumo model. This

formulation follows the description on page 2 of “Optimal control approach …” by Nagaiah, Kunisch and Plank, 2013, J Math Biol with w replaced by s. Note that this model introduces one additional parameter compared to the original 1994 Rogers-McCulloch model.

This is a model containing two nonlinear, ODEs for the evolution of the transmembrane potential v and one additional state variable s:

\[\]

rac{dv}{dt} = - I_{ion}(v, s)

rac{ds}{dt} = F(v, s)

where

\[ \begin{align}\begin{aligned}I_{ion}(v, s) = g v (1 - v/v_th)(1 - v/v_p) + \eta_1 v s\\ F(v, s) = \eta_2 (v/vp - \eta_3 s)\end{aligned}\end{align} \]

F(v, s, time=None)[source]#: Return right-hand side for state variable evolution.

I(v, s, time=None)[source]#: Return the ionic current.

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

num_states()[source]#: Return number of state variables.

class cbcbeat.SingleCellSolver(model, time, params=None)[source]#: Bases: CardiacODESolver

class cbcbeat.SpatialCoordinate(domain)[source]#

Bases: GeometricCellQuantity

UFL geometry representation: The coordinate in a domain.

In the context of expression integration, represents the domain coordinate of each quadrature point.

In the context of expression evaluation in a point, represents the value of that point.

count()[source]#

evaluate(x, mapping, component, index_values)[source]#: Return the value of the coordinate.

is_cellwise_constant()[source]#: Return whether this expression is spatially constant over each cell.

name = 'x'#

property ufl_shape#: Return the number of coordinates defined (i.e. the geometric dimension of the domain).

class cbcbeat.SplittingSolver(model, params=None, V_index=0)[source]#

Bases: BasicSplittingSolver

An optimised solver for the bidomain equations based on the operator splitting scheme described in Sundnes et al 2006, p. 78 ff.

The solver computes as solutions:

“vs” (dolfin.Function) representing the solution for the transmembrane potential and any additional state variables, and

“vur” (dolfin.Function) representing the transmembrane potential in combination with the extracellular potential and an additional Lagrange multiplier.

Arguments

model (cbcbeat.cardiacmodels.CardiacModel): a CardiacModel object describing the simulation set-up
params (dolfin.Parameters, optional): a Parameters object controlling solver parameters

Example of usage:

from cbcbeat import *

# Describe the cardiac model
mesh = UnitSquareMesh(100, 100)
time = Constant(0.0)
cell_model = FitzHughNagumoManual()
stimulus = Expression("10*t*x[0]", t=time, degree=1)
cm = CardiacModel(mesh, time, 1.0, 1.0, cell_model, stimulus)

# Extract default solver parameters
ps = SplittingSolver.default_parameters()

# Examine the default parameters
info(ps, True)

# Update parameter dictionary
ps["theta"] = 1.0 # Use first order splitting
ps["CardiacODESolver"]["scheme"] = "GRL1" # Use Generalized Rush-Larsen scheme

# Use monodomain equations as the PDE model
ps["pde_solver"] = "monodomain"
# Use direct linear solver of the PDEs
ps["MonodomainSolver"]["linear_solver_type"] = "direct"
# Use backward Euler for temporal discretization for the PDEs
ps["MonodomainSolver"]["theta"] = 1.0

solver = SplittingSolver(cm, params=ps)

# Extract the solution fields and set the initial conditions
(vs_, vs, vur) = solver.solution_fields()
vs_.assign(cell_model.initial_conditions())

# Solve
dt = 0.1
T = 1.0
interval = (0.0, T)
for (timestep, fields) in solver.solve(interval, dt):
    (vs_, vs, vur) = fields
    # Do something with the solutions

Assumptions

The cardiac conductivities do not vary in time

static default_parameters()[source]#

Initialize and return a set of default parameters for the splitting solver

Returns: The set of default parameters (dolfin.Parameters)

Example of usage:

info(SplittingSolver.default_parameters(), True)

cbcbeat.TensorConstant(domain, count=None)[source]#

class cbcbeat.TensorElement(family, cell=None, degree=None, shape=None, symmetry=None, quad_scheme=None, variant=None)[source]#

Bases: MixedElement

A special case of a mixed finite element where all elements are equal.

extract_subelement_component(i)[source]#: Extract direct subelement index and subelement relative component index for a given component index.

flattened_sub_element_mapping()[source]#

mapping()[source]#: Not implemented.

reconstruct(**kwargs)[source]#

shortstr()[source]#: Format as string for pretty printing.

symmetry()[source]#: Return the symmetry dict, which is a mapping \(c_0 \to c_1\) meaning that component \(c_0\) is represented by component \(c_1\). A component is a tuple of one or more ints.

variant()[source]#: Return the variant used to initialise the element.

class cbcbeat.TensorProductCell(*cells, **kwargs)[source]#

Bases: AbstractCell

cellname()[source]#: Return the cellname of the cell.

has_simplex_facets()[source]#: Return True if all the facets of this cell are simplex cells.

is_simplex()[source]#: Return True if this is a simplex cell.

num_edges()[source]#: The number of cell edges.

num_facets()[source]#: The number of cell facets.

num_vertices()[source]#: The number of cell vertices.

reconstruct(geometric_dimension=None)[source]#

sub_cells()[source]#: Return list of cell factors.

class cbcbeat.TensorProductElement(*elements, **kwargs)[source]#

Bases: FiniteElementBase

The tensor product of \(d\) element spaces:

\[V = V_1 \otimes V_2 \otimes ... \otimes V_d\]

Given bases \(\{\phi_{j_i}\}\) of the spaces \(V_i\) for \(i = 1, ...., d\), \(\{ \phi_{j_1} \otimes \phi_{j_2} \otimes \cdots \otimes \phi_{j_d} \}\) forms a basis for \(V\).

mapping()[source]#: Not implemented.

num_sub_elements()[source]#: Return number of subelements.

reconstruct(**kwargs)[source]#

shortstr()[source]#: Short pretty-print.

sobolev_space()[source]#: Return the underlying Sobolev space of the TensorProductElement.

sub_elements()[source]#: Return subelements (factors).

class cbcbeat.TensorProductMesh(meshes, ufl_id=None)[source]#

Bases: AbstractDomain

Symbolic representation of a mesh.

is_piecewise_linear_simplex_domain()[source]#

ufl_cell()[source]#

ufl_coordinate_element()[source]#

ufl_id()#: Return the ufl_id of this object.

ufl_meshes()[source]#

class cbcbeat.Tentusscher_2004_mcell(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

NOT_IMPLEMENTED

F(v, s, time=None)[source]#: Right hand side for ODE system

I(v, s, time=None)[source]#: Transmembrane current

I = -dV/dt

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

class cbcbeat.Tentusscher_panfilov_2006_M_cell(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

F(v, s, time=None)[source]#: Right hand side for ODE system

I(v, s, time=None)[source]#: Transmembrane current

I = -dV/dt

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

class cbcbeat.Tentusscher_panfilov_2006_epi_cell(params=None, init_conditions=None)[source]#

Bases: CardiacCellModel

F(v, s, time=None)[source]#: Right hand side for ODE system

I(v, s, time=None)[source]#: Transmembrane current

I = -dV/dt

static default_initial_conditions()[source]#: Set-up and return default initial conditions.

static default_parameters()[source]#: Set-up and return default parameters.

num_states()[source]#: Return number of state variables (in addition to the membrane potential).

cbcbeat.TestFunction(function_space, part=None)[source]#: UFL value: Create a test function argument to a form.

cbcbeat.TestFunctions(function_space)[source]#: UFL value: Create a TestFunction in a mixed space, and return a tuple with the function components corresponding to the subelements.

cbcbeat.TrialFunction(function_space, part=None)[source]#: UFL value: Create a trial function argument to a form.

cbcbeat.TrialFunctions(function_space)[source]#: UFL value: Create a TrialFunction in a mixed space, and return a tuple with the function components corresponding to the subelements.

exception cbcbeat.UFLException[source]#

Bases: Exception

Base class for UFL exceptions.

cbcbeat.VectorConstant(domain, count=None)[source]#

class cbcbeat.VectorElement(family, cell=None, degree=None, dim=None, form_degree=None, quad_scheme=None, variant=None)[source]#

Bases: MixedElement

A special case of a mixed finite element where all elements are equal.

mapping()[source]#: Not implemented.

reconstruct(**kwargs)[source]#

shortstr()[source]#: Format as string for pretty printing.

variant()[source]#: Return the variant used to initialise the element.

class cbcbeat.WithMapping(wrapee, mapping)[source]#

Bases: FiniteElementBase

Specify an alternative mapping for the wrappee. For example, to use identity mapping instead of Piola map with an element E, write remapped = WithMapping(E, “identity”)

mapping()[source]#: Not implemented.

reconstruct(**kwargs)[source]#

reference_value_shape()[source]#: Return the shape of the value space on the reference cell.

shortstr()[source]#

value_shape()[source]#: Return the shape of the value space on the global domain.

cbcbeat.acos(f)[source]#: UFL operator: Take the inverse cosine of f.

cbcbeat.action(form, coefficient=None)[source]#: UFL form operator: Given a bilinear form, return a linear form with an additional coefficient, representing the action of the form on the coefficient. This can be used for matrix-free methods.

cbcbeat.add_indent(increment=1)#: Add to indentation level.

cbcbeat.add_logfile(filename=None, mode='a', level=10)#: Add a log file.

cbcbeat.adjoint(form, reordered_arguments=None)[source]#

UFL form operator: Given a combined bilinear form, compute the adjoint form by changing the ordering (count) of the test and trial functions, and taking the complex conjugate of the result.

By default, new Argument objects will be created with opposite ordering. However, if the adjoint form is to be added to other forms later, their arguments must match. In that case, the user must provide a tuple *reordered_arguments*=(u2,v2).

cbcbeat.as_cell(cell)[source]#

Convert any valid object to a Cell or return cell if it is already a Cell.

Allows an already valid cell, a known cellname string, or a tuple of cells for a product cell.

cbcbeat.as_domain(domain)[source]#: Convert any valid object to an AbstractDomain type.

cbcbeat.as_matrix(expressions, indices=None)[source]#: UFL operator: As as_tensor(), but limited to rank 2 tensors.

cbcbeat.as_tensor(expressions, indices=None)[source]#

UFL operator: Make a tensor valued expression.

This works in two different ways, by using indices or lists.

1) Returns \(A\) such that \(A\) [indices] = expressions. If indices are provided, expressions must be a scalar valued expression with all the provided indices among its free indices. This operator will then map each of these indices to a tensor axis, thereby making a tensor valued expression from a scalar valued expression with free indices.

2) Returns \(A\) such that \(A[k,...]\) = expressions*[k]. If no indices are provided, *expressions must be a list or tuple of expressions. The expressions can also consist of recursively nested lists to build higher rank tensors.

cbcbeat.as_ufl(expression)[source]#: Converts expression to an Expr if possible.

cbcbeat.as_vector(expressions, index=None)[source]#: UFL operator: As as_tensor(), but limited to rank 1 tensors.

cbcbeat.asin(f)[source]#: UFL operator: Take the inverse sine of f.

cbcbeat.atan(f)[source]#: UFL operator: Take the inverse tangent of f.

cbcbeat.atan_2(f1, f2)[source]#: UFL operator: Take the inverse tangent with two the arguments f1 and f2.

cbcbeat.avg(v)[source]#: UFL operator: Take the average of v across a facet.

cbcbeat.begin(*message)#: Begin task: write message and increase indentation level.

cbcbeat.bessel_I(nu, f)[source]#: UFL operator: regular modified cylindrical Bessel function.

cbcbeat.bessel_J(nu, f)[source]#: UFL operator: cylindrical Bessel function of the first kind.

cbcbeat.bessel_K(nu, f)[source]#: UFL operator: irregular modified cylindrical Bessel function.

cbcbeat.bessel_Y(nu, f)[source]#: UFL operator: cylindrical Bessel function of the second kind.

cbcbeat.cell_avg(f)[source]#: UFL operator: Take the average of v over a cell.

cbcbeat.cofac(A)[source]#: UFL operator: Take the cofactor of A.

cbcbeat.conditional(condition, true_value, false_value)[source]#: UFL operator: A conditional expression, taking the value of true_value when condition evaluates to true and false_value otherwise.

cbcbeat.conj(f)[source]#: UFL operator: The complex conjugate of f

cbcbeat.cos(f)[source]#: UFL operator: Take the cosine of f.

cbcbeat.cosh(f)[source]#: UFL operator: Take the hyperbolic cosine of f.

cbcbeat.cross(a, b)[source]#: UFL operator: Take the cross product of a and b.

cbcbeat.curl(f)[source]#: UFL operator: Take the curl of f.

cbcbeat.debug(*message)#: Write debug message.

cbcbeat.deprecate(*message)#: Write deprecation message.

cbcbeat.derivative(form, coefficient, argument=None, coefficient_derivatives=None)[source]#

UFL form operator: Compute the Gateaux derivative of form w.r.t. coefficient in direction of argument.

If the argument is omitted, a new Argument is created in the same space as the coefficient, with argument number one higher than the highest one in the form.

The resulting form has one additional Argument in the same finite element space as the coefficient.

A tuple of Coefficient s may be provided in place of a single Coefficient, in which case the new Argument argument is based on a MixedElement created from this tuple.

An indexed Coefficient from a mixed space may be provided, in which case the argument should be in the corresponding subspace of the coefficient space.

If provided, coefficient_derivatives should be a mapping from Coefficient instances to their derivatives w.r.t. coefficient.

cbcbeat.det(A)[source]#: UFL operator: Take the determinant of A.

cbcbeat.dev(A)[source]#: UFL operator: Take the deviatoric part of A.

cbcbeat.diag(A)[source]#

UFL operator: Take the diagonal part of rank 2 tensor A or make a diagonal rank 2 tensor from a rank 1 tensor.

Always returns a rank 2 tensor. See also diag_vector.

cbcbeat.diag_vector(A)[source]#

UFL operator: Take the diagonal part of rank 2 tensor A and return as a vector.

cbcbeat package

Contents

1. cbcbeat package#

1.1. Subpackages#

1.2. Submodules#

1.3. cbcbeat.bidomainsolver module#

1.4. cbcbeat.cardiacmodels module#

1.5. cbcbeat.cellsolver module#

1.6. cbcbeat.dolfinimport module#

1.7. cbcbeat.gossplittingsolver module#

1.8. cbcbeat.gotran2cellmodel module#

1.9. cbcbeat.markerwisefield module#

1.10. cbcbeat.monodomainsolver module#

1.11. cbcbeat.splittingsolver module#

1.12. cbcbeat.utils module#

1.13. Module contents#

1. Example:#

1. Example:#

1. Example:#

1. Example:#

1. Example:#