element

This submodule defines operations and utilities for dealing with elements. This means both organizing them, as well as using that data for iterating over them for computations.

At the core of organization is the ArrayCom type, which holds information necessary to coordinate operations over all elements. For now that is just the number of all elements, though in the future, this is probably where OpenMPI data would be stored.

class mfv2d.element.ArrayCom(element_cnt: int)

Shared array coordination type.

Parameters:

element_cnt (int) – Number of elements in the array.

element_cnt: int

Storing Data

From ArrayCom object element arrays can be created. These contain data for each of the elements. The base for these objects is the ElementArrayBase type.

class mfv2d.element.ElementArrayBase(com: ArrayCom, values: tuple[npt.NDArray[_ElementDataType], ...], dtype: type[_ElementDataType])

Base of element arrays.

Implements common methods for all element arrays, such as assignment, getting values, etc.

Parameters:

• com (ArrayCom) – Array coordination object.

• values (tuple of arrays) – Tuple that holds arrays for each element.

• dtype (type) – Type of the elements which are stored in the arrays.

__getitem__(i: SupportsIndex, /) → npt.NDArray[_ElementDataType]

Return element’s values.

__iter__() → Iterator[npt.NDArray[_ElementDataType]]

Iterate over all elements.

__len__() → int

Return the number of elements.

__setitem__(i: SupportsIndex, val: array_like, /) → None

Set the value for the current element.

com: ArrayCom

dtype: type[_ElementDataType]

unique(axis: int | None = None) → npt.NDArray[_ElementDataType]

Find unique values within the array.

values: tuple[ndarray[tuple[Any, ...], dtype[_ElementDataType]], ...]

The first subtype of ElementArrayBase is the FixedElementArray. This stores data as arrays of equal shape for each element.

class mfv2d.element.FixedElementArray(com: ArrayCom, shape: int | Sequence[int], dtype: type[_ElementDataType], /)

Array with values and fixed shape per element spread over elements.

This is used for storing data which must have the same shape for each element, such as the corners of the element.

Parameters:

• com (ArrayCom) – Array coordination object.

• values (tuple of arrays) – Tuple that holds arrays for each element.

• shape (int or Sequence of ints) – Shape of the array used for all elements.

• dtype (type) – Type of the elements which are stored in the arrays.

__getitem__(i: SupportsIndex, /) → npt.NDArray[_ElementDataType]

Return element’s values.

__iter__() → Iterator[npt.NDArray[_ElementDataType]]

Iterate over all elements.

__len__() → int

Return the number of elements.

__setitem__(i: SupportsIndex, val: array_like, /) → None

Set the value for the current element.

com: ArrayCom

dtype: type[_ElementDataType]

shape: tuple[int, ...]

unique(axis: int | None = None) → npt.NDArray[_ElementDataType]

Find unique values within the array.

values: tuple[ndarray[tuple[Any, ...], dtype[_ElementDataType]], ...]

The second subtype is FlexibleElementArray. This type stores data as arrays with the same number of dimensions, but different shapes for each element. As such these arrays also contain a FixedElementArray, which acts as storage for the shapes of individual element arrays.

class mfv2d.element.FlexibleElementArray(com: ArrayCom, dtype: type[_ElementDataType], shapes: FixedElementArray[_IntegerDataType])

Array with values and variable shape per element.

This is used for storing data which can have different shapes, but same number of dimensions for each element, such as the element matrices or degrees of freedom.

Parameters:

• com (ArrayCom) – Array coordination object.

• dtype (type) – Type of the elements which are stored in the arrays.

• shapes (FixedElementArray) – Element array containing shapes of the elements.

• values (tuple of arrays) – Tuple that holds arrays for each element.

__getitem__(i: SupportsIndex, /) → npt.NDArray[_ElementDataType]

Return element’s values.

__iter__() → Iterator[npt.NDArray[_ElementDataType]]

Iterate over all elements.

__len__() → int

Return the number of elements.

__post_init__() → None

Allocate memory for values.

__setitem__(i: SupportsIndex, val: array_like, /) → None

Set the value for the current element.

com: ArrayCom

copy() → FlexibleElementArray[_ElementDataType, _IntegerDataType]

Create a copy of itself.

dtype: type[_ElementDataType]

resize_entry(i: int, new_size: int | Sequence[int], /) → None

Change the size of the entry, which also zeroes it.

shapes: FixedElementArray[_IntegerDataType]

unique(axis: int | None = None) → npt.NDArray[_ElementDataType]

Find unique values within the array.

values: tuple[ndarray[tuple[Any, ...], dtype[_ElementDataType]], ...]

The last subtype is ObjectElementArray, which is designed to hold one object of a specific type for each element. If the element entry is not set before accessing it, an exception is raised.

class mfv2d.element.ObjectElementArray(com: ArrayCom, dtype: type[_ElementObjectType], values: Iterable[_ElementObjectType] | None = None)

Element array that stores Python objects.

__delitem__(i: SupportsIndex, /) → None

Set the value for the current element.

__getitem__(i: SupportsIndex, /) → _ElementObjectType

Return element’s values.

__iter__() → Iterator[_ElementObjectType]

Iterate over all elements.

__len__() → int

Return the number of elements.

__setitem__(i: SupportsIndex, val: _ElementObjectType, /) → None

Set the value for the current element.

com: ArrayCom

dtype: type[_ElementObjectType]

values: list[_ElementObjectType | None]

Organizing Elements

Mesh elements are organized into an ElementCollection. This collection holds as ArrayCom object for all arrays related to this collection. Other important data it holds is related to parent-child relations, corners of individual elements and the orders of polynomials for these.

It also has some convenience methods to obtain indices of degrees of freedom on the boundary of an element, or an order of a specific element.

class mfv2d.element.ElementCollection(elements: Sequence[Element2D])

Element collection which contains relations and information about elements.

Parameters:

elements (Sequence of Element2D) – Sequence of elements which are to be used to construct the collection.

child_array: FlexibleElementArray[uint32, uint32]

Array which contains the indices of the children of each element. Number of children is stored in child_count_array.

child_count_array: FixedElementArray[uint32]

Array with the number of children of each element (shape (1,)). It is also the shapes array for the child_array array.

com: ArrayCom

Array coordination object.

corners_array: FixedElementArray[float64]

Array with the corners of the elements (shape (4, 2)).

orders_array: FixedElementArray[uint32]

Array with element orders in both directions (shape (2,)).

parent_array: FixedElementArray[uint32]

Indices of the parent elements (shape (1,)), which uses 1-based indexing, with the 0 value indicating that the element is a root element with no parent.

root_indices: ndarray[tuple[Any, ...], dtype[uint32]]

Array with indices of all top-level elements. Should be improved.

Executing per Element

One of key requirements for the solver is to be able to execute functions for each element. This is done by calling either the call_per_element_fix(), call_per_element_flex(), or call_per_leaf_obj() functions, which work for functions whose results are either fixed size, flexible size, or objects respectively.

mfv2d.element.call_per_element_fix(com: ~mfv2d.element.ArrayCom, dtype: type[~mfv2d.element._ElementDataType], shape: int | ~collections.abc.Sequence[int], fn: ~collections.abc.Callable[[~typing.Concatenate[int, ~_FuncArgs]], TypeAliasForwardRef('array_like')], *args: ~typing.~_FuncArgs, **kwargs: ~typing.~_FuncArgs) → FixedElementArray[_ElementDataType]

Call a function for each element and the return results.

Parameters:

• com (ArrayCom) – Array communication object.

• dtype (type) – Type of the result.

• shape (int of Sequence of int) – Shape of the results.

• fn (Callable) – Function to call for each element.

• *args – Arguments passed to the function.

• **kwargs – Keyword arguments passed to the function.

Returns:

Flexible array containing the results.

Return type:

FlexibleElementArray[dtype, np.uint32]

mfv2d.element.call_per_element_flex(com: ~mfv2d.element.ArrayCom, dims: int, dtype: type[~mfv2d.element._ElementDataType], fn: ~collections.abc.Callable[[~typing.Concatenate[int, ~_FuncArgs]], TypeAliasForwardRef('array_like')], *args: ~typing.~_FuncArgs, **kwargs: ~typing.~_FuncArgs) → FlexibleElementArray[_ElementDataType, uint32]

Call a function for each element and return the result as a flexible array.

Parameters:

• com (ArrayCom) – Array communication object.

• dims (int) – Number of dimensions of the result.

• dtype (type) – Type of the result.

• fn (Callable) – Function to call for each element.

• *args – Arguments passed to the function.

• **kwargs – Keyword arguments passed to the function.

Returns:

Flexible array containing the results.

Return type:

FlexibleElementArray[dtype, np.uint32]

mfv2d.element.call_per_leaf_obj(col: ~mfv2d.element.ElementCollection, dtype: type[~mfv2d.element._ElementObjectType], fn: ~collections.abc.Callable[[~typing.Concatenate[int, ~_FuncArgs]], ~mfv2d.element._ElementObjectType], *args: ~typing.~_FuncArgs, **kwargs: ~typing.~_FuncArgs) → ObjectElementArray[_ElementObjectType]

Call a function for each leaf element and return the result as a flexible array.

Note that this is only done for leaf elements (elements with no children).

Parameters:

• col (ElementCollection) – Element collection to use.

• dtype (type) – Type of the result.

• fn (Callable) – Function to call for each leaf element.

• *args – Arguments passed to the function.

• **kwargs – Keyword arguments passed to the function.

Returns:

Flexible array containing the results.

Return type:

ObjectElementArray[dtype]

There are also two specialized variant of the call_per_element_flex() function, named call_per_root_flex() and call_per_leaf_flex(). As can probably inferred from their names, these are called only for root or leaf elements.

mfv2d.element.call_per_root_flex(col: ~mfv2d.element.ElementCollection, dims: int, dtype: type[~mfv2d.element._ElementDataType], fn: ~collections.abc.Callable[[~typing.Concatenate[int, ~_FuncArgs]], TypeAliasForwardRef('array_like')], *args: ~typing.~_FuncArgs, **kwargs: ~typing.~_FuncArgs) → FlexibleElementArray[_ElementDataType, uint32]

Call a function for each root element and return the result as a flexible array.

Note that this is only done for root elements (elements with no parents).

Parameters:

• col (ElementCollection) – Element collection to use.

• dims (int) – Number of dimensions of the result.

• dtype (type) – Type of the result.

• fn (Callable) – Function to call for each leaf element.

• *args – Arguments passed to the function.

• **kwargs – Keyword arguments passed to the function.

Returns:

Flexible array containing the results.

Return type:

FlexibleElementArray[dtype, np.uint32]

mfv2d.element.call_per_leaf_flex(col: ~mfv2d.element.ElementCollection, dims: int, dtype: type[~mfv2d.element._ElementDataType], fn: ~collections.abc.Callable[[~typing.Concatenate[int, ~_FuncArgs]], TypeAliasForwardRef('array_like')], *args: ~typing.~_FuncArgs, **kwargs: ~typing.~_FuncArgs) → FlexibleElementArray[_ElementDataType, uint32]

Call a function for each leaf element and return the result as a flexible array.

Note that this is only done for leaf elements (elements with no children).

Parameters:

• col (ElementCollection) – Element collection to use.

• dims (int) – Number of dimensions of the result.

• dtype (type) – Type of the result.

• fn (Callable) – Function to call for each leaf element.

• *args – Arguments passed to the function.

• **kwargs – Keyword arguments passed to the function.

Returns:

Flexible array containing the results.

Return type:

FlexibleElementArray[dtype, np.uint32]

Element Utilities

There are some utility functions provided for computing number of elemet degrees of freedom or lagrange multipliers. These include their single-element variats _compute_element_dofs() as well as the exported compute_dof_sizes().

mfv2d.element._compute_element_dofs(ie: int, ordering: UnknownOrderings, elements: ElementCollection) → npt.NDArray[np.uint32]

Compute number of DoFs for each element.

Parameters:

• ie (int) – Index of the element.

• ordering (UnknownOrderings) – Orders of differential forms in the system.

• elements (ElementCollection) – Element collection to use.

Returns:

Array with count of degrees of freedom for of each differential form for the element.

Return type:

array of int

mfv2d.element.compute_dof_sizes(elements: ElementCollection, ordering: UnknownOrderings) → FixedElementArray[uint32]

Compute the number of DoFs for each element.

Parmeters

elementsElementCollection

Collection of elements for which to compute the sizes.

orderingUnknownOrderings

Orders of differential forms in the system.

returns:

Array with count of degrees of freedom for of each differential form for each element.

rtype:

FixedElementArray of np.uint32

Lagrange Multipliers and Constraints

Since continuity and strong boundary conditions are enforced throught Lagrange multipliers, there are types to represent these relations efficiently. The base building block is the ElementConstraint type, which describes the degrees of freedom and coefficients involved in a specific constraint. These may then be combined into Constraint type which associates one or more ElementConstraint with their right hand side value.

Note that ElementConstraint are also used to convey weak boundary condition information, however, in that case the coefficients in the ElementConstraint object represent contributions to the right hand side.

class mfv2d.element.ElementConstraint(i_e: int, dofs: npt.NDArray[np.uint32], coeffs: npt.NDArray[np.float64])

Type intended to enforce a constraint on an element.

Parameters:

• i_e (int) – Index of the element for which this constraint is applied.

• dofs ((n,) array) – Array with indices of the degrees of freedom of the element involved.

• coeffs ((n,) array) – Array with coefficients of degrees of freedom of the element involved.

coeffs: ndarray[tuple[Any, ...], dtype[float64]]

dofs: ndarray[tuple[Any, ...], dtype[uint32]]

i_e: int

class mfv2d.element.Constraint(rhs: float, *element_constraints: ElementConstraint)

Type used to specify constraints on degrees of freedom.

This type combines the individual ElementConstraint together with a right-hand side of the constraint.

Parameters:

• rhs (float) – The right-hand side of the constraint.

• *element_constraints (ElementConstraint) – Constraints to combine together.

element_constraints: tuple[ElementConstraint, ...]

rhs: float

Element Geometry

Some functions are also provided for computing geometry of elements. This includes computing the Jacobian matrix (with jacobian()), as well as computing the physical coordinates \((x, y)\) as functions of the reference domain coordinates \((\xi, \eta)\) (poly_x() and poly_y() respectively).

mfv2d.element.jacobian(corners: npt.NDArray[np.float64], nodes_1: npt.ArrayLike, nodes_2: npt.ArrayLike) → tuple[tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]], tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]]

Evaluate the Jacobian matrix entries.

The Jacobian matrix \(\mathbf{J}\) is defined such that:

\[\begin{split}\mathbf{J} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \\ \end{bmatrix}\end{split}\]

Which means that a coordinate transformation is performed by:

\[\begin{split}\begin{bmatrix} {dx} \\ {dy} \end{bmatrix} = \mathbf{J} \begin{bmatrix} {d\xi} \\ {d\eta} \end{bmatrix}\end{split}\]

Parameters:

• corners ((4, 2) array_like) – Array-like containing the corners of the element.

• nodes_1 (array_like) – The first computational component for the element where the Jacobian should be evaluated.

• nodes_2 (array_like) – The second computational component for the element where the Jacobian should be evaluated.

Returns:

• j00 (array) – The \((1, 1)\) component of the Jacobian corresponding to the value of \(\frac{\partial x}{\partial \xi}\).

• j01 (array) – The \((1, 2)\) component of the Jacobian corresponding to the value of \(\frac{\partial y}{\partial \xi}\).

• j10 (array) – The \((2, 1)\) component of the Jacobian corresponding to the value of \(\frac{\partial x}{\partial \eta}\).

• j11 (array) – The \((2, 2)\) component of the Jacobian corresponding to the value of \(\frac{\partial y}{\partial \eta}\).

mfv2d.element.poly_x(corner_x: npt.NDArray[np.float64], xi: npt.ArrayLike, eta: npt.ArrayLike) → npt.NDArray[np.float64]

Compute the x-coordiate of (xi, eta) points.

The relation for the x-coordinates is given by equation (1), where \(x_0\) is the bottom-left corner, \(x_1\) the bottom right, \(x_2\) the top right, and \(x_3\) is the top left corner.

(1)\[x(\xi, \eta) = x_0 \cdot \frac{(1 - \xi)(1 - \eta)}{4} + x_1 \cdot \frac{(1 + \xi)(1 - \eta)}{4} + x_2 \cdot \frac{(1 + \xi)(1 + \eta)}{4} + x_3 \cdot \frac{(1 - \xi)(1 + \eta)}{4}\]

Parameters:

• corner_x ((4,) array_like) – Array with the x-coordinates of the corners of the element.

• xi (array_like) – Array of the xi-coordinates on the reference domain where the x-coordinate should be evaluated.

• eta (array_like) – Array of the eta-coordinates on the reference domain where the x-coordinate should be evaluated.

Returns:

Array of x-coordinates evaluated at the given xi and eta points.

Return type:

array

mfv2d.element.poly_y(corner_y: npt.NDArray[np.float64], xi: npt.ArrayLike, eta: npt.ArrayLike) → npt.NDArray[np.float64]

Compute the y-coordiate of (xi, eta) points.

The relation for the y-coordinates is given by equation (2), where \(y_0\) is the bottom-left corner, \(x_1\) the bottom right, \(y_2\) the top right, and \(y_3\) is the top left corner.

(2)\[y(\xi, \eta) = y_0 \cdot \frac{(1 - \xi)(1 - \eta)}{4} + y_1 \cdot \frac{(1 + \xi)(1 - \eta)}{4} + y_2 \cdot \frac{(1 + \xi)(1 + \eta)}{4} + y_3 \cdot \frac{(1 - \xi)(1 + \eta)}{4}\]

Parameters:

• corner_y ((4,) array_like) – Array with the y-coordinates of the corners of the element.

• xi (array_like) – Array of the xi-coordinates on the reference domain where the x-coordinate should be evaluated.

• eta (array_like) – Array of the eta-coordinates on the reference domain where the y-coordinate should be evaluated.

Returns:

Array of y-coordinates evaluated at the given xi and eta points.

Return type:

array

Element Projection

This submodule also contains functions for projection of arbitrary functions on the element as either primal or dual degrees of freedom. Note that the dual projection is faster, since the primal has to be followed by a multiplication of an inverse of the mass matrix for the specific \(k\)-form.

mfv2d.element.element_primal_dofs(order: UnknownFormOrder, element_cache: ElementMassMatrixCache, function: Callable[[npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.ArrayLike]) → npt.NDArray[np.float64]

Compute the primal degrees of freedom of projection of a function on the element.

Primal degrees of freedom allow for reconstruction of the function’s \(L^2\) projection. This means that given these degrees of freedom, the following relation holds:

\[\left(\psi^{(k)}, f^{(k)} \right)_\Omega = \left(\psi^{(k)}, \bar{f}^{(k)} \right) = \int_{\Omega} \psi^{(k)} \sum\left( f_i \psi^{(k)}_i \right) dx \wedge dy\]

for all basis functions \(\psi^{(k)}\).

Parameters:

• order (UnknownFormOrder) – Order of the differential form to use for basis.

• corners (array_like) – Array of corners of the element.

• basis (Basis2D) – Basis function to use for the element.

• function (Callable) – Function to project.

Returns:

Array with primal degrees of freedom.

Return type:

array

mfv2d.element.element_dual_dofs(order: UnknownFormOrder, element_cache: ElementMassMatrixCache, function: Callable[[npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.ArrayLike]) → npt.NDArray[np.float64]

Compute the dual degrees of freedom (projection) of the function on the element.

Dual projection of a function is computed by taking an integral over the reference domain of the function with each of the two dimensional basis. This also means that these integrals are the values of the \(L^2\) projection of the function on the element.

For 0-forms are quite straight forward:

\[\left(f(x, y), \psi^{(0)}_i \right)_\Omega = \int_{\Omega} f(x, y) \psi^{(0)}_i(x, y) {dx} \wedge {dy} = \int_{\bar{\Omega}} f(x(\xi, \eta), y(\xi, \eta)) \psi^{(0)}_i(x(\xi, \eta), y(\xi, \eta)) \left| \mathbf{J} \right| d\xi \wedge d\eta\]

Similarly, 2-forms are also straight forward:

\[\left(f(x, y), \psi^{(2)}_i \right)_\Omega = \int_{\Omega} f(x, y) \psi^{(2)}_i(x, y) {dx} \wedge {dy} = \int_{\bar{\Omega}} f(x(\xi, \eta), y(\xi, \eta)) \psi^{(2)}_i(x(\xi, \eta), y(\xi, \eta)) \left| \mathbf{J} \right|^{-1} d\xi \wedge d\eta\]

Lastly, for 1-forms it is a bit more involved, since it has multiple components:

\[\left(f_x dy - f_y dx, {\psi^{(1)}_i}_x dy - {\psi^{(1)}_i}_y dx \right)_\Omega = \int_{\Omega} f_x {\psi^{(1)}_i}_x + f_y {\psi^{(1)}_i}_y {dx} \wedge {dy} = \int_{\bar{\Omega}} \left( \mathbf{J} \vec{f} \right) \cdot \left( \mathbf{J} \vec{\psi^{(1)}} \right) d\xi \wedge d\eta\]

Parameters:

• order (UnknownFormOrder) – Order of the differential form to use for basis.

• corners (array_like) – Array of corners of the element.

• basis (Basis2D) – Basis function to use for the element.

• function (Callable) – Function to project.

Returns:

Array with dual degrees of freedom.

Return type:

array

Reconstruction

To be able to return result for a solver as VTK file or even for being able to compute interior product, there has to be functionality to compute pointwise values for a \(k\)-form given its element degrees of freedom. This is handled by the reconstruct() function.

mfv2d.element.reconstruct(corners: npt.NDArray[np.float64], k: UnknownFormOrder, dofs: npt.ArrayLike, xi: npt.ArrayLike, eta: npt.ArrayLike, basis: Basis2D) → npt.NDArray[np.float64]

Reconstruct a k-form on the element from its primal degrees of freedom.

Parameters:

• corners (array_like) – Array of corners of the element.

• k (UnknownFormOrder) – Order of the differential form to use for basis.

• dofs (array_like) – Degrees of freedom of the k-form.

• xi (array_like) – Coordinates of the xi-coordinate in the reference domain.

• eta (array_like) – Coordinates of the eta-coordinate in the reference domain.

• basis (Basis2D) – Basis function to use for the element.

Returns:

Array with the point values of the k-form at the specified coordinates.

Return type:

array