Python API

Python API#

sumfem#

Finite element models written as summable components, on FEniCSx.

The generic framework has three layers, imported as subpackages:

sumfem.symbolic, the symbolic description, without mesh data: domains, fields, parameters, residuals and functionals, systems of equations, sequences of stages, time derivatives, boundary conditions, transfers and dimensional analysis.
sumfem.discrete, which binds them to meshes: discretizations, and the discrete systems and boundary conditions.
sumfem.solvers, which solves the discrete systems.

The names users need most are exported here too. Where a symbolic and a discrete class share a name, such as System or DirichletBC, the symbolic one is exported.

Battery physics and the other model families live in sumfem.models, and post-processing in sumfem.post_processing:

import sumfem

domain = sumfem.Domain("square", "triangle")
u = sumfem.Field("u", domain, ("Lagrange", 1))

The symbolic layer#

The symbolic description of variational problems, without mesh data.

Domain: a symbolic description of a type of mesh, with named markers.
Cartesian, Cylindrical, Spherical and Weighted: coordinate systems of domains, which weight their integrals, e.g. for axisymmetric or radially symmetric bodies.
Field, Parameter, KnownCoefficient: the named objects residuals and functionals are written with, with units.
Residual: a reusable contribution to a residual, summed with + into a ResidualSum.
Functional: a reusable contribution to a scalar functional, summed into a FunctionalSum, whose variation is a residual.
build_system(): turns a residual into blocked residual and Jacobian forms, optionally in dimensionless form.
Dt() and Clock: semi-discrete transient terms.
DirichletBC: boundary conditions by marker name.
check_units(), equation_scales() and pi_groups(): dimensional analysis with ufl_units.

class sumfem.symbolic.Cartesian[source]#

Integrals without a weight, the default of every domain.

The weight is the integer one, which UFL drops from products, and integrands are kept as they are, so forms keep their integrands.

check(gdim)[source]#

Any geometric dimension.

Parameters:: gdim (int)
Return type:: None

class sumfem.symbolic.CoordinateSystem(*args, **kwargs)[source]#

How the integrals over a domain are weighted.

A coordinate system may also write the vectors of integrands in its basis, with the method lift(integrand, x), see Cylindrical.lift(), and say whether it supports vector-valued unknowns with the class attribute vectors. Without them, integrands are kept as they are and only scalar unknowns are supported, see lift() and supports_vectors().

weight(x)[source]#

The weight of an integral.

Parameters:: x (Expr) – The spatial coordinate of the domain.
Return type:: Expr

check(gdim)[source]#

Check that a domain of geometric dimension gdim can use it.

Raises:: ComponentError – If it cannot.
Parameters:: gdim (int)
Return type:: None

class sumfem.symbolic.Cylindrical[source]#

An axisymmetric body, with the weight 2 pi r.

The mesh lies in the half-plane x[0] = r >= 0, with x[1] = z on a 2D mesh. On a 1D mesh in r the integrals are per unit length of an infinite cylinder. On the axis the weight vanishes, so the symmetry condition of scalars there is the natural one. A vector’s radial component vanishes on the axis, which is a Dirichlet condition.

Vectors are written in the basis (r, z, theta): a field of shape (gdim,) has no swirl, and one of shape (3,) has the swirl u_theta as its last component. Fields do not depend on theta.

lift(integrand, x)[source]#

See _Lift.

Parameters:

integrand (Expr)
x (Expr)

Return type:

Expr

gradient(f, x)[source]#

The gradient of a scalar or a vector in (r, z, theta).

Raises:

ComponentError – For other shapes.

Parameters:

f (Expr)
x (Expr)

Return type:

Expr

divergence(f, x)[source]#

The divergence of a vector or a tensor in (r, z, theta).

Raises:

ComponentError – For other shapes.

Parameters:

f (Expr)
x (Expr)

Return type:

Expr

class sumfem.symbolic.Spherical[source]#

A radially symmetric body, with the weight 4 pi r^2.

The mesh is an interval in x[0] = r >= 0. Vectors are written in the basis (r, theta, phi), and only have a radial component.

lift(integrand, x)[source]#

See _Lift.

Parameters:

integrand (Expr)
x (Expr)

Return type:

Expr

gradient(f, x)[source]#

The gradient of a scalar or a radial vector in (r, theta, phi).

Raises:

ComponentError – For other shapes.

Parameters:

f (Expr)
x (Expr)

Return type:

Expr

divergence(f, x)[source]#

The divergence of a radial vector or a diagonal tensor.

Raises:

ComponentError – For other shapes.

Parameters:

f (Expr)
x (Expr)

Return type:

Expr

class sumfem.symbolic.Weighted(weight_of)[source]#

Integrals with any weight of the spatial coordinate.

E.g. Weighted(lambda x: 2 * pi * (1 - cos(theta)) * x[0] ** 2) for a cone of half-angle theta around its axis, in the distance x[0] from its apex. Integrands are kept as they are, so only scalar unknowns are supported.

Parameters:: weight_of (Callable[[Expr], Expr | float])

weight_of#

The weight, as a function of the spatial coordinate.

Type:: Callable[[ufl.core.expr.Expr], ufl.core.expr.Expr | float]

check(gdim)[source]#

Any geometric dimension.

Parameters:: gdim (int)
Return type:: None

sumfem.symbolic.apply_coordinates(form)[source]#

A form with each integral in its domain’s coordinate system.

Each integrand is lifted, see CoordinateSystem.lift(), and weighted. Integrals over a mesh that is not a Domain are Cartesian. On interior facets the weight is taken from the "+" side.

Parameters:: form (Form | ZeroBaseForm)
Return type:: Form | ZeroBaseForm

exception sumfem.symbolic.ComponentError[source]#: Raised when a component, or a sum of components, is defined inconsistently.

exception sumfem.symbolic.SystemBuildError[source]#: Raised when a residual cannot be turned into a system of equations.

sumfem.symbolic.Dt(expression)[source]#

The time derivative of an expression of fields.

For a field this is field.dot. For any other expression, such as the momentum H * u of the shallow water equations, it is a TimeDerivative, which a time-stepping scheme replaces by a difference of the whole expression over time steps. That keeps conserved quantities conserved, which the chain rule u * Dt(H) + H * Dt(u) would not.

Parameters:: expression (Expr)
Return type:: Expr

class sumfem.symbolic.Clock(domain, t0=0.0, dt=None, name='clock')[source]#

The current time t and the time step dt as parameters.

Both are Parameter objects, so models can use them in expressions, for instance a time-dependent source ufl.sin(clock.t).

Parameters:

domain (Domain) – The domain the parameters are defined on.
t0 (float) – Initial time.
dt (float | None) – Default time step size.
name (str) – Prefix of the parameter names, "<name>.t" and "<name>.dt".

property t: Parameter#: Current time.

property dt: Parameter#: Time step size.

class sumfem.symbolic.Component(name)[source]#

What residuals and functionals share: a name, and named inputs.

Parameters:: name (str) – Unique name of the component. It prefixes the names of parameters and known coefficients the component creates, as "<name>.<parameter>".

property name: str#: Name of the component.

parameter(name, domain, shape=(), default=None, unit=None, doc='')[source]#

Return the parameter "<component>.<name>", creating it on first call.

Calling this again with the same name returns the same object, so it is safe to call while building a form.

Parameters:

name (str)
domain (Domain)
shape (tuple[int, ...])
default (Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None)
unit (Expr | None)
doc (str)

Return type:

Parameter

known_coefficient(name, domain, element, unit=None, scale=1.0, doc='')[source]#

Return the known coefficient "<component>.<name>", made on first call.

Parameters:

name (str)
domain (Domain)
element (AbstractFiniteElement | str | tuple)
unit (Expr | None)
scale (float)
doc (str)

Return type:

KnownCoefficient

input(name, value=None, *, domain, element=None, shape=(), default=None, unit=None, scale=1.0, doc='')[source]#

Resolve an input given to the constructor.

If value is given, it is converted with as_expr(). Passing a Field couples this component to the field. Otherwise a placeholder named "<component>.<name>" is created: a KnownCoefficient when element is given, else a Parameter of the given shape. unit applies to either, scale to a known coefficient and default to a parameter.

Parameters:

name (str)
value (Expr | Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None)
domain (Domain)
element (AbstractFiniteElement | str | tuple | None)
shape (tuple[int, ...])
default (Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None)
unit (Expr | None)
scale (float)
doc (str)

Return type:

Expr

property components: tuple[Component, ...]#: The summands. A single component is its own summand.

class sumfem.symbolic.Functional(name)[source]#

Base class for a reusable contribution to a scalar functional.

Subclasses override functional() with a form of rank zero, e.g. an energy, an objective or a quantity of interest. Functionals sum with + into a FunctionalSum. The form may be assembled, e.g. with assemble_scalar(), and to_residual() gives the residual of its stationary point.

Parameters:: name (str) – Unique name of the component. It prefixes the names of parameters and known coefficients the component creates.

functional()[source]#

Return this component’s contribution to the functional.

Return type:: Form

property components: tuple[Functional, ...]#: The summands. A single component is its own summand.

property form: Form[source]#: The functional, built on first access and then cached.

property inventory: Inventory[source]#: The named objects used by this component, validated.

property fields: tuple[Field, ...]#: Every field the functional uses, in creation order.

property parameters: Mapping[str, Parameter]#: Parameters used by the functional, by name.

property known_coefficients: Mapping[str, KnownCoefficient]#: Known coefficients used by the functional, by name.

property domains: tuple[Domain, ...]#: Domains used by the functional.

to_residual(fields=None, name=None)[source]#

The residual of the functional’s stationary point.

Its residual is sum_f derivative(J, f, f.v) over the given fields, so each field gets the equation of the functional’s variation in it.

Parameters:

fields (Sequence[Field] | None) – The fields to vary, e.g. the unknowns of one stage. Defaults to every field of the functional.
name (str | None) – Name of the residual. Defaults to the functional’s name.

Return type:

Residual

class sumfem.symbolic.FunctionalSum(components, name=None)[source]#

The sum of several functionals, flattened like a ResidualSum.

Parameters:

components (Sequence[_C]) – The components to sum.
name (str | None) – Name of the sum. Defaults to the summand names joined by +.

Raises:

ComponentError – If the same component appears twice.

functional()[source]#

Sum of the summands’ functionals.

Return type:: Form

class sumfem.symbolic.Jacobian(name)[source]#

Base class for a reusable bilinear form: a linearization or preconditioner.

Subclasses override jacobian() with a form in the test functions field.v and trial functions field.du. Jacobians sum with + into a JacobianSum, and a residual takes one as its Jacobian or preconditioner with Residual.with_jacobian() or Residual.with_preconditioner(), or from its hooks.

Parameters:: name (str) – Unique name of the component. It prefixes the names of parameters and known coefficients the component creates.

jacobian()[source]#

Return this component’s bilinear form.

Return type:: Form | ZeroBaseForm

property components: tuple[Jacobian, ...]#: The summands. A single component is its own summand.

property form: Form | ZeroBaseForm[source]#: The bilinear form, built on first access and then cached.

property inventory: Inventory[source]#: The named objects used by this component, validated.

property fields: tuple[Field, ...]#: Every field the form uses, in creation order.

property parameters: Mapping[str, Parameter]#: Parameters used by the form, by name.

property known_coefficients: Mapping[str, KnownCoefficient]#: Known coefficients used by the form, by name.

property domains: tuple[Domain, ...]#: Domains used by the form.

linearized(unknowns)[source]#

The form in a system with these unknowns.

Derived linearizations, of Residual.to_jacobian() without fields, are taken with respect to the unknowns.

Parameters:: unknowns (Sequence[Field])
Return type:: Form | ZeroBaseForm

replaced(mapping)[source]#

This Jacobian with terminals replaced, derivatives left unevaluated.

See replace().

Parameters:: mapping (Mapping[Any, Expr])
Return type:: Jacobian

class sumfem.symbolic.JacobianSum(components, name=None)[source]#

The sum of several Jacobians, flattened like a ResidualSum.

Parameters:

components (Sequence[_C]) – The components to sum.
name (str | None) – Name of the sum. Defaults to the summand names joined by +.

Raises:

ComponentError – If the same component appears twice.

jacobian()[source]#

Sum of the summands’ forms.

Return type:: Form | ZeroBaseForm

class sumfem.symbolic.ResidualSum(components, name=None)[source]#

The sum of several residuals.

Sums are flattened, so (a + b) + c and a + (b + c) both have the summands (a, b, c). Use a + b or sum([a, b, c]) rather than constructing this class directly.

Parameters:

components (Sequence[_C]) – The components to sum.
name (str | None) – Name of the sum. Defaults to the summand names joined by +.

Raises:

ComponentError – If the same component appears twice.

residual()[source]#

Sum of the summands’ residuals.

Return type:: Form | ZeroBaseForm

class sumfem.symbolic.DirichletBC(field, marker, value, entity_type='facet')[source]#

A Dirichlet condition on a field, on the entities of a named marker.

The condition is symbolic: a Discretization turns it into a {py:class}`dolfinx.fem.DirichletBC` on the mesh of the field’s domain.

Parameters:

field (Field)
marker (str | tuple[str, ...])
value (Value | ufl.core.expr.Expr)
entity_type (EntityType)

field#

The field the condition applies to.

Type:: Field

marker#

Name of the marker on the field’s mesh, e.g. "Inlet", or several names.

Type:: str | tuple[str, …]

value#

The prescribed value. A UFL expression written with parameters, fields and the coordinates of the domain is evaluated at runtime, before every solve and every Newton iteration, so it follows time and other changing values. It is evaluated in SI base units. Any other Value is evaluated once, in the field’s unit.

Type:: Value | ufl.core.expr.Expr

entity_type#

The entities the marker refers to, facets by default.

Type:: EntityType

property markers: tuple[str, ...]#: The marker names, as a tuple.

class sumfem.symbolic.Subdomain(name, parent, dim, marker, length_scale=None)[source]#

Marked entities of a domain, as a domain of their own.

Fields on a subdomain live on a submesh of a subset of the parent’s mesh entities, those with the subdomain’s marker. These are cells, e.g. the electrolyte of a cell, or facets, e.g. the concentration of a species on a surface, or the unknown of a PDE on a manifold. dx integrates over the subdomain, ds over its boundary and dS over its interior facets. Markers are named after the parent’s: those of dx are parent markers of the subdomain’s dimension, and those of ds and dS parent markers of one dimension lower. For a surface of a 3D domain, these are facet and ridge markers. Its coordinates are the parent’s, and so is its coordinate system.

Whether a form is integrated on the submesh or on the parent is decided per form, see place_subdomain_integrals().

Parameters:

name (str) – Name of the domain.
parent (Domain) – The domain the entities belong to.
dim (int) – Topological dimension of the entities, e.g. parent.tdim - 1 for facets.
marker (str) – The parent’s marker of the entities, e.g. "Membrane".
length_scale (float | None) – See Domain. Defaults to the parent’s.

property parent: Domain#: The domain the entities belong to.

property marker: str#: The parent’s marker of the entities.

property codim: int#: Codimension of the subdomain in its parent.

property coordinates: CoordinateSystem#: The parent’s coordinate system, since the coordinates are the same.

measure(integral_type, *markers, degree=None, metadata=None)[source]#

A measure of this domain, as Domain.measure().

The markers are the parent’s markers of the measure’s entities, which have the same dimension in the parent: e.g. for a surface of a 3D domain, facet markers for dx, ridge markers for ds and dS, and vertex markers for dr.

Parameters:

integral_type (Literal['dx', 'ds', 'dS', 'dr', 'dP'])
markers (str)
degree (int | None)
metadata (dict | None)

Return type:

Measure

dx(*markers, degree=None, metadata=None)#

Cell measure of this domain, optionally restricted to named parts.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

ds(*markers, degree=None, metadata=None)#

Measure of this domain’s boundary, optionally of named parts of it.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

dS(*markers, degree=None, metadata=None)#

Measure of this domain’s interior facets, optionally of named ones.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

dr(*markers, degree=None, metadata=None)#

Ridge measure of this domain, optionally restricted to named ridges.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

dP(*markers, degree=None, metadata=None)#

Vertex measure of this domain, optionally restricted to named vertices.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

on_parent(integral)[source]#

An integral of this domain, as one on the parent’s entities.

A cell integral moves to the parent’s cells for a subdomain of cells, and to its exterior facets for a subdomain of facets. For a subdomain of facets between two subdomains of cells, an interface of the parent, see Domain.interface(), an integral over the whole subdomain moves to the interior facets of the interface, on its "+" side, which holds the same values of the facets’ functions as the other. An exterior facet integral of a subdomain of cells moves to the parent’s facets seen from the subdomain’s side, under a one-sided marker, see Domain.one_sided(): e.g. the joint of a branching surface, seen from one branch. The functions and the geometry of the submesh in it, e.g. its spatial coordinate, are mapped by the entity map, and its facet normal is the parent’s.

Raises:: ComponentError – For other integrals, or subdomains of codimension above one.
Parameters:: integral (Integral)
Return type:: Integral

class sumfem.symbolic.Sequential(stages, max_iterations=1, tol=None)[source]#

Solve stages one after another, each with the latest values of the others.

A model is a stage that solves for every field it writes an equation for, and a System one built with options of its own, e.g. a time-stepping scheme’s substitutions. The other fields a stage uses are data, taken from the stages that solve for them. A Transfer is a stage too, which writes data from one field to another, e.g. a velocity solved on one mesh into a known coefficient of another.

By default the sequence runs once, which is weak coupling. With a tolerance it is repeated until the solution stops changing, which is strong coupling.

A stage may itself be a sequence, which runs to its own tolerance whenever the outer sequence reaches it, e.g. the iterations of an enthalpy method inside a pass that then solves for the flow.

Parameters:

stages (Sequence[Stage | Sequential]) – The models, systems, transfers and nested sequences, in the order they run.
max_iterations (int) – Number of passes through the stages without a tolerance, or the maximum number with one.
tol (float | None) – Relative change of the unknowns of its stages between two passes below which strong coupling has converged.

Raises:

ValueError – For an empty sequence, a non-positive number of iterations, or a non-positive tolerance.
TypeError – If a stage is not a model, a system, a transfer or a sequence.

property stages: tuple[Residual | System | Transfer | Sequential, ...]#: The stages and nested sequences, in the order they run.

property max_iterations: int#: Number of passes, or maximum number of passes with a tolerance.

property tol: float | None#: Tolerance of the relative change for strong coupling, if any.

leaves()[source]#

The stages of this and the nested sequences, in order.

Return type:: tuple[Residual | System | Transfer, …]

map(function)[source]#

The same sequence, with each stage replaced by function(stage).

Parameters:: function (Callable[[Residual | System | Transfer], Residual | System | Transfer])
Return type:: Sequential

description()[source]#

The structure of the sequence by stage names, e.g. for metadata.

Return type:: dict[str, Any]

class sumfem.symbolic.Transfer(name, source, target)[source]#

A stage that writes a field or known coefficient from another one.

Transfers couple by data, not by summing forms: for instance between non-matching meshes, where no form can hold both, see Interpolation. A transfer is a stage of a Sequential sequence, or is applied by hand: disc[transfer].apply(). It only describes what is transferred. How is up to its discrete counterpart: register one for a new kind of transfer, e.g. a projection, with @counterpart.register.

Parameters:

name (str) – Name of the transfer, unique among the stages it runs with.
source (Transferable) – The field or known coefficient to read.
target (Transferable) – The known coefficient or field to write, usually on another domain. A stage may solve for it, starting from the transferred values, e.g. a correction or an initial guess.

Raises:

ComponentError – If the source and the target differ in shape or in the dimension of their units, or are the same.

property name: str#: Name of the transfer.

property source: Field | KnownCoefficient#: The field or known coefficient to read.

property target: Field | KnownCoefficient#: The known coefficient or field to write.

property target_domain: Domain#: The domain of the target.

class sumfem.symbolic.Assignment(name, target, expression)[source]#

Write an expression of fields into the target, point by point.

The expression is evaluated at the target’s interpolation points, so it may only use terminals on the target’s domain. It may use the target itself, e.g. the update g + lambda (T - T_m) of a liquid fraction g, clipped to [0, 1], between solves for the temperature T.

Parameters:

name (str) – See Transfer.
target (Transferable) – See Transfer.
expression (ufl.core.expr.Expr) – The expression, of the target’s shape. Values are in SI base units, as every symbolic expression’s.

Raises:

ComponentError – If the expression and the target differ in shape.

property source: Expr#: The expression written into the target.

class sumfem.symbolic.Interpolation(name, source, target, markers=())[source]#

Interpolate into the target, from a source on any mesh.

The meshes need not match or be nested. The discrete layer locates the target’s interpolation points in the source mesh once, and evaluates the source there at every transfer, see Interpolation.

Parameters:

name (str) – See Transfer.
source (Transferable) – See Transfer.
target (Transferable) – See Transfer.
markers (Sequence[str]) – Cell markers of the target’s domain to interpolate on. Empty for all its cells.

property markers: tuple[str, ...]#: Cell markers of the target’s domain to interpolate on.

class sumfem.symbolic.Domain(name, cell, gdim=None, degree=1, length_scale=1.0, lagrange_variant=None, coordinates=None)[source]#

A symbolic description of a type of mesh.

A domain has a name, a cell type and named markers, but does not store any data.

A domain is named independently of the models that use it, and is shared between models by identity. Several models written on the same Domain are discretized on the same mesh.

Subdomains are referred to by marker names, see measure(). The discretization later maps the names to the tags of an actual mesh.

Parameters:

name (str) – Name of the domain, used to pair it with a mesh.
cell (CellType) – Cell type, e.g. "triangle" or "hexahedron".
gdim (int | None) – Geometric dimension. Defaults to the topological dimension.
degree (int) – Degree of the coordinate element.
length_scale (float) – Typical length of the domain in metres, used as reference length in dimensional analysis. Mesh coordinates are in metres.
lagrange_variant (basix.LagrangeVariant | str | None) – Lagrange variant of the coordinate element. It must match the mesh, e.g. "equispaced" for meshes read from gmsh. from_mesh() takes it from a mesh.
coordinates (CoordinateSystem | None) – The coordinate system, which weights the integrals over the domain, e.g. Cylindrical for an axisymmetric body. Defaults to Cartesian.

Raises:

ComponentError – For a name that is not an identifier, a length scale that is not positive, or a coordinate system the geometric dimension does not fit.

classmethod from_mesh(name, mesh, length_scale=1.0, coordinates=None)[source]#

A domain whose coordinate element matches an existing mesh.

Parameters:

name (str) – Name of the domain.
mesh (Mesh) – The mesh, which provides cell type, geometric dimension, coordinate degree and Lagrange variant.
length_scale (float) – See Domain.
coordinates (CoordinateSystem | None) – See Domain.

Return type:

Domain

property name: str#: Name of the domain.

property cell_name: str#: Name of the cell type, e.g. "triangle".

property length_scale: float#: Typical length of the domain in metres, for dimensional analysis.

property coordinates: CoordinateSystem#: The coordinate system, which weights the integrals.

property tdim: int#: Topological dimension.

property gdim: int#: Geometric dimension.

entity_dimension(entity_type)[source]#

Topological dimension of an entity type on this domain.

Raises:: ComponentError – If the entity type does not exist on this domain, for instance a ridge on an interval mesh.
Parameters:: entity_type (Literal['cell', 'facet', 'ridge', 'peak', 'vertex'])
Return type:: int

marker_ids(markers, entity_type)[source]#

Return the integer ids of named markers, assigning ids to new names.

The ids are numbered separately per entity dimension, so entity types that coincide share their markers: on a triangle mesh a ridge is a vertex, and "A" has the same id for both.

Parameters:

markers (str | Sequence[str]) – Name of a subdomain, e.g. "Electrode", or a sequence of names.
entity_type (Literal['cell', 'facet', 'ridge', 'peak', 'vertex']) – The entities the markers refer to.

Returns:

One id per marker, in the given order.

Return type:

tuple[int, …]

property markers: dict[int, dict[str, int]]#: Copy of the marker name to integer id maps, per entity dimension.

measure(integral_type, *markers, degree=None, metadata=None)[source]#

Create a measure over the whole domain or over named subdomains.

Subdomains are referred to by marker name, e.g. domain.measure("dx", "Electrode"), or with the shortcuts domain.dx("Electrode") and domain.ds("Inlet"). Entity types that coincide share their markers, see marker_ids(): on a triangle mesh dr("A") and dP("A") refer to the same vertices.

Parameters:

integral_type (Literal['dx', 'ds', 'dS', 'dr', 'dP']) – One of "dx", "ds", "dS", "dr" or "dP".
markers (str) – Names of the subdomains to integrate over. With no names, the measure covers the whole domain.
degree (int | None) – Quadrature degree. Shorthand for metadata={"quadrature_degree": degree}.
metadata (dict | None) – Form compiler metadata for the measure.

Return type:

Measure

dx(*markers, degree=None, metadata=None)#

Cell measure over the domain, or over the named cell markers.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

ds(*markers, degree=None, metadata=None)#

Exterior facet measure, optionally restricted to named facet markers.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

dS(*markers, degree=None, metadata=None)#

Interior facet measure, optionally restricted to named facet markers.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

dr(*markers, degree=None, metadata=None)#

Ridge measure (codimension 2), optionally restricted to named ridges.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

dP(*markers, degree=None, metadata=None)#

Vertex measure, optionally restricted to named vertex markers.

Parameters:

markers (str)
degree (int | None)
metadata (dict | None)

Return type:

ufl.Measure

one_sided(marker)[source]#

The subdomain a facet marker sees its facets from, if one-sided.

Facet integrals of a subdomain of cells move to its parent as integrals over facets seen from the subdomain’s side, see Subdomain.on_parent(), under markers named "<subdomain>.<marker>". Inside the parent, such facets may have cells on several sides, e.g. at the joint of a branching surface.

Returns:: The subdomain and its own facet marker, None for its whole boundary, or None if the marker is not one-sided.
Parameters:: marker (str)
Return type:: tuple[Subdomain, str | None] | None

interface(first, second, marker, degree=None, metadata=None)[source]#

Interior facet measure of an interface between two subdomains.

Its "+" side is always in first and its "-" side in second, so the functions of each subdomain can be restricted to their own side, e.g. u_0("+") and u_1("-") for an exchange flux between two materials.

Parameters:

first (Subdomain) – The subdomain of cells on the "+" side.
second (Subdomain) – The subdomain of cells on the "-" side.
marker (str) – This domain’s facet marker of the interface.
degree (int | None) – Quadrature degree, see measure().
metadata (dict | None) – Form compiler metadata, see measure().

Raises:

ComponentError – If a subdomain is not a subdomain of cells of this domain.

Return type:

Measure

interface_sides(marker)[source]#

The subdomains and facet marker of an interface measure’s marker.

See interface(). None if the marker is not an interface.

Parameters:: marker (str)
Return type:: tuple[Subdomain, Subdomain, str] | None

class sumfem.symbolic.Field(*args, **kw)[source]#

An unknown of a variational problem.

A field is the coefficient holding the solution, so models use it directly in expressions, e.g. ufl.grad(c). It also owns:

v: the test function of its equation,
du: the trial function, for hand-written Jacobians,
dot: a placeholder for its time derivative,
previous(k): its value k time steps back.

Models that write terms with the same field.v contribute to the same equation when they are summed. A model that uses the field without writing an equation for it is coupled to whichever model does.

The name must be passed positionally: ufl.Coefficient.__new__ reads the first positional argument.

Parameters:

name (str) – Unique name of the field.
domain (Domain) – The domain the field lives on.
element (ElementSpec) – The finite element, see ElementSpec.
unit (Unit | None) – Unit of the field, a sympy unit expression. None means dimensionless. Values given for the field, such as initial and boundary values, are in this unit. The discretization stores them in SI base units.
scale (float) – Typical magnitude of the field in its unit, used as reference value in dimensional analysis.
spatial (bool) – Whether the components of a vector-valued field are those of a vector in physical space, which change with the basis of the coordinates, e.g. a velocity or a displacement, or independent scalars, which do not, e.g. the concentrations of several species. The two differ only on domains whose coordinate system has a curved basis, e.g. Cylindrical: there a vector in space gains components in the angular directions, and its gradient the terms of the Christoffel symbols. In Cartesian coordinates, and with a weight alone, they are the same. Elements mapped by a Piola transformation, e.g. Raviart-Thomas or Nedelec, are always vectors in space.

Raises:

ComponentError – For spatial=False with an element mapped by a Piola transformation.

property name: str#: Name of the field.

property spatial: bool#: Whether the components are those of a vector in physical space.

property domain: Domain#: Domain the field lives on.

property element: AbstractFiniteElement#: Finite element of the field.

property unit: Expr#: Unit of the field, 1 if dimensionless.

property scale: float#: Typical magnitude of the field in its unit, for dimensional analysis.

property v: Argument[source]#

The test function of this field’s equation.

Its part number is the field’s count, unique among all fields.

property du: Argument[source]#: The trial function, for hand-written Jacobians.

property dot: FieldTerminal[source]#

Placeholder for the time derivative of the field.

A time-stepping scheme replaces it, for instance by a backward difference of the field and its previous() values.

previous(k=1)[source]#

The value of the field k time steps back.

Parameters:: k (int)
Return type:: FieldTerminal

class sumfem.symbolic.KnownCoefficient(*args, **kw)[source]#

Named, spatially varying data that a model reads but does not solve for.

It is also a ufl_units quantity, so dimensional analysis finds it in forms, with its typical magnitude as scale.

Parameters:

name (str) – Unique name. Known coefficients created through known_coefficient() are named "<component>.<name>".
domain (Domain) – The domain the data is defined on.
element (ElementSpec) – The finite element, see ElementSpec.
unit (Unit | None) – Unit of the data, a sympy unit expression. None means dimensionless. The discretization stores values in SI base units.
scale (float) – Typical magnitude of the data in its unit, for dimensional analysis.
doc (str) – Description of the data.

The name must be passed positionally: ufl.Coefficient.__new__ reads the first positional argument.

property name: str#: Unique name.

property element: AbstractFiniteElement#: Finite element of the data.

property doc: str#: Description of the data.

class sumfem.symbolic.Residual(name)[source]#

Base class for a reusable contribution to a residual.

Subclasses override residual() and write it with the fields they act on and their test functions field.v. Everything the residual uses is discovered from the form, so there is no registration of unknowns, parameters or known coefficients. A subclass that overrides one term only changes that term.

The residual is built on first use and then cached, so a component should be fully configured before it is summed or turned into a system.

Parameters:: name (str) – Unique name of the component. It prefixes the names of parameters and known coefficients the component creates, as "<name>.<parameter>".

residual()[source]#

Return this component’s contribution to the residual.

The form must be linear in test functions field.v and contain no trial functions. A contribution that is zero, e.g. for some inputs, is ufl.ZeroBaseForm((field.v,)): UFL drops the test function from a form whose terms vanish.

Return type:: Form | ZeroBaseForm

jacobian()[source]#

Return a hand-written linearization of residual(), or None.

A Jacobian, or a form with field.du as trial functions. Blocks for unknowns it does not mention are zero. With None, the Jacobian is derived automatically.

Return type:: Jacobian | Form | None

preconditioner()[source]#

Return this component’s contribution to a preconditioner, or None.

If any summand provides a preconditioner, summands without one contribute their Jacobian.

Return type:: Jacobian | Form | None

property components: tuple[Residual, ...]#: The summands. A single component is its own summand.

property form: Form | ZeroBaseForm[source]#: The residual, built on first access and then cached.

property jacobian_term: Jacobian | None[source]#: The hand-written Jacobian, if any, built on first access.

property preconditioner_term: Jacobian | None[source]#: The preconditioner contribution, if any, built on first access.

property inventory: Inventory[source]#: The named objects used by this component, validated.

property fields: tuple[Field, ...]#: Every field the component uses, in creation order.

property equations: tuple[Field, ...]#: Fields the component writes equation terms for, in creation order.

property parameters: Mapping[str, Parameter]#: Parameters used by the component, by name.

property known_coefficients: Mapping[str, KnownCoefficient]#: Known coefficients used by the component, by name.

property domains: tuple[Domain, ...]#: Domains used by the component.

to_jacobian(fields=None, name=None)[source]#

The derived linearization sum_f derivative(F, f, f.du).

It is derived after the substitutions of the system it is used in, e.g. of a time-stepping scheme, so the terms of time derivatives are linearized too.

Parameters:

fields (Sequence[Field] | None) – The fields to linearize in. Defaults to the unknowns of the system it is used in, or, on its own, the fields this residual writes equations for.
name (str | None) – Name of the Jacobian. Defaults to the residual’s name.

Return type:

Jacobian

with_jacobian(jacobian)[source]#

This residual, linearized by the given Jacobian.

Parameters:: jacobian (Jacobian)
Return type:: Residual

with_preconditioner(preconditioner)[source]#

This residual, with the given preconditioner.

It replaces the whole preconditioner of this residual. Its Jacobian is unchanged, including hand-written Jacobians of its summands.

Parameters:: preconditioner (Jacobian)
Return type:: Residual

copy(suffix)[source]#

Return a copy of this component that solves for new fields.

Each field the component writes an equation for is replaced by a new field named "<name><suffix>", with the same domain, element, unit and scale. Its test and trial functions, time derivative and previous values are replaced too. The new fields have their own UFL counts, so the copy can be summed with the original.

Each summand is copied as "<name><suffix>", together with the parameters and known coefficients it created. They keep their defaults, units and docs, under the new name. Everything else the component uses is shared with the original, such as fields it has no equation for and inputs given to its constructor.

The copy has the original’s forms with these replacements. It is not an instance of the original’s class, so attributes set by the original’s constructor are not available on it.

Parameters:: suffix (str) – Appended to the names of the copied fields, components, parameters and known coefficients, e.g. "_pos".
Raises:: ComponentError – If the suffix is empty.
Return type:: Residual

Example

The positive half-cell of a flow battery has the same equations as the negative one:

positive = negative.copy("_pos")
cell = negative + positive

class sumfem.symbolic.Parameter(*args, **kw)[source]#

A named value that is constant in space, with a unit.

The value lives on a Real function space, a single global degree of freedom per component. This means a residual can be differentiated with respect to a parameter, as needed for sensitivities and optimal control. The discretization owns the degree of freedom on one process and keeps all processes in sync.

A parameter is also a ufl_units quantity, so dimensional analysis finds it in forms. Its scale for that analysis is the largest magnitude of its default, or 1 without a non-zero default.

Parameters:

name (str) – Unique name. Parameters created through parameter() are named "<component>.<name>".
domain (Domain) – The domain the parameter is defined on.
shape (tuple[int, ...]) – Value shape, () for a scalar.
default (npt.ArrayLike | None) – Default value in unit, used when no other value is set. A scalar sets every component of a vector or tensor parameter.
unit (Unit | None) – Unit of the value, a sympy unit expression. None means dimensionless. The discretization stores values in SI base units.
doc (str) – Description, shown in documentation and archived metadata.

The name must be passed positionally: ufl.Coefficient.__new__ reads the first positional argument.

property name: str#: Unique name.

property shape: tuple[int, ...]#: Value shape.

property default: ndarray[tuple[Any, ...], dtype[Any]] | None#: Default value, or None if the value must be set explicitly.

property doc: str#: Description of the parameter.

class sumfem.symbolic.PiGroup(expression, value)[source]#

A dimensionless group of the quantities of a system.

Parameters:

expression (Expr)
value (float)

expression#

The group as a product of powers of quantity symbols.

Type:: sympy.core.expr.Expr

value#

Its value for the scales of the quantities.

Type:: float

class sumfem.symbolic.System(component, field_order=None, substitutions=mappingproxy({}), equation_scales=mappingproxy({}))[source]#

The equations of a residual, in a fixed block order.

A system is given by its component, the block order, substitutions and equation scales. Everything else is derived from these on first use. The unknowns are the fields the component writes equations for, its data are the other terminals of its forms, and its forms are blocked over a ufl.MixedFunctionSpace of the unknowns. build_system() creates systems, and also checks units and computes equation scales.

With one unknown the residual and Jacobian are single forms. With several unknowns the residual is a tuple with one form per unknown, and the Jacobian a tuple of rows with None for zero blocks. Each integral is weighted by the coordinate system of its domain, see coordinates.

component#

The residual.

Type:: sumfem.symbolic.components.Residual

field_order#

Block order, a permutation of the unknowns given as fields or names. None orders them by creation.

Type:: Sequence[sumfem.symbolic.terminals.Field | str] | None

substitutions#

Replacements in every form before the Jacobian is derived, such as those of a time-stepping scheme. Objects they introduce become part of the system.

Type:: Mapping[ufl.coefficient.Coefficient, ufl.core.expr.Expr]

equation_scales#

The size each equation is divided by, in SI base units, by the name of its unknown. Empty for dimensional equations.

Type:: Mapping[str, float]

Raises:

SystemBuildError – If the component has no equations, if field_order is not a permutation of the unknowns, if a hand-written Jacobian or preconditioner uses the trial function of a field that is not an unknown, if equation_scales does not give one positive scale per unknown, if a vector in space is an unknown on a domain whose coordinate system supports scalar unknowns only, or if a tensor-valued unknown lies on a domain whose coordinate system is not Cartesian.

Parameters:

component (Residual)
field_order (Sequence[Field | str] | None)
substitutions (Mapping[Coefficient, Expr])
equation_scales (Mapping[str, float])

property inventory: Inventory[source]#: The named objects of the component’s forms, after substitution.

property unknowns: tuple[Field, ...][source]#: The fields the component writes equations for, in block order.

property parameters: Mapping[str, Parameter]#: Parameters used by the system, by name.

property known_coefficients: Mapping[str, KnownCoefficient | Field][source]#

The data of the system, by name.

These are its known coefficients, and the fields it uses but does not solve for, such as fields solved for in another stage.

property domains: tuple[Domain, ...]#: Domains used by the system.

property test_functions: tuple[Argument, ...][source]#: The test function of each block, as used in the residual.

property trial_functions: tuple[Argument, ...][source]#: The trial function of each block, as used in the Jacobian.

property residual: Form | tuple[Form, ...]#: Residual form, or one form per unknown.

property jacobian: Form | tuple[tuple[Form | None, ...], ...]#: Jacobian form, or blocks J[i][j] = dF_i/du_j.

property preconditioner: Form | tuple[tuple[Form | None, ...], ...] | None#: Preconditioner form or blocks, or None.

property blocked: bool#: Whether the system has more than one unknown.

index(field)[source]#

Block index of an unknown, given as a field or by name.

Parameters:: field (Field | str)
Return type:: int

derivative(wrt)[source]#

Derivative of the residual with respect to a parameter or data.

The result is linear in a trial function on the space of wrt, giving the sensitivity forms dF_i/dm used by adjoint gradients.

Parameters:: wrt (Parameter | KnownCoefficient | Field | str) – A parameter or a known coefficient of this system, given as the object or by name.
Returns:: One form, or one form per unknown with None where the residual block does not depend on wrt.
Return type:: Form | tuple[Form | None, …] | None

exception sumfem.symbolic.UnitError[source]#: Raised when an equation is dimensionally inconsistent.

sumfem.symbolic.as_expr(value)[source]#

Convert a component input to a UFL expression.

UFL expressions, including a Field, Parameter and KnownCoefficient, are returned unchanged. Numbers and nested sequences of numbers become UFL literals.

Parameters:: value (Expr | Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
Return type:: Expr

sumfem.symbolic.build_system(component, field_order=None, substitutions=None, check_units=False, factorize_equations=False, time_scale=1.0)[source]#

Build the system of equations for a residual.

The component is not modified.

Parameters:

component (Residual) – A residual, or a sum of residuals. The unknowns are the fields it writes an equation for, with a term in field.v. Other fields used by the component are known coefficients of the system, as in a stage of a sequence.
field_order (Sequence[Field | str] | None) – Block order, a permutation of the unknowns, as objects or names. Defaults to field creation order.
substitutions (Mapping[Coefficient, Expr] | None) – Coefficients to replace in every form before the Jacobian is derived, such as {c.dot: (c - c.previous()) / dt} from a time-stepping scheme. Objects introduced by a substitution become part of the system.
check_units (bool) – Check that every equation is dimensionally consistent, see check_units().
factorize_equations (bool) – Make the equations dimensionless. Each equation is divided by its typical size, see equation_scales(), by dividing the test function field.v of its unknown. The unknowns keep their SI values. System.equation_scales records the sizes. This checks units too, since the sizes are only meaningful for consistent equations.
time_scale (float) – Reference time in seconds for time derivatives, used by the unit check and the equation sizes.

Raises:

SystemBuildError – If the component has no equations, if a hand-written Jacobian uses the trial function of a field that is not an unknown, if field_order is not a permutation of the unknowns, if a vector in space is an unknown on a domain whose coordinate system supports scalar unknowns only, or if a tensor-valued unknown lies on a domain whose coordinate system is not Cartesian.
UnitError – If units are checked, or equations factorized, and an equation is dimensionally inconsistent.

Return type:

System

sumfem.symbolic.check_units(system, time_scale=1.0)[source]#

Check that every term of every equation has the same dimension.

Parameters:

system (System) – The system.
time_scale (float) – Reference time in seconds, for time derivatives.

Raises:

UnitError – Naming the first dimensionally inconsistent equation.

Return type:

None

sumfem.symbolic.equation_dimensions(system, time_scale=1.0)[source]#

The dimension of each equation of a system, by field name.

Parameters:

system (System) – The system.
time_scale (float) – Reference time in seconds, for time derivatives.

Raises:

UnitError – If an equation is dimensionally inconsistent.

Return type:

dict[str, Expr]

sumfem.symbolic.equation_scales(system, time_scale=1.0)[source]#

The typical size of each equation of a system, by field name.

Each term of an equation gets a typical size from the scales of the quantities in it. The size of the equation is that of its largest term, in SI base units. Dividing the equation by it makes its dominant terms of order one.

Parameters:

system (System) – The system.
time_scale (float) – Reference time in seconds, for time derivatives.

Raises:

UnitError – If an equation is dimensionally inconsistent.

Return type:

dict[str, float]

sumfem.symbolic.pi_groups(system, time_scale=1.0)[source]#

The dimensionless groups of a system, by Buckingham Pi analysis.

Parameters:

system (System) – The system.
time_scale (float) – Reference time in seconds, for time derivatives.

Return type:

list[PiGroup]

The discrete layer#

Binding symbolic components to meshes, as their discrete counterparts.

MeshData: one mesh with its mesh tags, and names for its markers. It knows nothing of symbolic domains.
DomainMaps: pairs each symbolic domain with a MeshData, checks that they match, and creates the submeshes of subdomains.
discretize(): the Discretization of symbolic components on domain maps, holding one DOLFINx function per terminal, shared by all of them. disc[component] gives a component’s discrete counterpart.
System: the residual, Jacobian and preconditioner of a symbolic system, compiled and bound to the meshes and functions.
DirichletBC: a symbolic boundary condition on the meshes.
counterpart(): the discrete counterpart of a symbolic object, with options.

Solving the counterparts is left to sumfem.solvers.

A discrete counterpart has the name of its symbolic one, and is made by indexing a discretization, by counterpart(), or by its own constructor from the discretization and the symbolic object. sumfem exports the symbolic ones.

class sumfem.discrete.System(discretization, system)[source]#

The residual, Jacobian and preconditioner of a symbolic system, compiled and bound to the meshes and functions of a discretization.

The forms are evaluated at the current values of the functions, so a Solver assembles them again at every Newton iteration. A linear system is the case where one iteration suffices. The forms are compiled with the discretization’s compiler options, and cached per symbolic system and options.

Parameters:

discretization (Discretization) – The discretization holding the system’s terminals.
system (symbolic.System) – The symbolic system.

property symbolic: System#: The symbolic system.

property discretization: Discretization#: The discretization owning the meshes and functions.

property residual: Form | list[Form]#: Residual form, or one form per unknown.

property jacobian: Form | list[list[Form | None]]#: Jacobian form, or blocks with None for zero blocks.

property preconditioner: Form | list[list[Form | None]] | None#: Preconditioner forms, or None.

property solution: tuple[Function, ...]#: The function of each unknown, in block order.

class sumfem.discrete.Discretization(maps, max_quadrature_degree=None, form_compiler_options=None, jit_options=None)[source]#

The discrete state of symbolic components on domain maps.

Made by discretize(). It holds what its components need: the submesh of each subdomain they use, in the maps, and a function for each field, parameter, known coefficient and previous value in them. These are created when the components are declared, by discretize() or add(). A terminal not held raises, so nothing is created behind a lookup.

Indexing gives discrete counterparts, compiled on first access and kept: see __getitem__(). Every system of a discretization shares its functions and compiler options.

Declaring is collective: every process must declare the same components in the same order, since creating a submesh is.

Parameters:

maps (DomainMaps) – The meshes of the domains.
max_quadrature_degree (int | None) – Upper bound on the quadrature degree UFL estimates for each integral without an explicit degree, through FFCx’s max_quadrature_degree metadata. Use it to keep nonlinear integrands cheap. None for no bound.
form_compiler_options (dict | None) – Options for FFCx.
jit_options (dict | None) – Options for the JIT compilation of the generated code.

Raises:

ValueError – If the bound is negative.

add(*components)[source]#

Declare components, creating the submeshes and functions they need.

A component is a residual, functional or Jacobian, a symbolic system, a Dirichlet condition, a transfer, a sequence, a clock, a form, an expression, a single terminal, or a domain, e.g. a subdomain whose submesh is needed. A system holds the terminals a time-stepping scheme introduced, such as previous values and weights, which its residual does not. A parameter starts at its default value. Terminals and submeshes held already are kept.

Raises:

KeyError – If there is no mesh for a domain they use.
TypeError – For something that is not a component.

Parameters:

Return type:

None

property maps: DomainMaps#: The meshes of the domains.

property max_quadrature_degree: int | None#: Upper bound on estimated quadrature degrees, or None.

property form_compiler_options: dict#: Options for FFCx.

property jit_options: dict#: Options for the JIT compilation of the generated code.

property comm: Intracomm#: The MPI communicator of the meshes.

function_space(field)[source]#

The function space of a field, distinct from that of any other field.

Parameters:: field (Field)
Return type:: FunctionSpace

function(terminal)[source]#

The function of a field, parameter, known coefficient or history value.

Raises:

SystemBuildError – For the time-derivative placeholder field.dot, which a time-stepping scheme must replace.
KeyError – If the discretization does not hold the terminal.

Parameters:

terminal (Coefficient)

Return type:

Function

property fields: Mapping[str, Function]#: The function of each field held, by field name.

property parameters: Mapping[str, Function]#: The function of each parameter held, by parameter name.

property known_coefficients: Mapping[str, Function]#: The function of each known coefficient held, by name.

terminal(key)[source]#

The symbolic terminal with a name, among those held.

A terminal given as an object is returned as it is.

Raises:: KeyError – If no terminal held has the name.
Parameters:: key (str | Coefficient)
Return type:: Coefficient

set(values)[source]#

Set the values of fields, parameters or known coefficients.

Values are given in the unit of the terminal and stored in SI base units. A sympy quantity with its own units, e.g. 5 * sympy.physics.units.centimeter, is converted from those.

A symbolic UFL expression, of fields, parameters, known coefficients and spatial coordinates, sets a field or known coefficient by interpolation, e.g. an initial condition. It evaluates in SI base units, as every symbolic expression does.

Parameters:

Raises:

ValueError – If a sympy quantity has the wrong dimension, or an expression the wrong shape.
TypeError – For an expression given to a parameter.

Return type:

None

compile_expression(expr, space)[source]#

Compile a symbolic expression for interpolation into a space.

The expression is discretized with discretize() and compiled at the interpolation points of the space, with this discretization’s compiler options and entity maps. On a domain whose coordinate system writes vectors in 3D, see coordinates, the expression is evaluated in 3D, and a space of the mesh’s dimension receives its components in the mesh’s coordinates, e.g. those of a stress but its hoop component.

Raises:

ValueError – If the expression’s shape is not the space’s.

Parameters:

expr (Expr)
space (FunctionSpace)

Return type:

Expression

value(parameter)[source]#

The current value of a parameter, in the parameter’s unit.

Parameters:: parameter (Parameter | str)
Return type:: ndarray[tuple[Any, …], dtype[Any]]

has_value(parameter)[source]#

Whether a parameter has a value, from its default or set explicitly.

Parameters:: parameter (Parameter | str)
Return type:: bool

check_values(system)[source]#

Raise if a parameter of a system has no value.

Raises:: ValueError – Naming the parameters without a value.
Parameters:: system (System)
Return type:: None

history(field)[source]#

The functions of the field’s previous values held, by k.

Parameters:: field (Field)
Return type:: dict[int, Function]

initialize_history()[source]#

Set every previous value of every field to the field’s current value.

Return type:: None

advance_history()[source]#

Shift the previous values of every field one step back in time.

Afterwards previous(1) holds the current value of the field.

Return type:: None

assemble_scalar(form)[source]#

Assemble a symbolic functional, summed over all processes.

The functional is compiled on first use and kept, so assembling it again, e.g. at every time step, uses the current values of its terminals without compiling it again.

Parameters:: form (Form | Functional) – A form without test or trial functions over one domain, e.g. c * domain.ds("Membrane"), or a functional.
Return type:: float

class sumfem.discrete.DomainMaps(meshes, entity_maps=())[source]#

Pairs each symbolic domain with a MeshData, and checks that they match.

The maps also hold the submesh of each subdomain, created with its markers and its entity map from the parent’s mesh when a discretization first needs it, unless a mesh is given for the subdomain, and shared by every discretization on the maps. Integration entities of the markers are computed on first use.

Parameters:

meshes (Mapping[Domain | str, MeshLike]) – The mesh of each domain, by the domain or its name. A mesh given for a domain is checked against it here, one given by name when the domain is first used.
entity_maps (Iterable[dolfinx.mesh.EntityMap]) – Entity maps between meshes, for forms that couple fields on different meshes.

Raises:

ValueError – If there are no meshes, a domain is given twice, or a mesh does not match its domain.

property entity_maps: list[EntityMap]#: Entity maps between meshes, passed to every form.

property comm: Intracomm#: The MPI communicator of the meshes.

ensure_mesh(domain)[source]#

Make sure a domain has a mesh, creating the submesh of a subdomain if needed.

Collective: a submesh is created on every process together. Use mesh_data() or mesh() to get the mesh, which never create one, and so can be called on any process.

Raises:

KeyError – If there is no mesh for a domain that is no subdomain.
ValueError – If the mesh does not match the domain.

Parameters:

domain (Domain)

Return type:

None

mesh_data(domain)[source]#

The mesh of a domain, with its markers.

The submesh of a subdomain must have been created by ensure_mesh(), which a discretization does for the domains of what it holds.

Raises:

KeyError – If there is no mesh for the domain.
ValueError – If the mesh does not match the domain’s cell type, coordinate degree or geometric dimension.

Parameters:

domain (Domain)

Return type:

MeshData

mesh(domain)[source]#

The mesh of a domain.

Parameters:: domain (Domain)
Return type:: Mesh

integration_entities(domain, integral_type, marker)[source]#

Integration entities of a named marker, for a UFL integral type.

The entities are in the flattened layout DOLFINx assembles over, e.g. cells for a cell integral, (cell, local facet) pairs for an exterior facet integral, and (cell, local facet, cell, local facet) for an interior facet integral. The marker is one of:

a marker of the domain’s mesh, see MeshData.markers;
a one-sided marker "<subdomain>.<marker>" of an exterior facet integral moved from a subdomain to its parent, see one_sided();
an interface marker of an interior facet integral between two subdomains, see interface().

They are computed on first use and cached.

Raises:

KeyError – If there is no mesh for the domain.
ComponentError – If the facets of the marker cannot be integrated over with the integral type, e.g. an exterior facet marker with facets inside the mesh.

Parameters:

domain (Domain)
integral_type (str)
marker (str)

Return type:

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

check_interior_facets(domain, marker)[source]#

Raise unless the interior facets of a domain can be integrated over.

Each needs exactly two cells, and both on the process. A submesh of facets can have more, where three or more of its cells meet, which an interior facet integral does not cover. In parallel, the cells of a submesh of codimension one or more that meet at a facet on a process boundary share only a ridge of the parent, so the parent must be ghosted by vertex, see MeshData.ghost_by_vertex().

Parameters:

domain (Domain) – The domain.
marker (str | None) – A facet marker of the domain, or None for all facets.

Raises:

ComponentError – On every process, if a process has such a facet.

Return type:

None

class sumfem.discrete.MeshData(mesh, tags=mappingproxy({}), markers=mappingproxy({}))[source]#

One mesh with its mesh tags, and names for its markers.

It knows nothing of symbolic domains: DomainMaps pairs domains with meshes. For each entity dimension, a marker name maps to tag values. Create it from a locator function for each dimension and name with from_locators(), from existing mesh tags for each dimension and the integer tag values of each name directly, e.g. MeshData(mesh, {1: facet_tags}, {1: {"Inlet": (1,), "Walls": (2, 3)}}), or from the physical groups of a mesh read with dolfinx.io.gmsh with from_gmsh().

mesh#

The mesh.

Type:: dolfinx.mesh.Mesh

tags#

Mesh tags, by entity dimension.

Type:: collections.abc.Mapping[int, dolfinx.mesh.MeshTags]

markers#

The tag values of each named marker, by entity dimension.

Type:: collections.abc.Mapping[int, collections.abc.Mapping[str, tuple[int, …]]]

Raises:

ValueError – If mesh tags are given for another dimension than their own, or markers for a dimension without mesh tags.

Parameters:

mesh (Mesh)
tags (Mapping[int, MeshTags])
markers (Mapping[int, Mapping[str, tuple[int, ...]]])

classmethod from_gmsh(data)[source]#

Use the physical groups of a mesh read with dolfinx.io.gmsh.

Each physical group becomes a marker with the group’s name.

Parameters:: data (MeshData) – The mesh, its mesh tags and its physical groups, e.g. from dolfinx.io.gmsh.read_from_msh().
Return type:: MeshData

classmethod from_locators(mesh, locators, boundary=False)[source]#

Mark entities geometrically, as dolfinx.mesh.locate_entities().

Parameters:

mesh (Mesh) – The mesh.
locators (Mapping[int, Mapping[str, Callable[[ndarray[tuple[Any, ...], dtype[floating]]], ndarray[tuple[Any, ...], dtype[bool]]]]]) – For each entity dimension, a locator per marker name. Where locators overlap, the later marker wins.
boundary (bool) – Mark only entities on the boundary of the mesh, for dimensions below the cell dimension. Use this for boundary markers, so that interior facets touching the boundary are not marked.

Return type:

MeshData

entities(dim, marker)[source]#

Local indices of the entities of dimension dim with a marker.

Raises:

KeyError – If the mesh has no such marker.

Parameters:

dim (int)
marker (str)

Return type:

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

unmarked_boundary_facets(markers)[source]#

Local indices of the boundary facets that none of the markers covers.

A hyperbolic system needs a boundary flux on every boundary facet, so the markers of its boundary models should leave none.

Parameters:: markers (Iterable[str])
Return type:: ndarray[tuple[Any, …], dtype[int32]]

submesh(dim, marker)[source]#

The submesh of the entities of a dimension with a marker.

The submesh keeps the markers of its entities: the parent’s markers of dimension dim and lower become those of the submesh, e.g. a surface of a 3D mesh has the parent’s facet markers as cell markers and its ridge markers as facet markers.

Returns:

The submesh, and the entity map from its cells to the parent’s entities.

Parameters:

dim (int)
marker (str)

Return type:

tuple[MeshData, EntityMap]

ghost_by_vertex()[source]#

The mesh with every cell ghosted on each process sharing its vertices.

See ghost_by_vertex(). The markers carry over, onto the new ghost entities too. On one process the mesh is returned.

Return type:: MeshData

sumfem.discrete.ghost_by_vertex(mesh)[source]#

The mesh, with every cell ghosted on each process sharing its vertices.

Interior facet integrals on a submesh of facets need both cells of each of its facets on the process, but cells that share only a ridge are not ghosted by default. With this ghost layer, all cells around each owned vertex are on the process. The cells keep their owners. The new mesh has the global node numbers of mesh as its input node indices, see mesh.geometry.input_global_indices.

After J. S. Dokken, https://gist.github.com/jorgensd/1e3d89b8c46a05fce36eb2fe88cfe9a6.

Parameters:: mesh (Mesh)
Return type:: Mesh

class sumfem.discrete.DirichletBC(discretization, condition)[source]#

A DirichletBC discretized on its field’s mesh.

A UFL value is compiled into an expression, evaluated by pack(). Any other value is assigned once.

Parameters:

discretization (Discretization) – The discretization of the field.
condition (symbolic.DirichletBC) – The symbolic condition.

property symbolic: DirichletBC#: The symbolic condition.

property function: Function#: The function holding the prescribed values.

property dolfinx_bc: DirichletBC#: The DOLFINx boundary condition.

property expression: Expression | None#: The compiled expression of a symbolic value, or None.

pack()[source]#

Evaluate a symbolic value with the current values of its inputs.

Return type:: None

sumfem.discrete.counterpart(symbolic, discretization, **options)[source]#

sumfem.discrete.counterpart(condition, discretization, **options)

sumfem.discrete.counterpart(expr, discretization, **options)

sumfem.discrete.counterpart(system, discretization, **options)

sumfem.discrete.counterpart(residual, discretization, **options)

sumfem.discrete.counterpart(form, discretization)

sumfem.discrete.counterpart(functional, discretization, **options)

The discrete counterpart of a symbolic object on a discretization.

Dispatches on the type of the symbolic object. A discrete counterpart has the name of its symbolic one in sumfem.discrete, e.g. a symbolic System has a System, and a symbolic Interpolation a Interpolation. Expressions give expressions in the discrete functions, and functionals dolfinx.fem.Form.

Indexing a discretization, disc[symbolic], gives the counterpart with default options and keeps it. Call this function for other options, e.g. counterpart(interpolation, disc, padding=1e-10). Register the discrete counterpart of a new symbolic type, e.g. a transfer by projection, with @counterpart.register:

@counterpart.register
def _(projection: Projection, discretization, **options):
    return DiscreteProjection(discretization, projection, **options)

Parameters:

symbolic (Any) – The symbolic object.
discretization (Discretization) – The discretization holding its terminals.
options (Any) – Discrete options of the counterpart, e.g. the point search of an interpolation.

Raises:

NotImplementedError – If no discrete counterpart is registered for the type.

Return type:

Any

sumfem.discrete.discretize(components, maps, max_quadrature_degree=None, form_compiler_options=None, jit_options=None)[source]#

The discretization of components on domain maps.

Components that are discretized together share their functions, so declare every component that couples to another here, e.g. a flow and the transport it carries. Discretization.add() declares later additions.

Parameters:

components (Domain | Residual | Functional | Jacobian | System | DirichletBC | Transfer | Sequential | Clock | Form | Expr | Sequence[Domain | Residual | Functional | Jacobian | System | DirichletBC | Transfer | Sequential | Clock | Form | Expr]) – A component, or a list or tuple of them, see Discretization.add().
maps (DomainMaps) – The meshes of the domains.
max_quadrature_degree (int | None) – Upper bound on estimated quadrature degrees, see Discretization.
form_compiler_options (dict | None) – Options for FFCx.
jit_options (dict | None) – Options for the JIT compilation of the generated code.

Return type:

Discretization

Example

A flow, and the half-cell it feeds:

disc = discretize([flow, cell.model, *inlet], maps)
Solver(disc[flow], [outlet]).solve()
Simulation(disc, cell.model, inlet, clock=clock, scheme=BDF(1)).run()

The solvers layer#

Solving the discrete counterparts of sumfem.discrete.

Solver: Newton’s method for one discrete System, with its Dirichlet conditions.

class sumfem.solvers.Solver(system, bcs=(), petsc_options=None, petsc_options_prefix=None, kind=None)[source]#

Solve a System.

Parameters:

system (discrete_system.System) – The discrete system.
bcs (Iterable[DirichletBC | dolfinx.fem.DirichletBC]) – Dirichlet conditions, symbolic or as DOLFINx objects. Symbolic conditions must be on unknowns of the system.
petsc_options (Mapping[str, Any] | None) – PETSc options for the SNES solver, on top of DEFAULT_PETSC_OPTIONS, and of FACTORIZED_PETSC_OPTIONS for factorized equations. An LU solve without pc_factor_mat_solver_type gets the package of direct_solver_options().
petsc_options_prefix (str | None) – Options prefix of the PETSc objects. Defaults to a unique prefix.
kind (str | None) – PETSc matrix and vector kind, e.g. "nest" for field-split preconditioning.

Symbolic boundary values are evaluated before every solve and every Newton iteration, so they follow parameters, time and known coefficients, as in Irksome’s DOLFINx backend.

property system: System#: The discrete system being solved.

property nonlinear_problem: NonlinearProblem#: The DOLFINx problem, giving access to the PETSc objects.

property boundary_values: tuple[Function | Constant, ...]#: The functions or constants holding the prescribed boundary values.

property iterations: int#: Newton iterations of the last solve.

property converged: bool#: Whether the last solve converged.

property converged_to_tolerance: bool#

Whether the last solve stopped at a residual or step tolerance.

Solving again from its result then changes nothing, unlike after a solve stopped at its iteration count, e.g. a smoother.

pack_bcs()[source]#

Evaluate symbolic boundary values with the current values of their inputs.

Return type:: None

solve()[source]#

Solve the system, starting from the current values of the unknowns.

Returns:: The function of each unknown, in block order, holding the solution in SI units.
Raises:: ValueError – If a parameter of the system has no value.
Return type:: tuple[Function, …]

Post-processing#

Post-processing of DOLFINx functions: integrals, norms, lines and ranges.

sumfem.post_processing.create_line(p0, p1, num_points, comm=<mpi4py.MPI.Intracomm object>)[source]#

A mesh of interval cells on the line segment from p0 to p1.

The points are created on process 0, and DOLFINx distributes the cells.

Parameters:

p0 (ndarray[tuple[Any, ...], dtype[floating]]) – The first point, with as many coordinates as the line’s geometry.
p1 (ndarray[tuple[Any, ...], dtype[floating]]) – The last point.
num_points (int) – Number of points, including both ends.
comm (Intracomm) – The communicator of the line mesh.

Return type:

Mesh

sumfem.post_processing.evaluate_at_points(function, points)[source]#

Evaluate a function at points, each on the process that owns it.

dolfinx.geometry.determine_point_ownership() gives each point in the mesh one owning process, also on the boundary between processes. Points outside the mesh are evaluated nowhere. Collective.

Parameters:

function (Function) – The function.
points (Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – The points, of shape (n, gdim), the same on every process.

Returns:

The indices of the points evaluated on this process, increasing, and the values there, of shape (len(indices), value size).

Return type:

tuple[ndarray[tuple[Any, …], dtype[int32]], ndarray[tuple[Any, …], dtype[float64]]]

sumfem.post_processing.evaluate_over_line(function, start, end, num_points)[source]#

Evaluate a function at equidistant points on a line segment.

Each point is evaluated on the process that owns it, see evaluate_at_points(). Collective.

Parameters:

function (Function) – The function.
start (Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – The first point, with as many coordinates as the mesh.
end (Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – The last point.
num_points (int) – Number of points, including both ends.

Returns:

The points on this process, of shape (n, 3), and the values there, of shape (n, value size).

Return type:

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

sumfem.post_processing.function_range(function, entities=None, dim=None)[source]#

The global minimum and maximum of the values of a scalar function.

These are the extreme values of the degrees of freedom, which bound the function for Lagrange elements of degree one or zero.

Parameters:

function (Function) – The function.
entities (ndarray[tuple[Any, ...], dtype[int32]] | None) – Local indices of mesh entities whose degrees of freedom bound the range, e.g. the cells of a subdomain or the facets of a surface. None for the whole mesh.
dim (int | None) – Topological dimension of the entities. Defaults to the cells.

Raises:

ValueError – For a function with several components.

Return type:

tuple[float, float]

class sumfem.post_processing.NormType(*values)[source]#

The kind of a norm(). Plain strings with these values work too.

l2 = 'l2'#: The Euclidean norm of the degrees of freedom of a function.

L2 = 'L2'#: sqrt((f, f)).

H1 = 'H1'#: sqrt((f, f) + (grad f, grad f)).

H10 = 'H10'#: The seminorm sqrt((grad f, grad f)).

Hdiv = 'Hdiv'#: sqrt((f, f) + (div f, div f)), for vectors.

Hcurl = 'Hcurl'#: sqrt((f, f) + (curl f, curl f)), for vectors.

sumfem.post_processing.Squared#

the square of the integrand of a quantity f.

Type:: A custom norm

alias of Callable[[Expr], Expr]

sumfem.post_processing.error(function, exact, kind=NormType.L2, relative=False, coordinates=None)[source]#

The error of a function against an exact solution, in a norm().

Parameters:

function (Function) – The approximation.
exact (Expr | Callable[[SpatialCoordinate], Expr]) – The exact solution, as an expression on the function’s mesh or as a callable of the mesh’s spatial coordinate.
kind (NormType | str | Callable[[Expr], Expr]) – The norm, see NormType. Not l2.
relative (bool) – Divide by the norm of the exact solution.
coordinates (CoordinateSystem | None) – See integrate().

Return type:

float

sumfem.post_processing.integrate(expr, mesh=None, coordinates=None)[source]#

The integral of an integrand over its whole mesh, on every process.

Parameters:

expr (Expr | float) – A scalar integrand of DOLFINx functions, constants and the mesh’s coordinates, e.g. the square of an error.
mesh (Mesh | None) – The mesh, needed when the integrand does not name one, such as a number to measure the mesh’s volume.
coordinates (CoordinateSystem | None) – The coordinate system of the mesh, which weights the integral and writes its vectors, e.g. that of the domain the mesh belongs to. Cartesian by default.

Return type:

float

sumfem.post_processing.interpolate_to_line(expr, p0, p1, num_points, element)[source]#

Interpolate an expression into a function on a line segment.

The line is L(alpha) = p0 + alpha (p1 - p0) for alpha in [0, 1], meshed by create_line(). The expression is first interpolated into a space of the element on its own mesh, and then moved to the line by non-matching interpolation.

Parameters:

expr (Expr) – A UFL expression on a single mesh.
p0 (ndarray[tuple[Any, ...], dtype[floating]]) – The first point of the line.
p1 (ndarray[tuple[Any, ...], dtype[floating]]) – The last point of the line.
num_points (int) – Number of points on the line, including both ends.
element (ElementMetaData) – The element of the function on the line, and of the intermediate space on the expression’s mesh.

Returns:

The function on the line mesh.

Raises:

ValueError – If the expression is not on a single mesh.

Return type:

Function

sumfem.post_processing.norm(expr, kind=NormType.L2, mesh=None, coordinates=None)[source]#

A norm of an expression over its whole mesh, on every process.

Parameters:

expr (Expr) – A scalar, vector or tensor expression, e.g. an error.
kind (NormType | str | Callable[[Expr], Expr]) – Which norm, see NormType. For l2, expr must be a function. A callable gives a custom norm: it takes f and returns the integrand of the squared norm, e.g. lambda f: w * ufl.inner(f, f) for an L2 norm with the weight w.
mesh (Mesh | None) – See integrate().
coordinates (CoordinateSystem | None) – See integrate(). It does not weight l2.

Raises:

TypeError – For "l2" of anything but a function.
ValueError – For an unknown kind.

Return type:

float

Python API

Contents

Python API#

sumfem#

The symbolic layer#

The discrete layer#

The solvers layer#

Post-processing#