Finite element (dolfinx::fem)

Under development

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.

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::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.

  • 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 ony set for 2nd-order tensor blocked elements.

FiniteElement(std::vector<BasixElementData<geometry_type>> elements)

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.

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.

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.

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.

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. For blocked elements the block size is equal to the value size, except for symmetric rank-2 tensor blocked elements. For a symmetric rank-2 tensor blocked element the block size is 3 in 2D and 6 in 3D.

Returns:

Block size of the finite element space.

int value_size() const

Value size of the finite element field.

The value size is the number of components in the finite element field. It is the product of the value shape, e.g. is is 1 for a scalar function, 2 for a 2D vector, 9 for a second-order tensor in 3D, etc. For blocked elements, this function returns the value size for the full ‘blocked’ element.

Note

The return value of this function is inconsistent with value_shape() for rank-2 ‘symmetric’ elements. Due to issues elsewhere in the code base, rank-2 symmetric fields have value shape {3} (2D) or {6} rather than {2, 2} and {3, 3}, respectively. For symmetric rank-2 tensors this function returns 4 for 2D cases and 9 for 3D cases. This inconsistency will be fixed in the future.

Throws:

Exception – is thrown for a mixed element as mixed elements do 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.

The value shape describes the shape of the finite element field, e.g. {} for a scalar, {2} for a vector in 2D, {3, 3} for a rank-2 tensor in 3D, etc.

Throws:

Exception – is thrown for a mixed element as mixed elements do not have a value shape.

Returns:

The value shape.

int reference_value_size() const

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

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.

For blocked elements, this function returns the value shape for the ‘base’ element from which the blocked element is composed. For other elements, the return value is the same as FiniteElement::value_shape.

Throws:

Exception – is thrown for a mixed element as mixed elements do not have a 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.

For non-blocked elements, this function returns the same as FiniteElement::value_shape. For blocked and quadrature elements the returned shape will be {}.

Mixed elements do not have a reference value shape.

Throws:

Exception – is thrown for a mixed element as mixed elements do not have a 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?

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:

Throws – a std::runtime_error is there 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 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<U>, 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<U>, 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<U> 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<U> 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<U> 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<U> 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<U> 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<U> 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<U> 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<U> 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] Indicies 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] Indicies 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.

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.

FunctionSpace(FunctionSpace &&V) = default

Move constructor.

virtual ~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::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 DofMap> dofmap() const

The dofmap.

Functions

template<dolfinx::scalar T, std::floating_point U>
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)

Create function on given function space.

Parameters:

V – [in] The function space

inline Function(std::shared_ptr<const FunctionSpace<geometry_type>> V, std::shared_ptr<la::Vector<value_type>> x)

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.

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>>(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<const geometry_type, MDSPAN_IMPL_STANDARD_NAMESPACE::extents<std::size_t, 3, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>)> &f)

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

Parameters:

f – [in] Expression to be interpolated.

inline void interpolate(const std::function<std::pair<std::vector<value_type>, std::vector<std::size_t>>(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<const geometry_type, MDSPAN_IMPL_STANDARD_NAMESPACE::extents<std::size_t, 3, MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>)> &f, std::span<const std::int32_t> 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 Function<value_type, geometry_type> &u)

Interpolate a Function over all cells.

Parameters:

u – [in] Function to be interpolated.

Pre:

The mesh associated with this and the mesh associated with u must be the same mesh::Mesh.

inline void interpolate(const Function<value_type, geometry_type> &u0, std::span<const std::int32_t> cells0, std::span<const std::int32_t> 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. This argument can be empty when this and u0 share the same mesh. Otherwise the length of cells and the length of cells0 must be the same.

inline void interpolate(const Expression<value_type, geometry_type> &e)

Interpolate an Expression on all cells.

Parameters:

e – [in] Expression to be interpolated.

Pre:

If a mesh is associated with Function coefficients of e, it must be the same as the mesh::Mesh associated with this.

inline void interpolate(const Expression<value_type, geometry_type> &e0, std::span<const std::int32_t> cells0, std::span<const std::int32_t> cells1 = {})

Interpolate an Expression over a subset of cells.

Parameters:

  • e0 – [in] Expression to be interpolated. The Expression must have been created using the reference coordinates created by FiniteElement::interpolation_points for the element associated with this.

  • cells0 – [in] Cells in the mesh associated with e0 to interpolate from if e0 has Function coefficients. If no mesh can be associated with e0 then the mesh associated with this is used.

  • 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. This argument can be empty when this and u0 share the same mesh. Otherwise the length of cells and the length of cells0 must be the same.

inline void interpolate(const Function<value_type, geometry_type> &v, std::span<const std::int32_t> cells, const geometry::PointOwnershipData<U> &interpolation_data)

Interpolate a Function defined on a different mesh.

Parameters:

  • v – [in] Function to be interpolated.

  • cells – [in] Cells in the mesh associated with this to interpolate into.

  • interpolation_data – [in] Data required for associating the interpolation points of this with cells in v. Can be computed with fem::create_interpolation_data.

inline void eval(std::span<const geometry_type> x, std::array<std::size_t, 2> xshape, std::span<const std::int32_t> cells, std::span<value_type> u, std::array<std::size_t, 2> ushape) 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.

Public Members

std::string name = "u"

Name.

Constants

template<dolfinx::scalar T>
class Constant

Constant 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_type_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.

  • Kern – Element kernel.

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, 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::map<std::shared_ptr<const mesh::Mesh<geometry_type>>, std::span<const std::int32_t>> &entity_maps, std::shared_ptr<const mesh::Mesh<geometry_type>> mesh = nullptr)

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.

  • integrals – [in] The integrals in the form. For each integral type, there is a list of integral data.

  • coefficients – [in] Coefficients in the form.

  • constants – [in] Constants in the form.

  • needs_facet_permutations – [in] Set to true is any of the integration kernels require cell permutation data.

  • entity_maps – [in] If any trial functions, test functions, or coefficients in the form are not defined over the same mesh as the integration domain, entity_maps must be supplied. For each key (a mesh, different to the integration domain mesh) a map should be provided relating the entities in the integration domain mesh to the entities in the key mesh e.g. for a pair (msh, emap) in entity_maps, emap[i] is the entity in msh corresponding to entity i in the integration domain mesh.

  • mesh – [in] Mesh of the domain. This is required when there are no argument functions from which the mesh can be extracted, e.g. for functionals.

Pre:

The integral data in integrals must be sorted by domain (domain id).

Form(const Form &form) = delete

Copy constructor.

Form(Form &&form) = default

Move constructor.

virtual ~Form() = default

Destructor.

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

Extract common mesh for the form.

Returns:

The 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 uint8_t*)> kernel(IntegralType type, int i) const

Get the kernel function for integral i on given domain type.

Parameters:

  • type – [in] Integral type.

  • i – [in] Domain identifier (index).

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 int num_integrals(IntegralType type) const

Number of integrals on given domain type.

Parameters:

type – [in] Integral type.

Returns:

Number of integrals.

inline std::vector<int> active_coeffs(IntegralType type, std::size_t i) const

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

A form is split into multiple integrals (kernels) and each integral might container 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.

  • i – [in] Index of the integral.

inline std::vector<int> integral_ids(IntegralType type) const

Get the IDs for integrals (kernels) for given integral type.

The IDs correspond to the domain IDs which the integrals are defined for in the form. ID=-1 is the default integral over the whole domain.

Parameters:

type – [in] Integral type.

Returns:

List of IDs for given integral type.

inline std::span<const std::int32_t> domain(IntegralType type, int i) const

Get the list of mesh entity indices for the ith integral (kernel) of a given type.

For IntegralType::cell, returns a list of cell indices.

For IntegralType::exterior_facet, returns a list of (cell_index, local_facet_index) pairs. Data is flattened with row-major layout, shape=(num_facets, 2).

For IntegralType::interior_facet, returns list of tuples of the form (cell_index_0, local_facet_index_0, cell_index_1, / local_facet_index_1). Data is flattened with row-major layout, shape=(num_facets, 4).

Parameters:

  • type – [in] Integral domain type.

  • i – [in] Integral ID, i.e. (sub)domain index.

Returns:

List of active entities for the given integral (kernel).

inline std::vector<std::int32_t> domain(IntegralType type, int i, const mesh::Mesh<geometry_type> &mesh) const

Compute the list of entity indices in mesh for the ith integral (kernel) of a given type (i.e. cell, exterior facet, or interior facet).

Parameters:

  • type – Integral type.

  • i – Integral ID, i.e. the (sub)domain index.

  • mesh – The mesh the entities are numbered with respect to.

Returns:

List of active entities in mesh for the given integral.

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

Access coefficients.

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.

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

Access constants.

Dirichlet boundary conditions

template<dolfinx::scalar T, std::floating_point U = dolfinx::scalar_value_type_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<U>> 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 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 sorted.

  • V – [in] The function space to be constrained

Pre:

dofs must be sorted.

template<typename X>
inline DirichletBC(std::shared_ptr<const Constant<T>> g, X &&dofs, std::shared_ptr<const FunctionSpace<U>> 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 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 sorted.

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 sorted.

template<typename X>
inline DirichletBC(std::shared_ptr<const Function<T, U>> g, X &&V_g_dofs, std::shared_ptr<const FunctionSpace<U>> 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.

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<U>> 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

Access dof indices (local indices, unrolled), including ghosts, to which a Dirichlet condition is applied, and the index to the first non-owned (ghost) index. The array of indices is sorted.

Returns:

Sorted array of dof indices (unrolled) and index to the first entry in the dof index array that is not owned. 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.

If x includes ghosted entries (entries available on the calling rank but owned by another rank), ghosted entries constrained by a Dirichlet condition will also be set.

Parameters:

  • x – [inout] Array to modify for Dirichlet boundary conditions.

  • 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(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<const std::int32_t, MDSPAN_IMPL_STANDARD_NAMESPACE::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 3 Global 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 5 Unassembled 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

MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<const std::int32_t, MDSPAN_IMPL_STANDARD_NAMESPACE::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.

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<dolfinx::scalar T, std::floating_point U>
void assemble_vector(std::span<T> 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<dolfinx::scalar T, std::floating_point U>
void assemble_vector(std::span<T> b, const Form<T, U> &L)

Assemble linear form into a vector.

Parameters:

  • b – [inout] The vector to be assembled. It will not be zeroed before assembly.

  • L – [in] The linear forms to assemble into b

template<dolfinx::scalar T, std::floating_point U>
void apply_lifting(std::span<T> b, 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 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.

template<dolfinx::scalar T, std::floating_point U>
void apply_lifting(std::span<T> b, 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 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.

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.

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, std::span<const T> constants, const std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>> &coefficients, const std::vector<std::reference_wrapper<const DirichletBC<T, U>>> &bcs)

Assemble bilinear form into a matrix

Parameters:

  • mat_add – [in] The function for adding values into the matrix

  • a – [in] The bilinear from to assemble

  • constants – [in] Constants that appear in a

  • coefficients – [in] Coefficients that appear in a

  • 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, const std::vector<std::reference_wrapper<const DirichletBC<T, U>>> &bcs)

Assemble bilinear form into a matrix

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>
void set_diagonal(auto set_fn, std::span<const std::int32_t> rows, T diagonal = 1.0)

Sets a value to the diagonal of a matrix for specified rows.

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

Parameters:

  • set_fn – [in] The function for setting values to a matrix

  • rows – [in] The row blocks, in local indices, for which to add a value to the diagonal

  • diagonal – [in] The value to add to the diagonal for the specified rows

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

Sets a value to the diagonal of the matrix for rows with a Dirichlet boundary conditions applied.

This function is typically called after assembly. The assembly function zeroes Dirichlet rows and columns. This function adds the value only to rows that are locally owned, and therefore does not create a need for parallel communication. For block matrices, this function should normally be called only on the diagonal blocks, i.e. blocks for which the test and trial spaces are the same.

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

  • diagonal – [in] The value to add to the diagonal for rows with a boundary condition applied

Interpolation

template<std::floating_point T>
std::vector<T> dolfinx::fem::interpolation_coords(const fem::FiniteElement<T> &element, const mesh::Geometry<T> &geometry, std::span<const std::int32_t> 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:

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

Sparsity pattern construction

template<dolfinx::scalar T, std::floating_point U>
la::SparsityPattern dolfinx::fem::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

PETSc helpers

Warning

doxygennamespace: Cannot find namespace “dolfinx::fem::petsc” in doxygen xml output for project “DOLFINx” from directory: ../xml/

Misc

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

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. For integral_type==IntegralType::cell, entities should be cell indices. For other IntegralType, entities should be facet indices.

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.

template<dolfinx::scalar T, std::floating_point U>
std::vector<std::vector<std::array<std::shared_ptr<const FunctionSpace<U>>, 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<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, std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv, 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, std::function<void(std::span<std::int32_t>, std::uint32_t)> permute_inv, 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

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_type_t<T>>
Form<T, U> create_form_factory(const ufcx_form &ufcx_form, const std::vector<std::shared_ptr<const FunctionSpace<U>>> &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::map<std::shared_ptr<const mesh::Mesh<U>>, std::span<const std::int32_t>> &entity_maps, std::shared_ptr<const mesh::Mesh<U>> 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.

Parameters:

  • ufcx_form – [in] The UFCx form.

  • 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_type_t<T>>
Form<T, U> create_form(const ufcx_form &ufcx_form, const std::vector<std::shared_ptr<const FunctionSpace<U>>> &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::map<std::shared_ptr<const mesh::Mesh<U>>, std::span<const std::int32_t>> &entity_maps, std::shared_ptr<const mesh::Mesh<U>> 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_type_t<T>>
Form<T, U> create_form(ufcx_form *(*fptr)(), const std::vector<std::shared_ptr<const FunctionSpace<U>>> &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::map<std::shared_ptr<const mesh::Mesh<U>>, std::span<const std::int32_t>> &entity_maps, std::shared_ptr<const mesh::Mesh<U>> 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

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.

template<dolfinx::scalar T, std::floating_point U>
std::span<const std::uint32_t> get_cell_orientation_info(const Function<T, U> &coefficient)
template<int _bs, dolfinx::scalar T>
void pack(std::span<T> coeffs, std::int32_t cell, int bs, std::span<const T> v, std::span<const std::uint32_t> cell_info, const DofMap &dofmap, auto transform)

Pack a single coefficient for a single cell.

template<dolfinx::scalar T, std::floating_point U>
void pack_coefficient_entity(std::span<T> c, int cstride, const Function<T, U> &u, std::span<const std::uint32_t> cell_info, std::span<const std::int32_t> entities, std::size_t estride, FetchCells auto &&fetch_cells, std::int32_t offset)

Pack a single coefficient for a set of active entities.

Parameters:

  • c – [out] Coefficient to be packed.

  • cstride – [in] Total number of coefficient values to pack for each entity.

  • u – [in] Function to extract coefficient data from.

  • cell_info – [in] Array of bytes describing which transformation has to be applied on the cell to map it to the reference element.

  • entities – [in] Set of active entities.

  • estride – [in] Stride for each entity in active entities.

  • fetch_cells – [in] Function that fetches the cell index for an entity in active_entities.

  • offset – [in] The offset for c.

template<dolfinx::scalar T, std::floating_point U>
std::pair<std::vector<T>, int> allocate_coefficient_storage(const Form<T, U> &form, IntegralType integral_type, int id)

Allocate storage for coefficients of a pair (integral_type, / id) from a Form.

Parameters:

  • form – [in] The Form

  • integral_type – [in] Type of integral

  • id – [in] The id of the integration domain

Returns:

A storage container and the column stride

template<dolfinx::scalar T, std::floating_point U>
std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>> allocate_coefficient_storage(const Form<T, U> &form)

Allocate memory for packed coefficients of a Form.

Parameters:

form – [in] The Form

Returns:

Map from a form (integral_type, domain_id) pair to a (coeffs, cstride) pair

template<dolfinx::scalar T, std::floating_point U>
void pack_coefficients(const Form<T, U> &form, IntegralType integral_type, int id, std::span<T> c, int cstride)

Pack coefficients of a Form for a given integral type and domain id.

Parameters:

  • form – [in] The Form

  • integral_type – [in] Type of integral

  • id – [in] The id of the integration domain

  • c – [inout] The coefficient array

  • cstride – [in] The coefficient stride

template<dolfinx::scalar T, std::floating_point U = scalar_value_type_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, std::shared_ptr<const FunctionSpace<U>> argument_function_space = nullptr)

Create Expression from UFC.

template<dolfinx::scalar T, std::floating_point U = scalar_value_type_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, std::shared_ptr<const FunctionSpace<U>> argument_function_space = nullptr)

Create Expression from UFC input (with named coefficients and constants).

template<dolfinx::scalar T, std::floating_point U>
void pack_coefficients(const Form<T, U> &form, std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>> &coeffs)

Pack coefficients of a Form.

Warning

This is subject to change

Parameters:

  • form – [in] The Form

  • coeffs – [inout] A map from an (integral_type, domain_id) pair to a (coeffs, cstride) pair. coeffs is a storage container representing an array of shape (num_int_entities, cstride) in which to pack the coefficient data, where num_int_entities is the number of entities being integrated over and cstride is the number of coefficient data entries per integration entity. coeffs is flattened into row-major layout.

template<dolfinx::scalar T, std::floating_point U>
std::pair<std::vector<T>, int> pack_coefficients(const Expression<T, U> &e, std::span<const std::int32_t> entities, std::size_t estride)

Pack coefficients of a Expression u for a give list of active entities.

Parameters:

  • e – [in] The Expression

  • entities – [in] A list of active entities

  • estride – [in] Stride for each entity in active entities (1 for cells, 2 for facets)

Returns:

A pair of the form (coeffs, cstride)

template<typename U>
std::vector<typename U::scalar_type> pack_constants(const U &u)

Pack constants of u into a single array ready for assembly.

Warning

This function is subject to change.