Finite element (`dolfinx::fem`)#

Finite elements#

template<std::floating_point T> class FiniteElement#

Model of a finite element.

Provides the dof layout on a reference element, and various methods for evaluating and transforming the basis.

Value shapes

An element has two value shapes, and they are not interchangeable.

A basis is tabulated on the reference cell and pushed forward to a physical cell. FiniteElement::reference_value_shape is the shape Basix tabulates in, and is the shape of the data returned by FiniteElement::tabulate. FiniteElement::value_shape is the shape of the field after the push-forward. It is the shape a user of the space sees, and the one UFL reports for a Coefficient or Argument on the space, so it is the one to use for anything user-facing, for checking that two spaces are compatible, and for sizing buffers that hold physical values.

The two differ for three independent reasons:

Blocking. A blocked element repeats a scalar base element at each dof point, e.g. a vector Lagrange space. Its value shape is whatever the caller asked for, while its reference value shape is that of the scalar base element, {}.
Quadrature elements, whose reference value shape is {}.
The map. A Piola-mapped basis is pushed forward with the Jacobian J, which has shape (gdim, tdim). The value axes that J contracts therefore have extent gdim in physical space and tdim on the reference cell. These are equal unless gdim != tdim, i.e. on a manifold. See ::compute_value_shape.

Blocking

FiniteElement::block_size is the number of dofs collocated at a dof point. For a blocked element it equals FiniteElement::value_size, except for a symmetric rank-2 tensor, which stores only its independent components: a {3, 3} symmetric element has value size 9 but block size 6.

Code that needs the number of physical components in one block must therefore call FiniteElement::physical_base_value_size rather than compute value_size() / block_size(), which is wrong for symmetric elements.

Examples

All on a mesh of triangles, so tdim == 2. RT is Raviart-Thomas and sym marks a symmetric element. The columns are FiniteElement::value_shape, FiniteElement::block_size, FiniteElement::reference_value_shape and FiniteElement::physical_base_value_size. FiniteElement::value_size and FiniteElement::reference_value_size are the products of the respective shapes and are not tabulated.

element	gdim	vshape	bs	rvshape	bvsize
P1	2	{}	1	{}	1
P1	3	{}	1	{}	1
P1, shape {2}	2	{2}	2	{}	1
P1, shape {2}	3	{2}	2	{}	1
P1, shape {5}	2	{5}	5	{}	1
P1, sym {3, 3}	2	{3, 3}	6	{}	1
RT 1	2	{2}	1	{2}	2
RT 1	3	{3}	1	{2}	3
Regge 0	2	{2, 2}	1	{2, 2}	4
Regge 0	3	{3, 3}	1	{2, 2}	9
mixed	any	throws	1	throws	throws

Note first that gdim changes nothing for P1 or P1, shape {2}: both are identity mapped, so their physical and reference value shapes agree on a manifold exactly as they do on a 2D mesh, and the blocked shape {2} stays {2} on a gdim == 3 mesh. Contrast the RT 1 and Regge 0 pairs, where the Piola push-forward does introduce gdim.

Some rows deserve further comment:

P1, shape {5}. A blocked element’s value shape is chosen by the caller and has nothing to do with gdim. A 5-vector, or a {3, 3} tensor field, on a mesh of any geometric dimension is legal, and gdim must never be substituted into such a shape.
P1, sym {3, 3}. The only case in which the block size is not the value size: the value size is 9, but only the six independent components are stored, one dof each, and the other three are recovered by symmetry.
RT 1 with gdim == 3, i.e. a triangle embedded in 3D. The contravariant Piola push-forward multiplies by a (3, 2) Jacobian, so the field has three components in physical space while Basix still tabulates two. This is the case that separates the two value shapes for a non-blocked element.
mixed. A mixed element has no value shape of its own, so every shape and size accessor throws; only FiniteElement::block_size is defined, and is 1. Extract a sub-element with ::extract_sub_element to get at its shapes. FiniteElement::is_mixed
reports this case, and is implemented as “has no reference value

shape”.

Public Types

using geometry_type = T #: Geometry type of the Mesh that the FunctionSpace is defined on.

Public Functions

FiniteElement(const basix::FiniteElement<geometry_type> &element, std::size_t gdim, const std::optional<std::vector<std::size_t>> &value_shape = std::nullopt, bool symmetric = false)#

Create a finite element from a Basix finite element.

Parameters:

element – [in] Basix finite element.
gdim – [in] Geometric dimension of the mesh the element will be used on.
value_shape – [in] Value shape for blocked element, e.g. {3} for a vector in 3D or {2, 2} for a rank-2 tensor in 2D. Can only be set for blocked scalar element. For other elements and scalar elements it should be std::nullopt.
symmetric – [in] Is the element a symmetric tensor? Should only set for 2nd-order tensor blocked elements.

FiniteElement(std::vector<BasixElementData<geometry_type>> elements, std::size_t gdim)#

Create a mixed finite element from Basix finite elements.

See FiniteElement(const std::vector<std::shared_ptr<constFiniteElement<geometry_type>>>&) for a discussion of mixed elements.

Parameters:

elements – [in] List of (Basix finite element, block size, symmetric) tuples, one for each element in the mixed element.
gdim – [in] Geometric dimension of the mesh the element will be used on, applied to every sub-element.

FiniteElement(const std::vector<std::shared_ptr<const FiniteElement<geometry_type>>> &elements)#

Create a mixed finite element from a list of finite elements.

This constructs a mixed element \(E_0 \times E_1 \times \ldots \times E_{n-1}\). The *i*th sub-element \(E_i\) can be accessed by ::extract_sub_element. Functions defined on mixed element spaces cannot be interpolated into directly. It is necessary to first extract a sub-Function (view), which can then be interpolated into.

A mixed element can be constructed from one element. In this case the FiniteElement behaves like a mixed element and cannot be interpolated into. The underlying element can be accessed using ::extract_sub_element.

Parameters:: elements – [in] Finite elements to compose the mixed element from.

FiniteElement(mesh::CellType cell_type, std::span<const geometry_type> points, std::array<std::size_t, 2> pshape, std::vector<std::size_t> value_shape = {}, bool symmetric = false)#

Create a quadrature element.

Note

A quadrature element is identity mapped, so value_shape is both its reference and its physical value shape and no geometric dimension is required.

Parameters:

cell_type – [in] Cell type.
points – [in] Quadrature points.
pshape – [in] Shape of points array.
value_shape – [in] Value shape for the element.
symmetric – [in] Is the element a symmetric tensor?

FiniteElement(const FiniteElement &element) = delete#: Copy constructor.

FiniteElement(FiniteElement &&element) = default#: Move constructor.

~FiniteElement() = default#: Destructor.

FiniteElement &operator=(const FiniteElement &element) = delete#: Copy assignment.

FiniteElement &operator=(FiniteElement &&element) = default#: Move assignment.

bool operator==(const FiniteElement &e) const#

Check if two elements are equivalent.

Two elements are equivalent if they wrap the same Basix element and have the same physical value shape and symmetry, i.e. if degrees of freedom can be copied directly from one to the other.

Note

Equality can be checked only for non-mixed elements. For a mixed element, this function will throw an exception.

Returns:: True is the two elements are the same.

bool operator!=(const FiniteElement &e) const#

Check if two elements are not equivalent.

Note

Equality can be checked only for non-mixed elements. For a mixed element, this function will raise an exception.

Returns:: True is the two elements are not the same.

mesh::CellType cell_type() const noexcept#: Cell shape that the element is defined on.

const std::string &signature() const noexcept#

String identifying the finite element.

Warning

The function is provided for convenience, but it should not be relied upon for determining the element type. Use other functions, commonly returning enums, to determine element properties.

Returns:: Element signature

int space_dimension() const noexcept#

Dimension of the finite element function space (the number of degrees-of-freedom for the element).

For ‘blocked’ elements, this function returns the dimension of the full element rather than the dimension of the base element.

Returns:: Dimension of the finite element space.

int block_size() const noexcept#

Block size of the finite element function space.

For non-blocked elements, this is always 1. For blocked elements, this is the number of DOFs collocated at each DOF point, which equals FiniteElement::value_size except for a symmetric rank-2 tensor. A symmetric rank-2 tensor stores only its independent components, so a {2, 2} symmetric element has value size 4 and block size 3, and a {3, 3} one has value size 9 and block size 6.

Returns:: Block size of the finite element space.

int value_size() const#

Value size of the finite element field in physical space.

The value size is the number of components of the finite element field once the basis has been pushed forward to a physical cell. It is the product of FiniteElement::value_shape, e.g. 1 for a scalar function, 2 for a 2D vector, 9 for a second-order tensor in 3D, etc. For blocked elements this is the value size of the full ‘blocked’ element.

Throws:: std::runtime_error – for a mixed element, which does not have a value shape.
Returns:: The value size.

std::span<const std::size_t> value_shape() const#

Value shape of the finite element field in physical space.

The value shape describes the shape of the finite element field once the basis has been pushed forward to a physical cell, e.g. {} for a scalar, {2} for a vector in 2D, {3, 3} for a rank-2 tensor in 3D, etc.

It differs from FiniteElement::reference_value_shape for blocked and quadrature elements, and for a Piola-mapped element on a manifold: Raviart-Thomas on a triangle embedded in 3D has reference value shape {2} and value shape {3}, because the push-forward contracts the reference value axis with a Jacobian of shape (gdim, tdim). See ::compute_value_shape.

Throws:: std::runtime_error – for a mixed element, which does not have a value shape.
Returns:: The value shape.

int physical_base_value_size() const#

Number of physical components in one block of the finite element field.

A blocked element repeats a scalar base element FiniteElement::block_size times, so one block of its field is a single scalar and this is 1. A non-blocked element has a single block, so this is FiniteElement::value_size.

This is the size of the push-forward of one (non-blocked) basis function, and hence the extent a buffer needs when it holds physical values one block at a time. It is the physical counterpart of FiniteElement::reference_value_size, and equals it unless the element is Piola mapped on a manifold, where it is gdim rather than tdim.

Note

This is not value_size() / block_size(): for a symmetric rank-2 tensor element that expression gives 4/3 or 9/6 rather than the correct value of 1.

Throws:: std::runtime_error – for a mixed element, which does not have a value shape.
Returns:: Number of physical components per block.

int reference_value_size() const#

Value size of the base (non-blocked) finite element field on the reference cell.

The reference value size is the product of the reference value shape, e.g. it is 1 for a scalar element, 2 for a 2D (non-blocked) vector, 9 for a (non-blocked) second-order tensor in 3D, etc. It is the number of components produced by FiniteElement::tabulate.

For blocked elements, this function returns the value size for the ‘base’ element from which the blocked element is composed.

Throws:: std::runtime_error – for a mixed element, which does not have a reference value shape.
Returns:: The value size.

std::span<const std::size_t> reference_value_shape() const#

Value shape of the base (non-blocked) finite element field on the reference cell.

This is the shape Basix tabulates in. For blocked and quadrature elements the returned shape will be {}. For other elements it is the same as FiniteElement::value_shape except on a manifold, where a Piola-mapped element has tdim reference components and gdim physical ones.

Mixed elements do not have a reference value shape.

Throws:: std::runtime_error – for a mixed element, which does not have a reference value shape.
Returns:: The value shape.

const std::vector<std::vector<std::vector<int>>> &entity_dofs() const noexcept#: Local DOFs associated with each sub-entity of the cell.

const std::vector<std::vector<std::vector<int>>> &entity_closure_dofs() const noexcept#: Local DOFs associated with the closure of each sub-entity of the cell.

bool symmetric() const#

Does the element represent a symmetric 2-tensor?

A symmetric element has a square rank-2 FiniteElement::value_shape, but stores only the independent components, so its FiniteElement::block_size is d * (d + 1) / 2 rather than d * d. See the examples in the class documentation.

void tabulate(std::span<geometry_type> values, std::span<const geometry_type> X, std::array<std::size_t, 2> shape, int order) const#

Evaluate derivatives of the basis functions up to given order at points in the reference cell.

Parameters:

values – [inout] Array that will be filled with the tabulated basis values. Must have shape (num_derivatives, num_points, num_dofs, reference_value_size) (row-major storage)
X – [in] The reference coordinates at which to evaluate the basis functions. Shape is (num_points, topological dimension) (row-major storage).
shape – [in] Shape of X.
order – [in] Number of derivatives (up to and including this order) to tabulate for.

std::pair<std::vector<geometry_type>, std::array<std::size_t, 4>> tabulate(std::span<const geometry_type> X, std::array<std::size_t, 2> shape, int order) const#

Evaluate all derivatives of the basis functions up to given order at given points in reference cell.

Parameters:

X – [in] The reference coordinates at which to evaluate the basis functions. Shape is (num_points, topological dimension) (row-major storage).
shape – [in] Shape of X.
order – [in] Number of derivatives (up to and including this order) to tabulate for.

Returns:

Basis function values and array shape (row-major storage).

int num_sub_elements() const noexcept#

Number of sub elements (for a mixed or blocked element).

Returns:: Number of sub elements.

bool is_mixed() const noexcept#

Check if element is a mixed element.

A mixed element is composed of two or more elements of different types. A blocked element, e.g. a Lagrange element with block size >= 1 is not considered mixed.

Returns:: True if element is mixed.

const std::vector<std::shared_ptr<const FiniteElement<geometry_type>>> &sub_elements() const noexcept#: Get subelements (if any).

std::shared_ptr<const FiniteElement<geometry_type>> extract_sub_element(const std::vector<int> &component) const#: Extract sub finite element for component.

const basix::FiniteElement<geometry_type> &basix_element() const#

Return underlying Basix element (if it exists).

Throws:: std::invalid_argument – if there is no Basix element.

basix::maps::type map_type() const#: Get the map type used by the element.

bool interpolation_ident() const noexcept#

Check if interpolation into the finite element space is an identity operation given the evaluation on an expression at specific points, i.e. the degree-of-freedom are equal to point evaluations. The function will return true for Lagrange elements.

Returns:: True if interpolation is an identity operation

bool map_ident() const noexcept#

Check if the push forward/pull back map from the values on reference to the values on a physical cell for this element is the identity map.

Returns:: True if the map is the identity

std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> interpolation_points() const#

Points on the reference cell at which an expression needs to be evaluated in order to interpolate the expression in the finite element space.

For Lagrange elements the points will just be the nodal positions. For other elements the points will typically be the quadrature points used to evaluate moment degrees of freedom.

Returns:: Interpolation point coordinates on the reference cell, returning the (0) coordinates data (row-major) storage and (1) the shape (num_points, tdim).

std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> interpolation_operator() const#

Interpolation operator (matrix) Pi that maps a function evaluated at the points provided by FiniteElement::interpolation_points to the element degrees of freedom, i.e. dofs = Pi f_x. See the Basix documentation for basix::FiniteElement::interpolation_matrix for how the data in f_x should be ordered.

Returns:: The interpolation operator Pi, returning the data for Pi (row-major storage) and the shape (num_dofs, num_points * value_size)

std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> create_interpolation_operator(const FiniteElement &from) const#

Create a matrix that maps degrees of freedom from one element to this element (interpolation).

Note

Does not support mixed elements.

Parameters:: from – [in] The element to interpolate from.
Returns:: Matrix operator that maps the from degrees-of-freedom to the degrees-of-freedom of this element. The (0) matrix data (row-major storage) and (1) the shape (num_dofs of this element, num_dofs of from) are returned.
Pre:: The two elements must use the same mapping between the reference and physical cells.

bool needs_dof_transformations() const noexcept#

Check if DOF transformations are needed for this element.

DOF transformations will be needed for elements which might not be continuous when two neighbouring cells disagree on the orientation of a shared sub-entity, and when this cannot be corrected for by permuting the DOF numbering in the dofmap.

For example, Raviart-Thomas elements will need DOF transformations, as the neighbouring cells may disagree on the orientation of a basis function, and this orientation cannot be corrected for by permuting the DOF numbers on each cell.

Returns:: True if DOF transformations are required.

bool depends_on_facet_orientations() const noexcept#

Check if the basis changes sign on the facets whose two cells disagree on the normal (see mesh::Topology::create_facet_orientations).

True for H(div)-conforming elements on a manifold, i.e. elements mapped by the contravariant Piola map that need DOF transformations (see T_apply(), which applies the signs), and for mixed elements with such a sub-element. Discontinuous elements mapped by the contravariant Piola map do not depend on it: with no continuity between cells to preserve, a sign changes the basis but not the space.

Returns:: True if the basis depends on the facet signs.

bool needs_dof_permutations() const noexcept#

Check if DOF permutations are needed for this element.

DOF permutations will be needed for elements which might not be continuous when two neighbouring cells disagree on the orientation of a shared subentity, and when this can be corrected for by permuting the DOF numbering in the dofmap.

For example, higher order Lagrange elements will need DOF permutations, as the arrangement of DOFs on a shared sub-entity may be different from the point of view of neighbouring cells, and this can be corrected for by permuting the DOF numbers on each cell.

Returns:: True if DOF transformations are required.

template<typename U> inline std::function<void(std::span, std::span<const std::uint32_t>, std::int32_t, int)> dof_transformation_fn(doftransform ttype, bool scalar_element = false) const#

Return a function that applies a DOF transformation operator to some data (see T_apply()).

The transformation is applied from the left-hand side, i.e.

\[ u \leftarrow T u. \]

If the transformation for the (sub)element is a permutation only, the returned function will do change the ordering for the (sub)element as it is assumed that permutations are incorporated into the degree-of-freedom map.

See the documentation for T_apply() for a description of the transformation for a single element type. This function generates a function that can apply the transformation to a mixed element.

The signature of the returned function has four arguments:

[in,out] data The data to be transformed. This data is flattened with row-major layout, shape=(num_dofs, block_size)
[in] cell_info Permutation data for the cell. The size of this is num_cells. For elements where no transformations are required, an empty span can be passed in.
[in] cell The cell number.
[in] n The block_size of the input data.

Parameters:

ttype – [in] The transformation type. Typical usage is:
- doftransform::standard Transforms basis function data from the reference element to the conforming ‘physical’ element, e.g. \(\phi = T \tilde{\phi}\).
- doftransform::transpose Transforms degree-of-freedom data from the conforming (physical) ordering to the reference ordering, e.g. \(\tilde{u} = T^{T} u\).
- doftransform::inverse: Transforms basis function data from the the conforming (physical) ordering to the reference ordering, e.g. \(\tilde{\phi} = T^{-1} \phi\).
- doftransform::inverse_transpose: Transforms degree-of-freedom data from the reference element to the conforming (physical) ordering, e.g. \(u = T^{-t} \tilde{u}\).
scalar_element – [in] Indicates whether the scalar transformations should be returned for a vector element.

template<typename U> inline std::function<void(std::span, std::span<const std::uint32_t>, std::int32_t, int)> dof_transformation_right_fn(doftransform ttype, bool scalar_element = false) const#

Return a function that applies DOF transformation to some transposed data (see T_apply_right()).

The transformation is applied from the right-hand side, i.e.

\[ u^{t} \leftarrow u^{t} T. \]

If the transformation for the (sub)element is a permutation only, the returned function will do change the ordering for the (sub)element as it is assumed that permutations are incorporated into the degree-of-freedom map.

The signature of the returned function has four arguments:

[in,out] data The data to be transformed. This data is flattened with row-major layout, shape=(num_dofs, block_size)
[in] cell_info Permutation data for the cell. The size of this is num_cells. For elements where no transformations are required, an empty span can be passed in.
[in] cell The cell number
[in] block_size The block_size of the input data

Parameters:

ttype – [in] Transformation type. See dof_transformation_fn().
scalar_element – [in] Indicate if the scalar transformations should be returned for a vector element.

template<typename U> inline void T_apply(std::span data, std::uint32_t cell_permutation, int n) const#

Transform basis functions from the reference element ordering and orientation to the globally consistent physical element ordering and orientation.

Consider that the value of a finite element function \(f_{h}\) at a point is given by

\[ f_{h} = \phi^{T} c, \]

where \(f_{h}\) has shape \(r \times 1\), \(\phi\) has shape \(d \times r\) and holds the finite element basis functions, and \(c\) has shape \(d \times 1\) and holds the degrees-of-freedom. The basis functions and degree-of-freedom are with respect to the physical element orientation. If the degrees-of-freedom on the physical element orientation are given by

\[\phi = T \tilde{\phi}, \]

where \(T\) is a \(d \times d\) matrix, it follows from \(f_{h} = \phi^{T} c = \tilde{\phi}^{T} T^{T} c\) that

\[ \tilde{c} = T^{T} c. \]

This function applies \(T\) to data. The transformation is performed in-place. The operator \(T\) is orthogonal for many elements, but not all.

This function calls the corresponding Basix function.

Parameters:

data – [inout] Data to transform. The shape is (m, n), where m is the number of dgerees-of-freedom and the storage is row-major.
cell_permutation – [in] Permutation data for the cell
n – [in] Number of columns in data.

template<typename U> inline void Tt_inv_apply(std::span data, std::uint32_t cell_permutation, int n) const#

Apply the inverse transpose of the operator applied by T_apply().

The transformation

\[ v = T^{-T} u \]

is performed in-place.

Parameters:

data – [inout] The data to be transformed. This data is flattened with row-major layout, shape=(num_dofs, block_size).
cell_permutation – [in] Permutation data for the cell.
n – [in] Block_size of the input data.

template<typename U> inline void Tt_apply(std::span data, std::uint32_t cell_permutation, int n) const#

Apply the transpose of the operator applied by T_apply().

The transformation

\[ u \leftarrow T^{T} u \]

is performed in-place.

Parameters:

data – [inout] The data to be transformed. This data is flattened with row-major layout, shape=(num_dofs, block_size).
cell_permutation – [in] Permutation data for the cell.
n – [in] The block size of the input data.

template<typename U> inline void Tinv_apply(std::span data, std::uint32_t cell_permutation, int n) const#

Apply the inverse of the operator applied by T_apply().

The transformation

\[ v = T^{-1} u \]

is performed in-place.

Parameters:

data – [inout] The data to be transformed. This data is flattened with row-major layout, shape=(num_dofs, block_size).
cell_permutation – [in] Permutation data for the cell.
n – [in] Block size of the input data.

template<typename U> inline void T_apply_right(std::span data, std::uint32_t cell_permutation, int n) const#

Right(post)-apply the operator applied by T_apply().

Computes

\[ v^{T} = u^{T} T \]

in-place.

Parameters:

data – [inout] The data to be transformed. This data is flattened with row-major layout, shape=(num_dofs, block_size).
cell_permutation – [in] Permutation data for the cell
n – [in] Block size of the input data

template<typename U> inline void Tinv_apply_right(std::span data, std::uint32_t cell_permutation, int n) const#

Right(post)-apply the inverse of the operator applied by T_apply().

Computes

\[ v^{T} = u^{T} T^{-1} \]

in-place.

Parameters:

data – [inout] Data to be transformed. This data is flattened with row-major layout, shape=(num_dofs, block_size).
cell_permutation – [in] Permutation data for the cell
n – [in] Block size of the input data

template<typename U> inline void Tt_apply_right(std::span data, std::uint32_t cell_permutation, int n) const#

Right(post)-apply the transpose of the operator applied by T_apply().

Computes

\[ v^{T} = u^{T} T^{T} \]

in-place.

Parameters:

data – [inout] Data to be transformed. The data is flattened with row-major layout, shape=(num_dofs, block_size).
cell_permutation – [in] Permutation data for the cell
n – [in] Block size of the input data.

template<typename U> inline void Tt_inv_apply_right(std::span data, std::uint32_t cell_permutation, int n) const#

Right(post)-apply the transpose inverse of the operator applied by T_apply().

Computes

\[ v^{T} = u^{T} T^{-T} \]

in-place.

Parameters:

data – [inout] Data to be transformed. This data is flattened with row-major layout, shape=(num_dofs, block_size).
cell_permutation – [in] Permutation data for the cell.
n – [in] Block size of the input data.

void permute(std::span<std::int32_t> doflist, std::uint32_t cell_permutation) const#

Permute indices associated with degree-of-freedoms on the reference element ordering to the globally consistent physical element degree-of-freedom ordering.

Given an array \(\tilde{d}\) that holds an integer associated with each degree-of-freedom and following the reference element degree-of-freedom ordering, this function computes

\[ d = P \tilde{d},\]

where \(P\) is a permutation matrix and \(d\) holds the integers in \(\tilde{d}\) but permuted to follow the globally consistent physical element degree-of-freedom ordering. The permutation is computed in-place.

Parameters:

doflist – [inout] Indices associated with the degrees-of-freedom. Size=num_dofs.
cell_permutation – [in] Permutation data for the cell.

void permute_inv(std::span<std::int32_t> doflist, std::uint32_t cell_permutation) const#

Perform the inverse of the operation applied by permute.

Given an array \(d\) that holds an integer associated with each degree-of-freedom and following the globally consistent physical element degree-of-freedom ordering, this function computes

\[ \tilde{d} = P^{T} d, \]

where \(P^{T}\) is a permutation matrix and \(\tilde{d}\) holds the integers in \(d\) but permuted to follow the reference element degree-of-freedom ordering. The permutation is computed in-place.

Parameters:

doflist – [inout] Indices associated with the degrees-of-freedom. Size=num_dofs.
cell_permutation – [in] Permutation data for the cell.

std::function<void(std::span<std::int32_t>, std::uint32_t)> dof_permutation_fn(bool inverse = false, bool scalar_element = false) const#

Return a function that applies a degree-of-freedom permutation to some data.

The returned function can apply permute to mixed-elements.

The signature of the returned function has three arguments:

[in,out] doflist The numbers of the DOFs, a span of length num_dofs
[in] cell_permutation Permutation data for the cell
[in] block_size The block_size of the input data

Parameters:

inverse – [in] Indicates if the inverse transformation should be returned.
scalar_element – [in] Indicates is the scalar transformations should be returned for a vector element.

enum class dolfinx::fem::doftransform : std::uint8_t#

DOF transformation type.

Values:

enumerator standard#: Standard.

enumerator transpose#: Transpose.

enumerator inverse#: Inverse.

enumerator inverse_transpose#: Transpose inverse.

Function spaces and functions#

Finite element functions, expressions and constants

Function spaces#

template<std::floating_point T> class FunctionSpace#

This class represents a finite element function space defined by a mesh, a finite element, and a local-to-global map of the degrees-of-freedom.

Template Parameters:: T – The floating point (real) type of the mesh geometry and the finite element basis.

Public Types

using geometry_type = T #: Geometry type of the Mesh that the FunctionSpace is defined on.

Public Functions

inline FunctionSpace(std::shared_ptr<const mesh::Mesh<geometry_type>> mesh, std::shared_ptr<const FiniteElement<geometry_type>> element, std::shared_ptr<const DofMap> dofmap)#

Create function space for given mesh, element and degree-of-freedom map.

Parameters:

mesh – [in] Mesh that the space is defined on.
element – [in] Finite element for the space.
dofmap – [in] Degree-of-freedom map for the space.

Throws:

std::runtime_error – If an element depends on the facet signs (see FiniteElement::depends_on_facet_orientations) and they have not been computed on the mesh (see mesh::Topology::has_facet_orientations).

inline FunctionSpace(std::shared_ptr<const mesh::Mesh<geometry_type>> mesh, std::vector<std::shared_ptr<const FiniteElement<geometry_type>>> elements, std::vector<std::shared_ptr<const DofMap>> dofmaps)#

Create function space for given mesh, elements and degree-of-freedom maps.

Parameters:

mesh – [in] Mesh that the space is defined on.
elements – [in] Finite elements for the space, one for each cell type. The elements must be ordered to be consistent with mesh::topology::cell_types.
dofmaps – [in] Degree-of-freedom maps for the space, one for each element. The dofmaps must be ordered in the same way as the elements.

Throws:

std::runtime_error – If an element depends on the facet signs (see FiniteElement::depends_on_facet_orientations) and they have not been computed on the mesh (see mesh::Topology::has_facet_orientations).

FunctionSpace(FunctionSpace &&V) = default#: Move constructor.

~FunctionSpace() = default#: Destructor.

FunctionSpace &operator=(FunctionSpace &&V) = default#: Move assignment operator.

inline FunctionSpace sub(const std::vector<int> &component) const#

Create a subspace (view) for a specific component.

Note

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

Parameters:: component – [in] Subspace component.
Returns:: A subspace.

inline bool contains(const FunctionSpace &V) const#

Check whether V is subspace of this, or this itself.

Parameters:: V – [in] The space to be tested for inclusion
Returns:: True if V is contained in or is equal to this FunctionSpace

inline std::pair<FunctionSpace, std::vector<std::vector<std::int32_t>>> collapse() const#

Collapse a subspace and return a new function space and a map from new to old dofs

Returns:: The new function space and a map from new to old dofs

inline std::vector<int> component() const#

Get the component with respect to the root superspace.

Returns:: The component with respect to the root superspace, i.e. W.sub(1).sub(0) == [1, 0].

inline bool symmetric() const#: Indicate whether this function space represents a symmetric 2-tensor.

inline std::vector<geometry_type> tabulate_dof_coordinates(bool transpose) const#

Tabulate the physical coordinates of all dofs on this process.

Todo:: Remove - see function in interpolate.h

Parameters:: transpose – [in] If false the returned data has shape (num_points, 3), otherwise it is transposed and has shape (3, num_points).
Returns:: The dof coordinates [([x0, y0, z0], [x1, y1, z1], ...) if transpose is false, and otherwise the returned data is transposed. Storage is row-major.

inline std::shared_ptr<const mesh::Mesh<geometry_type>> mesh() const#: The mesh.

inline std::shared_ptr<const FiniteElement<geometry_type>> element() const#: The finite element.

inline std::shared_ptr<const FiniteElement<geometry_type>> elements(int cell_type_idx) const#: The finite elements.

inline std::shared_ptr<const DofMap> dofmap() const#: The dofmap.

inline const std::vector<std::shared_ptr<const DofMap>> &dofmaps() const#: The dofmaps.

Functions#

template<dolfinx::scalar T, std::floating_point U = dolfinx::scalar_value_t<T>> class Function#

This class represents a function \( u_h \) in a finite element function space \( V_h \), given by

\[ u_h = \sum_{i=1}^{n} U_i \phi_i, \]

where \( \{\phi_i\}_{i=1}^{n} \) is a basis for \( V_h \), and \( U \) is a vector of expansion coefficients for \( u_h \).

Template Parameters:

T – The function scalar type.
U – The mesh geometry scalar type.

Public Types

using value_type = T #: Field type for the Function, e.g. double, std::complex<float>, etc.

using geometry_type = U #: Geometry type of the Mesh that the Function is defined on.

Public Functions

inline explicit Function(std::shared_ptr<const FunctionSpace<geometry_type>> V, std::string name = "u")#

Create function on given function space.

Parameters:

V – [in] The function space
name – [in] The function name.

inline Function(std::shared_ptr<const FunctionSpace<geometry_type>> V, std::shared_ptr<la::Vector<value_type>> x, std::string name = "u")#

Create function on given function space with a given vector.

Warning

This constructor is intended for internal library use only.

Parameters:

V – [in] The function space.
x – [in] The vector.
name – [in] The function name.

Function(Function &&v) = default#: Move constructor.

~Function() = default#: Destructor.

Function &operator=(Function &&v) = default#: Move assignment.

inline Function sub(int i) const#

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

Parameters:: i – [in] Index of subfunction
Returns:: The sub-function

inline Function collapse() const#

Collapse a subfunction (view into a Function) to a stand-alone Function.

Returns:: New collapsed Function.

inline std::shared_ptr<const FunctionSpace<geometry_type>> function_space() const#

Access the function space.

Returns:: The function space.

inline std::shared_ptr<const la::Vector<value_type>> x() const#: Underlying vector (const version).

inline std::shared_ptr<la::Vector<value_type>> x()#: Underlying vector.

inline void interpolate(const std::function<std::pair<std::vector<value_type>, std::vector<std::size_t>>(md::mdspan<const geometry_type, md::extents<std::size_t, 3, md::dynamic_extent>>)> &f, mesh::CellRange auto &&cells)#

Interpolate an expression f(x) over a set of cells.

Parameters:

f – [in] Expression function to be interpolated.
cells – [in] Cells to interpolate on.

inline void interpolate(const std::function<std::pair<std::vector<value_type>, std::vector<std::size_t>>(md::mdspan<const geometry_type, md::extents<std::size_t, 3, md::dynamic_extent>>)> &f)#

Interpolate an expression f(x) on the whole domain.

Parameters:: f – [in] Expression to be interpolated.

inline void interpolate(const Function<value_type, geometry_type> &u0, mesh::CellRange auto &&cells0, mesh::CellRange auto &&cells1)#

Interpolate a Function over a subset of cells.

The Function being interpolated from and the Function being interpolated into can be defined on different sub-meshes, i.e. views into a subset a cells.

Parameters:

u0 – [in] Function to be interpolated.
cells0 – [in] Cells to interpolate from. These are the indices of the cells in the mesh associated with u0.
cells1 – [in] Cell indices associated with the mesh of this that will be interpolated to. If cells0[i] is the index of a cell in the mesh associated with u0, then cells1[i] is the index of the same cell but in the mesh associated with this.

Pre:

cells0 and cells1 must have the same length.

inline void interpolate(const Function<value_type, geometry_type> &u, mesh::CellRange auto &&cells)#

Interpolate a Function over a subset of cells.

The Functions must be defined on the same mesh.

Parameters:

u – [in] Function to be interpolated.
cells – [in] Cells to interpolate from.

inline void interpolate(const Function<value_type, geometry_type> &u)#

Interpolate a Function over all cells.

The Functions must be defined on the same mesh.

Parameters:: u – [in] Function to be interpolated.

inline void interpolate(const Function<value_type, geometry_type> &u, mesh::CellRange auto &&cells, double tol, int maxit, const geometry::PointOwnershipData &interpolation_data)#

Interpolate a Function defined on a different mesh.

Parameters:

u – [in] Function to be interpolated.
cells – [in] Cells in the mesh associated with this to interpolate into.
tol – [in] Tolerance for convergence in Newton method for non-affine pullbacks. If the mesh geometry is affine this argument is ignored.
maxit – [in] Maximum number of Newton iterations in non-affine pull-back. If the mesh geometry is affine this argument is ignored.
interpolation_data – [in] Data required for associating the interpolation points of this with cells in u. Can be computed with fem::create_interpolation_data.

inline void eval(std::span<const geometry_type> x, std::array<std::size_t, 2> xshape, mesh::CellRange auto &&cells, std::span<value_type> u, std::array<std::size_t, 2> ushape, double tol, int maxit) const#

Evaluate the Function at points.

Parameters:

x – [in] The coordinates of the points. It has shape (num_points, 3) and storage is row-major.
xshape – [in] Shape of x.
cells – [in] Cell indices such that cells[i] is the index of the cell that contains the point x(i). Negative cell indices can be passed, in which case the corresponding point is ignored.
u – [out] Values at the points. Values are not computed for points with a negative cell index. This argument must be passed with the correct size. Storage is row-major.
ushape – [in] Shape of u.
tol – [in] Tolerance for convergence in Newton method for non-affine pullbacks. If the mesh geometry is affine this argument is ignored.
maxit – [in] Maximum number of Newton iterations in non-affine pull-back. If the mesh geometry is affine this argument is ignored.

Public Members

std::string name#: Name.

Constants#

template<dolfinx::scalar T> class Constant#

Constant (in space) value which can be attached to a Form.

Constants may be scalar (rank 0), vector (rank 1), or tensor-valued.

Template Parameters:: T – Scalar type of the Constant.

Public Types

using value_type = T #: Field type.

Public Functions

inline explicit Constant(value_type c)#

Create a rank-0 (scalar-valued) constant.

Parameters:: c – [in] Value of the constant.

inline explicit Constant(std::span<const value_type> c)#

Create a rank-1 (vector-valued) constant.

Parameters:: c – [in] Value of the constant.

inline Constant(std::span<const value_type> c, std::span<const std::size_t> shape)#

Create a rank-d constant.

Parameters:

c – [in] Value of the Constant (row-majors storage)
shape – [in] Shape of the Constant

Public Members

std::vector<value_type> value#: Values, stored as a row-major flattened array.

std::vector<std::size_t> shape#: Shape.

Forms#

template<dolfinx::scalar T, std::floating_point U = dolfinx::scalar_value_t<T>> class Form#

A representation of finite element variational forms.

A note on the order of trial and test spaces: FEniCS numbers argument spaces starting with the leading dimension of the corresponding tensor (matrix). In other words, the test space is numbered 0 and the trial space is numbered 1. However, in order to have a notation that agrees with most existing finite element literature, in particular

\[ a = a(u, v) \]

the spaces are numbered from right to left

\[ a: V_1 \times V_0 \rightarrow \mathbb{R} \]

This is reflected in the ordering of the spaces that should be supplied to generated subclasses. In particular, when a bilinear form is initialized, it should be initialized as a(V_1, V_0) = ..., where V_1 is the trial space and V_0 is the test space. However, when a form is initialized by a list of argument spaces (the variable function_spaces in the constructors below), the list of spaces should start with space number 0 (the test space) and then space number 1 (the trial space).

Template Parameters:

T – Scalar type in the form.
U – Float (real) type used for the finite element and geometry.

Public Types

using scalar_type = T #: Scalar type.

using geometry_type = U #: Geometry type.

Public Functions

template<typename X> inline Form(const std::vector<std::shared_ptr<const FunctionSpace<geometry_type>>> &V, X &&integrals, std::shared_ptr<const mesh::Mesh<geometry_type>> mesh, const std::vector<std::shared_ptr<const Function<scalar_type, geometry_type>>> &coefficients, const std::vector<std::shared_ptr<const Constant<scalar_type>>> &constants, bool needs_facet_permutations, const std::vector<std::reference_wrapper<const mesh::EntityMap>> &entity_maps)#

Create a finite element form.

Note

User applications will normally call a factory function rather using this interface directly.

Note

For the single domain case, pass an empty entity_maps.

Parameters:

V – [in] Function spaces for the form arguments, e.g. test and trial function spaces.
integrals – [in] Integrals in the form, where integrals[{type, i, kernel_index}] gives the integral_data of type type at position i in the flattened, sorted-by-subdomain-id list of integrals, for kernel kernel_index. The subdomain ids can contain duplicate entries referring to different kernels over the same subdomain.
coefficients – [in] Coefficients in the form.
constants – [in] Constants in the form.
mesh – [in] Mesh of the domain to integrate over (the ‘integration domain’).
needs_facet_permutations – [in] Set to true if any of the integration kernels require cell permutation data.
entity_maps – [in] A list of EntityMaps. For every mesh other than mesh on which a trial function, test function, or coefficient is defined, entity_maps must contain an EntityMap relating that mesh and mesh.

Form(Form &&form) = default#: Move constructor

Note

Valid because ::_integrals is a std::map, whose elements keep a stable address across a move, so the std::spans cached in ::_edata and ::_cdata remain valid after the move.

~Form() = default#: Destructor.

Form &operator=(Form &&form) = default#: Move assignment

Note

Valid for the same reason as the move constructor: move assigning a std::map transfers its nodes, so the entity vectors aliased by the std::spans cached in ::_edata and ::_cdata keep their addresses.

inline int rank() const#

Rank of the form.

bilinear form = 2, linear form = 1, functional = 0, etc.

Returns:: The rank of the form.

inline std::shared_ptr<const mesh::Mesh<geometry_type>> mesh() const#

Common mesh for the form (the ‘integration domain’).

Returns:: The integration domain mesh.

inline const std::vector<std::shared_ptr<const FunctionSpace<geometry_type>>> &function_spaces() const#

Function spaces for all arguments.

Returns:: Function spaces.

inline std::function<void(scalar_type*, const scalar_type*, const scalar_type*, const geometry_type*, const int*, const std::uint8_t*, void*)> kernel(IntegralType type, int idx, int kernel_idx) const#

Get the kernel function for an integral.

Parameters:

type – [in] Integral type.
idx – [in] Integral index in the flattened list of integral kernels (see ::domain).
kernel_idx – [in] Index of the kernel (we may have multiple kernels for a given idx in mixed-topology meshes).

Returns:

Function to call for tabulate_tensor.

inline std::set<IntegralType> integral_types() const#

Get types of integrals in the form.

Returns:: Integrals types.

inline std::vector<int> active_coeffs(IntegralType type, int idx) const#

Indices of coefficients that are active for a given integral (kernel).

A form is split into multiple integrals (kernels) and each integral might contain only a subset of all coefficients in the form. This function returns an indicator array for a given integral kernel that signifies which coefficients are present.

Parameters:

type – [in] Integral type.
idx – [in] Integral index in the flattened list of integral kernels (see ::domain).

Returns:

Indices of the coefficients that are active (present) in the given integral.

inline int num_integrals(IntegralType type, int kernel_idx) const#

Get number of integrals (kernels) for a given integral type and kernel index.

For a form containing two integrals integral_a and integral_b with subdomain ids (1, 4) and (3, 4, 5) respectively, the integrals are stored as a flattened list, sorted by subdomain id:

auto form_integrals = {integral_a, integral_b, integral_a,
                       integral_b, integral_b};
auto form_integral_ids = {1, 3, 4, 4, 5};

Parameters:

type – [in] Integral type.
kernel_idx – [in] Index of the kernel (we may have multiple kernels for a integral type in mixed-topology meshes).

Returns:

Number of integrals (kernels) of the given type and kernel index.

inline std::span<const std::int32_t> domain(IntegralType type, int idx, int kernel_idx) const#

Mesh entity indices to integrate over for a given integral (kernel).

These are the entities in the mesh returned by ::mesh that are integrated over by a given integral (kernel).

For IntegralType::cell, returns a list of cell indices.
For IntegralType::exterior_facet, returns a list with shape (num_facets, 2) (row-major storage), where for row i, [i, 0] is the cell index and [i, 1] is the local facet index relative to the cell.
For IntegralType::interior_facet, returns a list with shape (num_facets, 4) (row-major storage), where for row i, [i, 0] is the index of one attached cell, [i, 1] is the local facet index relative to that cell, [i, 2] is the index of the other attached cell, and [i, 3] is the local facet index relative to that cell.

Parameters:

type – [in] Integral type.
idx – [in] Integral index in the flattened list of integral kernels. For a form containing two integrals integral_a and integral_b with subdomain ids (1, 4) and (3, 4, 5) respectively, the integrals are stored as a flattened list, sorted by subdomain id:
```
auto form_integrals = {integral_a, integral_b, integral_a,
                       integral_b, integral_b};
auto form_integral_ids = {1, 3, 4, 4, 5};
```
kernel_idx – [in] Index of the kernel within the domain (we may have multiple kernels for a given id in mixed-topology meshes).

Returns:

Entity indices, with respect to the mesh::Mesh returned by ::mesh, to integrate over.

inline std::span<const std::int32_t> domain_arg(IntegralType type, int rank, int idx, int kernel_idx) const#

Argument function mesh integration entity indices.

Integration can be performed over cells/facets involving functions that are defined on different meshes but which share common cells, i.e. meshes can be ‘views’ into a common mesh. Meshes can share some cells but a common cell will have a different index in each mesh::Mesh. Consider:

auto mesh = this->mesh();
auto entities = this->domain(type, idx, kernel_idx);
auto entities0 = this->domain_arg(type, rank, idx, kernel_idx);

Assembly is performed over entities, where entities[i] is an entity index (e.g., cell index) in mesh. entities0 holds the corresponding entity indices but in the mesh associated with the argument function (test/trial function) space. entities[i] and entities0[i] point to the same mesh entity, but with respect to different mesh views. In some cases, such as when integrating over the interface between two domains that do not overlap, an entity may exist in one domain but not another. In this case, the entity is marked with -1.

Parameters:

type – [in] Integral type.
rank – [in] Argument index, e.g. 0 for the test function space, 1 for the trial function space.
idx – [in] Integral identifier.
kernel_idx – [in] Index of the kernel (we may have multiple kernels for a given id in mixed-topology meshes).

Returns:

Entity indices in the argument function space mesh that is integrated over.

For cell integrals it has shape (num_cells,).
For exterior/interior facet integrals, it has shape (num_facets, 2) (row-major storage), where [i, 0] is the index of a cell and [i, 1] is the local index of the facet relative to the cell.

inline std::span<const std::int32_t> domain_coeff(IntegralType type, int idx, int c) const#

Coefficient function mesh integration entity indices.

This method is equivalent to ::domain_arg, but returns mesh entity indices for coefficient Functions.

Parameters:

type – [in] Integral type.
idx – [in] Integral identifier.
c – [in] Coefficient index.

Returns:

Entity indices in the coefficient function space mesh that is integrated over.

For cell integrals it has shape (num_cells,).
For exterior/interior facet integrals, it has shape (num_facets, 2) (row-major storage), where [i, 0] is the index of a cell and [i, 1] is the local index of the facet relative to the cell.

inline const std::vector<std::shared_ptr<const Function<scalar_type, geometry_type>>> &coefficients() const#

Access coefficients.

Returns:: Coefficients in the form.

inline bool needs_facet_permutations() const#

Get bool indicating whether permutation data needs to be passed into these integrals.

Returns:: True if cell permutation data is required

inline std::vector<int> coefficient_offsets() const#

Offset for each coefficient expansion array on a cell.

Used to pack data for multiple coefficients in a flat array. The last entry is the size required to store all coefficients.

Returns:: Coefficient offsets.

inline const std::vector<std::shared_ptr<const Constant<scalar_type>>> &constants() const#

Access constants.

Returns:: Constants in the form.

enum class dolfinx::fem::IntegralType : std::int8_t#

Type of integral.

Values:

enumerator cell#: Cell.

enumerator exterior_facet#: Exterior facet.

enumerator interior_facet#: Interior facet.

enumerator vertex#: Vertex.

enumerator ridge#: Ridge.

Dirichlet boundary conditions#

template<dolfinx::scalar T, std::floating_point U = dolfinx::scalar_value_t<T>> class DirichletBC#

Object for setting (strong) Dirichlet boundary conditions

\[u = g \ \text{on} \ G,\]

where \(u\) is the solution to be computed, \(g\) is a function and \(G\) is a sub domain of the mesh.

A DirichletBC is specified by the function \(g\), the function space (trial space) and degrees of freedom to which the boundary condition applies.

Public Functions

template<typename S, typename X> inline DirichletBC(const S &g, X &&dofs, std::shared_ptr<const FunctionSpace> V)#

Create a representation of a Dirichlet boundary condition constrained by a scalar- or vector-valued constant.

Note

Can be used only with point-evaluation elements.

Note

The indices in dofs are for blocks, e.g. a block index corresponds to 3 degrees-of-freedom if the dofmap associated with g has block size 3.

Note

The size of g must be equal to the block size if V. Use the Function version if this is not the case, e.g. for some mixed spaces.

Parameters:

g – [in] The boundary condition value (T or convertible to std::span<const T>)
dofs – [in] Degree-of-freedom block indices to be constrained. The indices must be strictly increasing.
V – [in] The function space to be constrained

Pre:

dofs must be strictly increasing, i.e. sorted and free of duplicates. Checked in Debug builds only.

template<typename X> inline DirichletBC(std::shared_ptr<const Constant<T>> g, X &&dofs, std::shared_ptr<const FunctionSpace> V)#

Create a representation of a Dirichlet boundary condition constrained by a fem::Constant.

Note

Can be used only with point-evaluation elements.

Note

The indices in dofs are for blocks, e.g. a block index corresponds to 3 degrees-of-freedom if the dofmap associated with g has block size 3.

Note

The size of g must be equal to the block size if V. Use the Function version if this is not the case, e.g. for some mixed spaces.

Parameters:

g – [in] The boundary condition value.
dofs – [in] Degree-of-freedom block indices to be constrained.
V – [in] The function space to be constrained

Pre:

dofs must be strictly increasing, i.e. sorted and free of duplicates. Checked in Debug builds only.

template<typename X> inline DirichletBC(std::shared_ptr<const Function<T, U>> g, X &&dofs)#

Create a representation of a Dirichlet boundary condition where the space being constrained is the same as the function that defines the constraint Function, i.e. share the same fem::FunctionSpace.

Note

The indices in dofs are for blocks, e.g. a block index corresponds to 3 degrees-of-freedom if the dofmap associated with g has block size 3.

Parameters:

g – [in] The boundary condition value.
dofs – [in] Degree-of-freedom block indices to be constrained.

Pre:

dofs must be strictly increasing, i.e. sorted and free of duplicates. Checked in Debug builds only.

template<typename X> inline DirichletBC(std::shared_ptr<const Function<T, U>> g, X &&V_g_dofs, std::shared_ptr<const FunctionSpace> V)#

Create a representation of a Dirichlet boundary condition where the space being constrained and the function that defines the constraint values do not share the same fem::FunctionSpace.

A typical example is when applying a constraint on a subspace. The (sub)space and the constrain function must have the same finite element.

Note

The indices in dofs are unrolled and not for blocks.

Parameters:

g – [in] The boundary condition value
V_g_dofs – [in] Two arrays of degree-of-freedom indices (std::array<std::vector<std::int32_t>, 2>). First array are indices in the space where boundary condition is applied (V), second array are indices in the space of the boundary condition value function g. The dof indices are unrolled, i.e. are not by dof block.
V – [in] The function (sub)space on which the boundary condition is applied

Pre:

The two degree-of-freedom arrays in V_g_dofs must be sorted by the indices in the first array, which must itself be free of duplicates. Checked in Debug builds only.

DirichletBC(const DirichletBC &bc) = default#

Copy constructor

Parameters:: bc – [in] The object to be copied

DirichletBC(DirichletBC &&bc) = default#

Move constructor

Parameters:: bc – [in] The object to be moved

~DirichletBC() = default#: Destructor.

DirichletBC &operator=(const DirichletBC &bc) = default#

Assignment operator

Parameters:: bc – [in] Another DirichletBC object

DirichletBC &operator=(DirichletBC &&bc) = default#: Move assignment operator.

inline std::shared_ptr<const FunctionSpace> function_space() const#

The function space to which boundary conditions are applied

Returns:: The function space

inline std::variant<std::shared_ptr<const Function<T, U>>, std::shared_ptr<const Constant<T>>> value() const#

Return boundary value function g

Returns:: The boundary values Function

inline std::pair<std::span<const std::int32_t>, std::int32_t> dof_indices() const#

Dof indices constrained by this condition, and the position of the first ghost index.

Indices are process-local and unrolled (block size expanded), into the dofmap of function_space(). For a condition on a sub-space they are indices in the parent space, so they apply directly to arrays over that space.

The array is strictly increasing, following from the dofs constructor precondition that the owned/ghost split relies on. That precondition is checked in Debug builds only.

Note

The returned span is valid for the lifetime of this object.

Returns:: Dof indices, and the position pos of the first non-owned (ghost) index. Entries dofs[:pos] are owned and entries dofs[pos:] are ghosts.

inline void set(std::span<T> x, std::optional<std::span<const T>> x0, T alpha = 1) const#

Set entries in an array that are constrained by Dirichlet boundary conditions.

Entries in x that are constrained by a Dirichlet boundary conditions are set to alpha * (x_bc - x0), where x_bc is the (interpolated) boundary condition value.

For elements with point-wise evaluated degrees-of-freedom, e.g. Lagrange elements, x_bc is the value of the boundary condition at the degree-of-freedom. For elements with moment degrees-of-freedom, x_bc is the value of the boundary condition interpolated into the finite element space.

x may hold only owned entries or also include ghosts (entries available on the calling rank but owned by another rank); a constrained degree-of-freedom beyond the end of x is silently skipped, so an owned-only x sets only owned entries.

Note

If x0 is provided, it must be at least as long as x (checked only in Debug builds).

Parameters:

x – [inout] Array to modify for Dirichlet boundary conditions. May include ghost entries.
x0 – [in] Optional array used in computing the value to set. If not provided it is treated as zero.
alpha – [in] Scaling to apply.

inline void mark_dofs(std::span<std::int8_t> markers) const#

Set markers[i] = true if dof i has a boundary condition applied.

Value of markers[i] is not changed otherwise.

Parameters:: markers – [inout] Entry makers[i] is set to true if dof i in V0 had a boundary condition applied, i.e. dofs which are fixed by a boundary condition. Other entries in markers are left unchanged.

Degree-of-freedom maps#

graph::AdjacencyList<std::int32_t> dolfinx::fem::transpose_dofmap(md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>> dofmap, std::int32_t num_cells)#

Create an adjacency list that maps a global index (process-wise) to the ‘unassembled’ cell-wise contributions.

It is built from the usual (cell, local index) -> global index dof map. An ‘unassembled’ vector is the stacked cell contributions, ordered by cell index. If the usual dof map is:

Cell: 0 1 2 3Global index: [ [0, 3, 5], [3, 2, 4], [4, 3, 2], [2, 1, 0]]

the ‘transpose’ dof map will be:

Global index: 0 1 2 3 4 5Unassembled index: [ [0, 11], [10], [4, 8, 9], [1, 3, 7], [5, 6], [2] ]

Parameters:

dofmap – [in] The standard dof map that for each cell (node) gives the global (process-wise) index of each local (cell-wise) index.
num_cells – [in] The number of cells (nodes) in dofmap to consider. The first num_cells are used. This is argument is typically used to exclude ghost cell contributions.

Returns:

Map from global (process-wise) index to positions in an unaassembled array. The links for each node are sorted.

class DofMap#

Degree-of-freedom map.

This class handles the mapping of degrees of freedom. It builds a dof map based on an ElementDofLayout on a specific mesh topology. It will reorder the dofs when running in parallel. Sub-dofmaps, both views and copies, are supported.

Public Functions

template<typename E, typename U> inline DofMap(E &&element, std::shared_ptr<const common::IndexMap> index_map, int index_map_bs, U &&dofmap, int bs)#

Create a DofMap from the layout of dofs on a reference element, an IndexMap defining the distribution of dofs across processes and a vector of indices.

Parameters:

element – [in] The layout of the degrees of freedom on an element
index_map – [in] The map describing the parallel distribution of the degrees of freedom.
index_map_bs – [in] The block size associated with the index_map.
dofmap – [in] Adjacency list with the degrees-of-freedom for each cell.
bs – [in] The block size of the dofmap.

DofMap(DofMap &&dofmap) = default#: Move constructor.

DofMap &operator=(DofMap &&dofmap) = default#: Move assignment.

bool operator==(const DofMap &map) const#

Equality operator.

Returns:: Returns true if the data for the two dofmaps are equal

inline std::span<const std::int32_t> cell_dofs(std::int32_t c) const#

Local-to-global mapping of dofs on a cell.

Parameters:: c – [in] The cell index
Returns:: Local-global dof map for the cell (using process-local indices)

int bs() const noexcept#: Return the block size for the dofmap.

DofMap extract_sub_dofmap(std::span<const int> component) const#

Extract subdofmap component.

Parameters:: component – [in] The component indices
Returns:: The dofmap for the component

std::pair<DofMap, std::vector<std::int32_t>> collapse(MPI_Comm comm, const mesh::Topology &topology, std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)> &&reorder_fn = nullptr) const#

Create a “collapsed” dofmap (collapses a sub-dofmap).

Parameters:

comm – [in] MPI Communicator
topology – [in] Mesh topology that the dofmap is defined on
reorder_fn – [in] Graph re-ordering function to apply to the dof data

Returns:

The collapsed dofmap

md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>> map() const#

Get dofmap data.

Returns:: The adjacency list with dof indices for each cell

inline const ElementDofLayout &element_dof_layout() const#: Layout of dofs on an element.

int index_map_bs() const#: Block size associated with the index_map.

Public Members

std::shared_ptr<const common::IndexMap> index_map#: Index map that describes the parallel distribution of the dofmap.

class ElementDofLayout#

The class represents the degree-of-freedom (dofs) for an element. Dofs are associated with a mesh entity. This class also handles sub-space dofs, which are views into the parent dofs.

Public Functions

ElementDofLayout(int block_size, const std::vector<std::vector<std::vector<int>>> &entity_dofs, const std::vector<std::vector<std::vector<int>>> &entity_closure_dofs, const std::vector<int> &parent_map, const std::vector<ElementDofLayout> &sub_layouts)#

Constructor

Parameters:

block_size – [in] The number of dofs co-located at each point.
entity_dofs – [in] The dofs on each entity, in the format: entity_dofs[entity_dim][entity_number] = [dof0, dof1, …]
entity_closure_dofs – [in] The dofs on the closure of each entity, in the format: entity_closure_dofs[entity_dim][entity_number] = [dof0, dof1, …]
parent_map – [in] TODO
sub_layouts – [in] (std::vector<ElementDofLayout>)TODO

ElementDofLayout copy() const#: Copy the DOF layout, discarding any parent information.

ElementDofLayout(const ElementDofLayout &dofmap) = default#: Copy constructor.

ElementDofLayout(ElementDofLayout &&dofmap) = default#: Move constructor.

~ElementDofLayout() = default#: Destructor.

ElementDofLayout &operator=(const ElementDofLayout &dofmap) = default#: Copy assignment.

ElementDofLayout &operator=(ElementDofLayout &&dofmap) = default#: Move assignment.

bool operator==(const ElementDofLayout &layout) const#

Equality operator

Note

The block sizes of the layouts are not compared

Returns:: Returns true if the layout data is the same. Sub- and parent dofmap data is not compared.

int num_dofs() const#

Return the dimension of the local finite element function space on a cell (number of dofs on element)

Returns:: Dimension of the local finite element function space.

const std::vector<int> &entity_dofs(int dim, int entity_index) const#

Local-local mapping of dofs on entity of cell

Parameters:

dim – [in] The entity dimension
entity_index – [in] The local entity index on the cell

Returns:

Cell-local degree-of-freedom indices

const std::vector<int> &entity_closure_dofs(int dim, int entity_index) const#

Local-local closure dofs on entity of cell

Parameters:

dim – [in] The entity dimension
entity_index – [in] The local entity index on the cell

Returns:

Cell-local degree-of-freedom indices

const std::vector<std::vector<std::vector<int>>> &entity_dofs_all() const#: Direct access to all entity dofs (dof = _entity_dofs[dim][entity][i]).

const std::vector<std::vector<std::vector<int>>> &entity_closure_dofs_all() const#: Direct access to all entity closure dofs (dof = _entity_dofs[dim][entity][i])

int num_sub_dofmaps() const#: Get number of sub-dofmaps.

const ElementDofLayout &sub_layout(std::span<const int> component) const#: Get sub-dofmap given by list of components, one for each level.

std::vector<int> sub_view(std::span<const int> component) const#: Get view for a sub-layout, defined by the component list (as for sub_layour()), into this dofmap. I.e., the dofs in this dofmap that are the sub-dofs.

int block_size() const#: Block size.

bool is_view() const#

True iff dof map is a view into another map

Returns:: bool True if the dof map is a sub-dof map (a view into another map).

Assembly#

Functions

template<dolfinx::scalar T> std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>> make_coefficients_span(const std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>> &coeffs)#: Create a map of std::spans from a map of std::vectors.

template<dolfinx::scalar T, std::floating_point U> T assemble_scalar(const Form<T, U> &M, std::span<const T> constants, const std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>> &coefficients)#

Assemble functional into scalar.

The caller supplies the form constants and coefficients for this version, which has efficiency benefits if the data can be re-used for multiple calls.

Note

Caller is responsible for accumulation across processes.

Parameters:

M – [in] The form (functional) to assemble
constants – [in] The constants that appear in M
coefficients – [in] The coefficients that appear in M

Returns:

The contribution to the form (functional) from the local process

template<dolfinx::scalar T, std::floating_point U> T assemble_scalar(const Form<T, U> &M)#

Assemble functional into scalar.

Note

Caller is responsible for accumulation across processes.

Parameters:: M – [in] The form (functional) to assemble.
Returns:: The contribution to the form (functional) from the local process.

template<typename V, std::floating_point U, dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type> void assemble_vector(V &&b, const Form<T, U> &L, std::span<const T> constants, const std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>> &coefficients)#

Assemble linear form into a vector.

The caller supplies the form constants and coefficients for this version, which has efficiency benefits if the data can be re-used for multiple calls.

Parameters:

b – [inout] The vector to be assembled. It will not be zeroed before assembly.
L – [in] The linear forms to assemble into b.
constants – [in] The constants that appear in L.
coefficients – [in] The coefficients that appear in L.

template<typename V, std::floating_point U, dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type> void assemble_vector(V &&b, const Form<T, U> &L)#

Assemble linear form into a vector.

Parameters:

b – [inout] Vector to be assembled. It will not be zeroed before assembly.
L – [in] Linear forms to assemble into b.

template<dolfinx::scalar T, std::floating_point U> std::vector<std::int8_t> bc_dof_markers(const FunctionSpace &V, const std::vector<std::reference_wrapper<const DirichletBC<T, U>>> &bcs)#

Mark the dofs of V (owned and ghost, unrolled) constrained by the boundary conditions in bcs that are defined on V or a subspace of it.

Returns:: Dof markers, or an empty array if no boundary condition applies.

template<dolfinx::scalar T, std::floating_point U> std::array<std::span<const std::int8_t>, 2> bc_dof_markers_pair(const Form<T, U> &a, const std::vector<std::reference_wrapper<const DirichletBC<T, U>>> &bcs, std::array<std::vector<std::int8_t>, 2> &storage)#

Mark the dofs of the test and trial spaces of a that the boundary conditions in bcs constrain.

Markers depend only on the space, so a form whose test and trial spaces are the same marks its dofs once and shares the array between its rows and its columns. The two returned spans are then the same array rather than equal copies.

Parameters:

a – [in] Bilinear form whose spaces are marked.
bcs – [in] Boundary conditions. Only those defined on a space or a subspace of it mark that space.
storage – [out] Backing storage for the returned spans, which are views into it and are valid for as long as it is.

Returns:

Markers on the test space, indexing the rows, and on the trial space, indexing the columns.

template<dolfinx::scalar T, std::floating_point U> std::pair<std::vector<std::vector<std::int8_t>>, std::vector<std::vector<T>>> bc_lifting_data(const std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>> &a, const std::vector<std::vector<std::reference_wrapper<const DirichletBC<T, U>>>> &bcs1)#

Constrained dof markers and boundary condition values on the trial space of each form in a, as used by apply_lifting().

Returns:: Markers and values for each block j, both empty if a[j] is null or bcs1[j] is empty.

template<typename V, std::floating_point U = scalar_value_t<typename std::remove_cvref_t<V>::value_type>, dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type> void apply_lifting(V &&b, const std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>> &a, const std::vector<std::span<const T>> &constants, const std::vector<std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>>> &coeffs, const std::vector<std::span<const std::int8_t>> &bc_markers1, const std::vector<std::span<const T>> &bc_values1, const std::vector<std::span<const T>> &x0, T alpha)#

Modify the right-hand side vector to account for constraints (Dirichlet boundary condition constraints), with the constrained dofs and their values given as arrays.

Computes

\[ b \leftarrow b - \alpha A_{j}^{(1)} (g_{j} - x_{j}) \]

for each block j, as described in the apply_lifting() overload that takes Dirichlet boundary conditions. That overload builds bc_markers1 and bc_values1 from the boundary conditions and calls this function.

Note

Ghost contributions are not accumulated (not sent to owner). Caller is responsible for reverse-scatter to update the ghosts.

Parameters:

b – [inout] The vector to modify inplace.
a – [in] List of bilinear forms, where a[j] is the form that generates the matrix \(A_{j}\). All forms in a must share the same test function space. The trial function spaces can differ.
constants – [in] Constant data appearing in the forms a.
coeffs – [in] Coefficient data appearing in the forms a.
bc_markers1 – [in] Constrained dof markers on the trial space of a[j], owned and ghost (unrolled): bc_markers1[j][i] is non-zero if dof i is constrained. An empty bc_markers1[j] means block j has no constraints and is skipped.
bc_values1 – [in] Boundary condition values \(g_{j}\), with bc_values1[j][i] the value for dof i. Read only where bc_markers1[j][i] is non-zero. Must have the same length as bc_markers1[j].
x0 – [in] The vectors \(x_{j}\). If empty, \(x_{j}\) is treated as zero. Otherwise must have the same length as a.
alpha – [in] Scalar used in the modification of b.

template<typename V, std::floating_point U = scalar_value_t<typename std::remove_cvref_t<V>::value_type>, dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type> void apply_lifting(V &&b, const std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>> &a, const std::vector<std::span<const T>> &constants, const std::vector<std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>>> &coeffs, const std::vector<std::vector<std::reference_wrapper<const DirichletBC<T, U>>>> &bcs1, const std::vector<std::span<const T>> &x0, T alpha)#

Modify the right-hand side vector to account for constraints (Dirichlet boundary condition constraints). This modification is known as ‘lifting’.

Consider the discrete algebraic system

\[\begin{split}\begin{bmatrix} A_{0} & A_{1} \end{bmatrix} \begin{bmatrix} u_{0} \\ u_{1} \end{bmatrix} = b, \end{split}\]

where \(A_{i}\) is a matrix. Partitioning each vector \(u_{i}\) into ‘unknown’ ( \(u_{i}^{(0)}\)) and prescribed ( \(u_{i}^{(1)}\)) groups,

\[\begin{split}\begin{bmatrix} A_{0}^{(0)} & A_{0}^{(1)} & A_{1}^{(0)} & A_{1}^{(1)} \end{bmatrix} \begin{bmatrix} u_{0}^{(0)} \\ u_{0}^{(1)} \\ u_{1}^{(0)} \\ u_{1}^{(1)} \end{bmatrix} = b. \end{split}\]

If \(u_{i}^{(1)} = \alpha(g_{i} - x_{i})\), where \(g_{i}\) is the Dirichlet boundary condition value, \(x_{i}\) is provided and \(\alpha\) is a constant, then

\[\begin{split}\begin{bmatrix} A_{0}^{(0)} & A_{0}^{(1)} & A_{1}^{(0)} & A_{1}^{(1)} \end{bmatrix} \begin{bmatrix} u_{0}^{(0)} \\ \alpha(g_{0} - x_{0}) \\ u_{1}^{(0)} \\ \alpha(g_{1} - x_{1}) \end{bmatrix} = b. \end{split}\]

Rearranging,

\[\begin{split}\begin{bmatrix} A_{0}^{(0)} & A_{1}^{(0)} \end{bmatrix} \begin{bmatrix} u_{0}^{(0)} \\ u_{1}^{(0)} \end{bmatrix} = b - \alpha A_{0}^{(1)} (g_{0} - x_{0}) - \alpha A_{1}^{(1)} (g_{1} - x_{1}). \end{split}\]

The modified \(b\) vector is

\[ b \leftarrow b - \alpha A_{0}^{(1)} (g_{0} - x_{0}) - \alpha A_{1}^{(1)} (g_{1} - x_{1}) \]

More generally,

\[ b \leftarrow b - \alpha A_{i}^{(1)} (g_{i} - x_{i}). \]

Note

Ghost contributions are not accumulated (not sent to owner). Caller is responsible for reverse-scatter to update the ghosts.

Note

Boundary condition values are not set in b by this function. Use DirichletBC::set to set values in b.

Note

Convenience overload for callers that have boundary conditions. It rebuilds the constrained dof markers and values on every call, and should not be called internally by the library; call the overload taking bc_markers1 and bc_values1 instead.

Parameters:

b – [inout] The vector to modify inplace.
a – [in] List of bilinear forms, where a[i] is the form that generates the matrix \(A_{i}\). All forms in a must share the same test function space. The trial function spaces can differ.
constants – [in] Constant data appearing in the forms a.
coeffs – [in] Coefficient data appearing in the forms a.
x0 – [in] The vector \(x_{i}\) above. If empty it is set to zero.
bcs1 – [in] Boundary conditions that provide the \(g_{i}\) values. bcs1[i] is the list of boundary conditions on \(u_{i}\).
alpha – [in] Scalar used in the modification of b.

template<typename V, std::floating_point U = scalar_value_t<typename std::remove_cvref_t<V>::value_type>, dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type> void apply_lifting(V &&b, const std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>> &a, const std::vector<std::vector<std::reference_wrapper<const DirichletBC<T, U>>>> &bcs1, const std::vector<std::span<const T>> &x0, T alpha)#

Modify the right-hand side vector to account for constraints (Dirichlet boundary conditions constraints). This modification is known as ‘lifting’.

See apply_lifting() for a detailed explanation of the lifting. The difference between this function and apply_lifting() is that apply_lifting() requires packed form constant and coefficient data to be passed to the function, whereas this function packs the constant and coefficient form data and then calls apply_lifting().

Note

Ghost contributions are not accumulated (not sent to owner). Caller is responsible for reverse-scatter to update the ghosts.

Note

Boundary condition values are not set in b by this function. Use DirichletBC::set to set values in b.

Note

Convenience overload for callers that have boundary conditions. It rebuilds the constrained dof markers and values on every call, and should not be called internally by the library; call the overload taking bc_markers1 and bc_values1 instead.

Parameters:

b – [inout] The vector to modify inplace.
a – [in] List of bilinear forms, where a[i] is the form that generates the matrix \(A_{i}\) (see apply_lifting()). All forms in a must share the same test function space. The trial function spaces can differ.
x0 – [in] The vector \(x_{i}\) described in apply_lifting(). If empty it is set to zero.
bcs1 – [in] Boundary conditions that provide the \(g_{i}\) values described in apply_lifting(). bcs1[i] is the list of boundary conditions on \(u_{i}\).
alpha – [in] Scalar used in the modification of b (see described in apply_lifting()).

template<dolfinx::scalar T, std::floating_point U> void assemble_matrix(la::MatSet<T> auto mat_add, const Form<T, U> &a, std::span<const T> constants, const std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>> &coefficients, std::span<const std::int8_t> dof_marker0, std::span<const std::int8_t> dof_marker1)#

Assemble bilinear form into a matrix. Matrix must already be initialised. Does not zero or finalise the matrix.

Note

This function can be used to insert kernels into different objects by replacing the mat_add function appropriately.

Template Parameters:

T – scalar type
U – geometry scalar type

Parameters:

mat_add – [in] The function for adding values into the matrix.
a – [in] The bilinear form to assemble.
constants – [in] Constants that appear in a.
coefficients – [in] Coefficients that appear in a.
dof_marker0 – [in] Boundary condition markers for the rows. If bc[i] is true then rows i in A will be zeroed. The index i is a local index.
dof_marker1 – [in] Boundary condition markers for the columns. If bc[i] is true then rows i in A will be zeroed. The index i is a local index.

template<dolfinx::scalar T, std::floating_point U> void assemble_matrix(auto mat_add, const Form<T, U> &a, const std::vector<std::reference_wrapper<const DirichletBC<T, U>>> &bcs)#

Assemble bilinear form into a matrix.

Note

Convenience overload for callers that have boundary conditions. It rebuilds the dof markers on every call, and should not be called internally by the library; call the overload taking dof_marker0 and dof_marker1 instead.

Parameters:

mat_add – [in] The function for adding values into the matrix.
a – [in] The bilinear from to assemble.
bcs – [in] Boundary conditions to apply. For boundary condition dofs the row and column are zeroed. The diagonal entry is not set.

template<dolfinx::scalar T, std::floating_point U> void assemble_matrix(auto mat_add, const Form<T, U> &a, std::span<const std::int8_t> dof_marker0, std::span<const std::int8_t> dof_marker1)#

Assemble bilinear form into a matrix. Matrix must already be initialised. Does not zero or finalise the matrix.

Parameters:

mat_add – [in] The function for adding values into the matrix.
a – [in] The bilinear form to assemble.
dof_marker0 – [in] Boundary condition markers for the rows. If bc[i] is true then rows i in A will be zeroed. The index i is a local index.
dof_marker1 – [in] Boundary condition markers for the columns. If bc[i] is true then rows i in A will be zeroed. The index i is a local index.

template<dolfinx::scalar T, std::floating_point U> void set_diagonal(auto set_fn, const FunctionSpace &V, const std::vector<std::reference_wrapper<const DirichletBC<T, U>>> &bcs, T diagonal = T(1))#

Set a value on the diagonal of the locally owned rows that a Dirichlet boundary condition constrains.

Set only locally owned rows to prevent accumulation when finalising A. A constrained degree-of-freedom that is a ghost on the calling rank is left untouched here and is set by the rank that owns it, so no communication is needed from this function.

This function is typically called after assembly, which zeroes Dirichlet rows and columns. For block matrices, it should normally be called only on the diagonal blocks, i.e. blocks for which the test and trial spaces are the same.

Note

Convenience overload for callers holding V and bcs rather than the row list, which it rebuilds on every call. Library code should cache the rows across repeated calls and set them with la::set_diagonal.

Note

Each row is set exactly once, even where several boundary conditions constrain the same degree-of-freedom, so set_fn may add rather than insert. Every condition sets the same diagonal value, so their order in bcs does not matter here.

Parameters:

set_fn – [in] The function for setting values to a matrix.
V – [in] The function space for the rows and columns of the matrix. It is used to extract only the Dirichlet boundary conditions that are define on V or subspaces of V.
bcs – [in] The Dirichlet boundary conditions. Only conditions defined on V or a subspace of it contribute, and of those only their locally owned dofs.
diagonal – [in] Value to set on the diagonal of each owned constrained row.

Interpolation#

template<std::floating_point T> std::vector<T> dolfinx::fem::interpolation_coords(const fem::FiniteElement<T> &element, const mesh::Geometry<T> &geometry, mesh::CellRange auto &&cells)#

Compute the evaluation points in the physical space at which an expression should be computed to interpolate it in a finite element space.

Parameters:

element – [in] Element to be interpolated into.
geometry – [in] Mesh geometry.
cells – [in] Indices of the cells in the mesh to compute interpolation coordinates for.

Returns:

Coordinates in the physical space at which to evaluate an expression. The shape is (3, num_points) and storage is row-major.

Functions

template<std::floating_point T> mesh::Mesh<T> interpolate_geometry(std::shared_ptr<mesh::Mesh<T>> mesh, const CoordinateElement<T> &new_cmap, const std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)> &reorder_fn = nullptr)#

Take an existing mesh and create a new mesh with its geometry interpolated into a new coordinate element.

The topology is shared between mesh and the returned mesh.

Useful for creating a higher-order mesh from a lower-order one for computation, or vice-versa, for IO.

If new_cmap is discontinuous, the new geometry is discontinuous: each cell has its own coordinate nodes, which are not shared with neighbouring cells.

Note

A discontinuous geometry has no coordinate degrees-of-freedom associated with sub-entities of a cell, so functions that extract the geometry of a sub-entity, e.g. mesh::entities_to_geometry, do not support it.

Parameters:

mesh – [in] Input mesh.
new_cmap – [in] Coordinate element for the new geometry.
reorder_fn – [in] Optional graph reordering function applied to the new geometry dofmap.

Returns:

A new mesh sharing the topology of mesh and with a geometry described by new_cmap.

Sparsity pattern construction#

Functions

template<dolfinx::scalar T, std::floating_point U> std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace>, 2>>> extract_function_spaces(const std::vector<std::vector<const Form<T, U>*>> &a)#

Extract test (0) and trial (1) function spaces pairs for each bilinear form for a rectangular array of forms.

Parameters:: a – [in] A rectangular block on bilinear forms.
Returns:: Rectangular array of the same shape as a with a pair of function spaces in each array entry. If a form is null, then the returned function space pair is (null, null).

template<dolfinx::scalar T, std::floating_point U> la::SparsityPattern create_sparsity_pattern(const Form<T, U> &a)#

Create a sparsity pattern for a given form.

Note

The pattern is not finalised, i.e. the caller is responsible for calling SparsityPattern::assemble.

Parameters:: a – [in] A bilinear form
Returns:: The corresponding sparsity pattern

template<dolfinx::scalar T, std::floating_point U> void build_sparsity_pattern(la::SparsityPattern &pattern, const Form<T, U> &a)#

Build a sparsity pattern for a given form.

Note

The pattern is not finalised, i.e. the caller is responsible for calling SparsityPattern::assemble.

Parameters:

pattern – [in] The sparsity pattern to add to
a – [in] A bilinear form

PETSc helpers#

template<std::floating_point T> Mat dolfinx::fem::petsc::create_matrix(const Form<PetscScalar, T> &a, std::optional<std::string> type = std::nullopt)#

Create a matrix.

Parameters:

a – [in] A bilinear form
type – [in] The PETSc matrix type to create

Throws:

std::invalid_argument – If type is MATIS and the mesh of a has ghost cells. The check is local, but the ghost mode is a mesh property, so every process throws or none does.

Returns:

A sparse matrix with a layout and sparsity that matches the bilinear form. The caller is responsible for destroying the Mat object.

template<std::floating_point T> Mat dolfinx::fem::petsc::create_matrix_block(const std::vector<std::vector<const Form<PetscScalar, T>*>> &a, std::optional<std::string> type = std::nullopt)#

Initialise a monolithic matrix for an array of bilinear forms.

Parameters:

a – [in] Rectangular array of bilinear forms. The a(i, j) form will correspond to the (i, j) block in the returned matrix
type – [in] The type of PETSc Mat. If empty the PETSc default is used.

Throws:

std::invalid_argument – If type is MATIS and the mesh has ghost cells.

Returns:

A sparse matrix with a layout and sparsity that matches the bilinear forms. The caller is responsible for destroying the Mat object.

template<std::floating_point T> Mat dolfinx::fem::petsc::create_matrix_nest(const std::vector<std::vector<const Form<PetscScalar, T>*>> &a, std::optional<std::vector<std::vector<std::optional<std::string>>>> types)#: Create nested (MatNest) matrix.

Note

The caller is responsible for destroying the Mat object.

Vec dolfinx::fem::petsc::create_vector_block(const std::vector<std::pair<std::reference_wrapper<const common::IndexMap>, int>> &maps)#

Initialise monolithic vector. Vector is not zeroed.

The caller is responsible for destroying the Vec object

Vec dolfinx::fem::petsc::create_vector_nest(const std::vector<std::pair<std::reference_wrapper<const common::IndexMap>, int>> &maps)#: Create nested (VecNest) vector. Vector is not zeroed.

template<std::floating_point T> void dolfinx::fem::petsc::assemble_vector(Vec b, const Form<PetscScalar, T> &L, std::span<const PetscScalar> constants, const std::map<std::pair<IntegralType, int>, std::pair<std::span<const PetscScalar>, int>> &coeffs)#

Assemble linear form into an already allocated PETSc vector.

Ghost contributions are not accumulated (not sent to owner). Caller is responsible for calling VecGhostUpdateBegin/End.

Parameters:

b – [inout] The PETsc vector to assemble the form into. The vector must already be initialised with the correct size. The process-local contribution of the form is assembled into this vector. It is not zeroed before assembly.
L – [in] The linear form to assemble
constants – [in] The constants that appear in L
coeffs – [in] The coefficients that appear in L

template<std::floating_point T> void dolfinx::fem::petsc::assemble_vector(Vec b, const Form<PetscScalar, T> &L)#

Assemble linear form into an already allocated PETSc vector.

Ghost contributions are not accumulated (not sent to owner). Caller is responsible for calling VecGhostUpdateBegin/End.

Parameters:

b – [inout] Vector to assemble the form into. The vector must already be initialised with the correct size. The process-local contribution of the form is assembled into this vector. It is not zeroed before assembly.
L – [in] Linear form to assemble.

template<std::floating_point T> void dolfinx::fem::petsc::apply_lifting(Vec b, const std::vector<std::optional<std::reference_wrapper<const Form<PetscScalar, T>>>> &a, const std::vector<std::span<const PetscScalar>> &constants, const std::vector<std::map<std::pair<IntegralType, int>, std::pair<std::span<const PetscScalar>, int>>> &coeffs, const std::vector<std::span<const std::int8_t>> &bc_markers1, const std::vector<std::span<const PetscScalar>> &bc_values1, const std::vector<Vec> &x0, PetscScalar alpha)#

Modify RHS vector to account for Dirichlet boundary conditions, with the constrained dofs and their values given as arrays.

Modify b such that:

b <- b - alpha * A_j (g_j - x0_j)

where j is a block (nest) index. For a non-blocked problem j = 0. The forms in [a] must have the same test space as L (from which b was built), but the trial space may differ. If x0 is not supplied, then it is treated as zero.

Ghost contributions are not accumulated (not sent to owner). Caller is responsible for calling VecGhostUpdateBegin/End.

Parameters:

b – [inout] Vector to modify by lifting.
a – [in] Bilinear forms, one per block j. A std::nullopt entry skips that block.
constants – [in] Constants that appear in each form in a, one entry per block j.
coeffs – [in] Coefficients that appear in each form in a, one entry per block j.
bc_markers1 – [in] Constrained dof markers on the trial space V_j of each block j (owned and ghost, unrolled). An empty entry means block j has no constraints.
bc_values1 – [in] Boundary condition values g_j on V_j, read where bc_markers1[j] is non-zero. Same length as bc_markers1[j].
x0 – [in] Vectors used in the lifting, one per block j. If empty, x0_j is treated as zero for every block. Otherwise must have the same length as a.
alpha – [in] Scaling to apply.

template<std::floating_point T> void dolfinx::fem::petsc::apply_lifting(Vec b, const std::vector<std::optional<std::reference_wrapper<const Form<PetscScalar, T>>>> &a, const std::vector<std::vector<std::reference_wrapper<const DirichletBC<PetscScalar, T>>>> &bcs1, const std::vector<Vec> &x0, PetscScalar alpha)#

Modify RHS vector to account for Dirichlet boundary conditions.

Modify b such that:

b <- b - alpha * A_j (g_j - x0_j)

where j is a block (nest) index. For a non-blocked problem j = 0. The boundary conditions bcs1 are on the trial spaces V_j. The forms in [a] must have the same test space as L (from which b was built), but the trial space may differ. If x0 is not supplied, then it is treated as zero.

Ghost contributions are not accumulated (not sent to owner). Caller is responsible for calling VecGhostUpdateBegin/End.

Note

Convenience overload for callers that have boundary conditions. It rebuilds the constrained dof markers and values on every call, and should not be called internally by the library; call the overload taking bc_markers1 and bc_values1 instead.

Parameters:

b – [inout] Vector to modify by lifting.
a – [in] Bilinear forms, one per block j. A std::nullopt entry skips that block.
bcs1 – [in] Boundary conditions on the trial space V_j for each block j.
x0 – [in] Vectors used in the lifting, one per block j. If empty, x0_j is treated as zero for every block. Otherwise must have the same length as a.
alpha – [in] Scaling to apply.

template<std::floating_point T> void dolfinx::fem::petsc::set_bc(Vec b, const std::vector<std::reference_wrapper<const DirichletBC<PetscScalar, T>>> &bcs, std::optional<const Vec> x0, PetscScalar alpha = 1)#

Entries in b that are constrained by a Dirichlet boundary conditions are set to alpha * (x_bc - x0), where x_bc is the (interpolated) boundary condition value.

Parameters:

b – [in] The vector to apply the boundary condition to. The local (owned) part of this vector is modified. The user is responsible for scattering the changes to the ghost part of the vector if necessary.
bcs – [in] The boundary conditions to apply.
x0 – [in] Optional vector used in computing the value to set. If not provided it is treated as zero. The local (owned) part of this vector is used.
alpha – [in] Scaling to apply.

template<std::floating_point T> void dolfinx::fem::petsc::assemble_residual(const Vec x, Vec b, const Form<PetscScalar, T> &F, const Form<PetscScalar, T> &J, const std::vector<std::reference_wrapper<const DirichletBC<PetscScalar, T>>> &bcs, Function<PetscScalar, T> &u)#

Assemble the residual \(F(x)\) of a nonlinear problem into b, with Dirichlet conditions applied.

Intended as the body of the residual callback of nls::petsc::SNESSolver, which passes the point to evaluate at x and the vector to assemble into b:

solver.set_F([&](const Vec x, Vec b)
             { assemble_residual(x, b, F, J, bcs, u); }, b_layout);

Entries of b constrained by bcs are set to x - g, so that a Newton update drives x to the boundary condition value g.

Parameters:

x – [in] Point at which to evaluate the residual, e.g. a line search trial point. Must be ghosted. Its ghost values are updated before use.
b – [out] Vector to assemble into, which is the one the solver passed to the callback and not necessarily the one registered with set_F. Zeroed first, and its ghost values are updated on return.
F – [in] Residual form.
J – [in] Jacobian form, used to lift bcs.
bcs – [in] Dirichlet boundary conditions.
u – [out] Function that F and J hold as a coefficient. Its degrees-of-freedom are set to x before assembly.

template<std::floating_point T> void dolfinx::fem::petsc::assemble_jacobian(const Vec x, Mat Jmat, Mat Pmat, const Form<PetscScalar, T> &J, const std::vector<std::reference_wrapper<const DirichletBC<PetscScalar, T>>> &bcs, Function<PetscScalar, T> &u, const Form<PetscScalar, T> *P = nullptr)#

Assemble the Jacobian \(dF/dx\) of a nonlinear problem into Jmat, and a preconditioner into Pmat.

Intended as the body of the Jacobian callback of nls::petsc::SNESSolver, which passes the point to evaluate at x and the matrices to assemble into Jmat and Pmat:

solver.set_J([&](const Vec x, Mat Jmat, Mat Pmat)
             { assemble_jacobian(x, Jmat, Pmat, J, bcs, u); },
             A_layout);

Rows and columns constrained by bcs are zeroed, and for a form whose test and trial spaces are the same a unit diagonal is set on the constrained rows, matching the residual assembled by assemble_residual.

Parameters:

x – [in] Point at which to evaluate the Jacobian, e.g. a line search trial point. Must be ghosted. Its ghost values are updated before use.
Jmat – [out] Matrix to assemble the Jacobian into, which is the one the solver passed to the callback and not necessarily the one registered with set_J. Zeroed first.
Pmat – [out] Matrix to assemble the preconditioner into. Zeroed first. Unused, and may be nullptr, if P is not given.
J – [in] Jacobian form.
bcs – [in] Dirichlet boundary conditions.
u – [out] Function that J and P hold as a coefficient. Its degrees-of-freedom are set to x before assembly.
P – [in] Preconditioner form. If not given, Pmat is left alone and PETSc preconditions with the Jacobian.

Pre:

P, if given, must have the same function spaces as J, so that the constrained dof markers built for J apply to it too.

Creating forms, function spaces and expressions#

Functions

std::vector<std::string> get_coefficient_names(const ufcx_form &ufcx_form)#

Get the name of each coefficient in a UFC form

Parameters:: ufcx_form – [in] The UFC form
Returns:: The name of each coefficient

std::vector<std::string> get_constant_names(const ufcx_form &ufcx_form)#

Get the name of each constant in a UFC form.

Parameters:: ufcx_form – [in] The UFC form
Returns:: The name of each constant

template<dolfinx::scalar T, std::floating_point U = scalar_value_t<T>> Form<T, U> create_form_factory(const std::vector<std::reference_wrapper<const ufcx_form>> &ufcx_forms, const std::vector<std::shared_ptr<const FunctionSpace>> &spaces, const std::vector<std::shared_ptr<const Function<T, U>>> &coefficients, const std::vector<std::shared_ptr<const Constant<T>>> &constants, const std::map<IntegralType, std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> &subdomains, const std::vector<std::reference_wrapper<const mesh::EntityMap>> &entity_maps, std::shared_ptr<const mesh::Mesh> mesh = nullptr)#

Create a Form from UFCx input with coefficients and constants passed in the required order.

Use fem::create_form to create a fem::Form with coefficients and constants associated with the name/string.

Note

Argument and coefficient elements are checked against the form only where UFCx supplies a nonzero Basix element hash. The hash of a blocked element is that of its scalar sub-element, so a block size (value shape) mismatch escapes the check, as does any mismatch in an element that has no Basix hash, e.g. a mixed element. Passing a space or coefficient that differs from the compiled form in either respect is undefined behaviour: the generated kernel reads a coefficient array packed to the wrong size.

Parameters:

ufcx_forms – [in] A list of UFCx forms, one for each cell type.
spaces – [in] Vector of function spaces. The number of spaces is equal to the rank of the form.
coefficients – [in] Coefficient fields in the form.
constants – [in] Spatial constants in the form.
subdomains – [in] Subdomain markers. The data can be computed using fem::compute_integration_domains.
entity_maps – [in] The entity maps for the form. Empty for single domain problems.
mesh – [in] The mesh of the domain.

Pre:

Each value in subdomains must be sorted by domain id.

template<dolfinx::scalar T, std::floating_point U = scalar_value_t<T>> Form<T, U> create_form(const ufcx_form &ufcx_form, const std::vector<std::shared_ptr<const FunctionSpace>> &spaces, const std::map<std::string, std::shared_ptr<const Function<T, U>>> &coefficients, const std::map<std::string, std::shared_ptr<const Constant<T>>> &constants, const std::map<IntegralType, std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> &subdomains, const std::vector<std::reference_wrapper<const mesh::EntityMap>> &entity_maps, std::shared_ptr<const mesh::Mesh> mesh = nullptr)#

Create a Form from UFC input with coefficients and constants resolved by name.

Parameters:

ufcx_form – [in] UFC form
spaces – [in] Function spaces for the Form arguments.
coefficients – [in] Coefficient fields in the form (by name).
constants – [in] Spatial constants in the form (by name).
subdomains – [in] Subdomain markers. The data can be computed using fem::compute_integration_domains.
entity_maps – [in] The entity maps for the form. Empty for single domain problems.
mesh – [in] Mesh of the domain. This is required if the form has no arguments, e.g. a functional.

Pre:

Each value in subdomains must be sorted by domain id.

Returns:

A Form

template<dolfinx::scalar T, std::floating_point U = scalar_value_t<T>> Form<T, U> create_form(ufcx_form *(*fptr)(), const std::vector<std::shared_ptr<const FunctionSpace>> &spaces, const std::map<std::string, std::shared_ptr<const Function<T, U>>> &coefficients, const std::map<std::string, std::shared_ptr<const Constant<T>>> &constants, const std::map<IntegralType, std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>> &subdomains, const std::vector<std::reference_wrapper<const mesh::EntityMap>> &entity_maps, std::shared_ptr<const mesh::Mesh> mesh = nullptr)#

Create a Form using a factory function that returns a pointer to a ufcx_form.

Coefficients and constants are resolved by name/string.

Parameters:

fptr – [in] Pointer to a function returning a pointer to ufcx_form.
spaces – [in] Function spaces for the Form arguments.
coefficients – [in] Coefficient fields in the form (by name),
constants – [in] Spatial constants in the form (by name),
subdomains – [in] Subdomain markers. The data can be computed using fem::compute_integration_domains.
entity_maps – [in] The entity maps for the form. Empty for single domain problems.
mesh – [in] Mesh of the domain. This is required if the form has no arguments, e.g. a functional.

Pre:

Each value in subdomains must be sorted by domain id.

Returns:

A Form

Functions

template<std::floating_point T> ElementDofLayout create_element_dof_layout(const fem::FiniteElement<T> &element, const std::vector<int> &parent_map = {})#: Create an ElementDofLayout from a FiniteElement.

DofMap create_dofmap(MPI_Comm comm, const ElementDofLayout &layout, mesh::Topology &topology, const std::function<void(std::span<std::int32_t>, std::uint32_t)> &permute_inv, const std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)> &reorder_fn)#

Create a dof map on mesh.

Parameters:

comm – [in] MPI communicator
layout – [in] Dof layout on an element
topology – [in] Mesh topology
permute_inv – [in] Function to un-permute dofs. nullptr when transformation is not required.
reorder_fn – [in] Graph reordering function called on the dofmap

Returns:

A new dof map

std::vector<DofMap> create_dofmaps(MPI_Comm comm, const std::vector<ElementDofLayout> &layouts, mesh::Topology &topology, const std::function<void(std::span<std::int32_t>, std::uint32_t)> &permute_inv, const std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)> &reorder_fn)#

Create a set of dofmaps on a given topology.

Note

The number of layouts must match the number of cell types in the topology

Parameters:

comm – [in] MPI communicator
layouts – [in] Dof layout on each element type
topology – [in] Mesh topology
permute_inv – [in] Function to un-permute dofs. nullptr when transformation is not required.
reorder_fn – [in] Graph reordering function called on the dofmaps

Returns:

The list of new dof maps

template<std::floating_point T> FunctionSpace<T> create_functionspace(std::shared_ptr<mesh::Mesh<T>> mesh, std::shared_ptr<const fem::FiniteElement<T>> e, std::function<std::vector<int>(const graph::AdjacencyList<std::int32_t>&)> reorder_fn = nullptr)#: NEW Create a function space from a fem::FiniteElement.

Functions

template<dolfinx::scalar T, std::floating_point U = scalar_value_t<T>> Expression<T, U> create_expression(const ufcx_expression &e, const std::vector<std::shared_ptr<const Function<T, U>>> &coefficients, const std::vector<std::shared_ptr<const Constant<T>>> &constants, const std::vector<std::reference_wrapper<const mesh::EntityMap>> &entity_maps, std::shared_ptr<const FunctionSpace> argument_space = nullptr)#: Create Expression from UFC.

template<dolfinx::scalar T, std::floating_point U = scalar_value_t<T>> Expression<T, U> create_expression(const ufcx_expression &e, const std::map<std::string, std::shared_ptr<const Function<T, U>>> &coefficients, const std::map<std::string, std::shared_ptr<const Constant<T>>> &constants, const std::vector<std::reference_wrapper<const mesh::EntityMap>> &entity_maps, std::shared_ptr<const FunctionSpace> argument_space = nullptr)#: Create Expression from UFC input (with named coefficients and constants).

Integration domains#

Functions

template<int num_cells> std::array<std::int32_t, 2 * num_cells> get_cell_facet_pairs(std::int32_t f, std::span<const std::int32_t> cells, const graph::AdjacencyList<std::int32_t> &c_to_f)#

Helper function to get an array of (cell, local_facet) pairs corresponding to a given facet index.

Parameters:

f – [in] Facet index
cells – [in] List of cells incident to the facet
c_to_f – [in] Cell to facet connectivity

Returns:

Vector of (cell, local_facet) pairs

template<int num_cells> std::array<std::int32_t, 2 * num_cells> get_cell_entity_pairs(std::int32_t e, std::span<const std::int32_t> cells, const graph::AdjacencyList<std::int32_t> &c_to_e)#

Helper function to get an array of of (cell, local_entity) pairs corresponding to a given entity index.

Note

If the entity is connected to multiple cells, the first one is picked.

Parameters:

e – [in] entity index
cells – [in] List of cells incident to the entity
c_to_e – [in] Cell to entity connectivity

Returns:

Vector of (cell, local_entity) pairs

std::vector<std::int32_t> compute_integration_domains(IntegralType integral_type, const mesh::Topology &topology, std::span<const std::int32_t> entities)#

Given an integral type and a set of entities, computes and return data for the entities that should be integrated over.

This function returns a list data, for each entity in entities, that is used in assembly. For cell integrals it is simply the cell cell indices. For exterior facet integrals, a list of (cell_index, local_facet_index) pairs is returned. For interior facet integrals, a list of (cell_index0, local_facet_index0, cell_index1, local_facet_index1) tuples is returned. The data computed by this function is typically used as input to fem::create_form.

Note

Owned mesh entities only are returned. Ghost entities are not included.

Parameters:

integral_type – [in] Integral type.
topology – [in] Mesh topology.
entities – [in] List of mesh entities. Depending on the IntegralType these are associated with different entities: IntegralType::cell: cells IntegralType::exterior_facet: facets IntegralType::interior_facet: facets IntegralType::vertex: vertices

Pre:

For facet integrals, the topology facet-to-cell and cell-to-facet connectivity must be computed before calling this function.

Returns:

List of integration entity data, depending on the IntegralType the data per entity has different layouts IntegralType::cell: cell IntegralType::exterior_facet: (cell, local_facet) IntegralType::interior_facet: (cell, local_facet) IntegralType::vertex: (cell, local_vertex)

Finite element (dolfinx::fem)#

Finite elements#

Function spaces and functions#

Function spaces#

Functions#

Constants#

Forms#

Dirichlet boundary conditions#

Degree-of-freedom maps#

Assembly#

Interpolation#

Sparsity pattern construction#

PETSc helpers#

Creating forms, function spaces and expressions#

Integration domains#

Finite element (`dolfinx::fem`)#