mimetic2d

Most of the module is in the process of getting refactored into others. Previously, code was written with element objects, which were used to hold element relations. These are now all handled by mfv2d.element.ElementCollection.

Incidence Functions

While these are not used internally, when it comes to testing these functions play a very important role. They are not fast nor are they efficient, their main purpose is to be used in tests to validate that operations done in C are correct. If they were to be rewritten to allow different polynomial order in each dimension, it would allow the code to deal with mixed-order elements.

Generating Matrices

These directly generate the incidence matrices.

mfv2d.mimetic2d.incidence_10(order: int) → npt.NDArray[np.float64]

Incidence matrix from 0.forms to 1-forms.

This applies the exterior derivative operation to primal 0-forms and maps them into 1-forms. The negative transpose is the equivalent operation for the dual 1-forms, the derivatives of which are consequently dual 2-forms.

This is done by mapping degrees of freedom of the original primal 0-form or dual 1-form into those of the derivative primal 1-forms or dual 2-forms respectively.

\[\vec{\mathcal{N}}^{(1)}(f) = \mathbb{E}^{(1,0)} \vec{\mathcal{N}}^{(0)}(f)\]

\[\tilde{\mathcal{N}}^{(2)}(f) = -\left(\mathbb{E}^{(1,0)}\right)^{T} \tilde{\mathcal{N}}^{(1)}(f)\]

Returns:

Incidence matrix \(\mathbb{E}^{(1,0)}\).

Return type:

array

mfv2d.mimetic2d.incidence_21(order: int) → npt.NDArray[np.float64]

Incidence matrix from 1-forms to 2-forms.

This applies the exterior derivative operation to primal 1-forms and maps them into 2-forms. The negative transpose is the equivalent operation for the dual 0-forms, the derivatives of which are consequently dual 1-forms.

This is done by mapping degrees of freedom of the original primal 1-form or dual 0-form into those of the derivative primal 2-forms or dual 1-forms respectively.

\[\vec{\mathcal{N}}^{(2)}(f) = \mathbb{E}^{(2,1)} \vec{\mathcal{N}}^{(1)}(f)\]

\[\tilde{\mathcal{N}}^{(1)}(f) = -\left(\mathbb{E}^{(2,1)}\right)^{T} \tilde{\mathcal{N}}^{(0)}(f)\]

Returns:

Incidence matrix \(\mathbb{E}^{(2,1)}\).

Return type:

array

Applying Matrices

Since C code does not explicitly compute incidnece matrices unless necessary, these are a lot more commonly used. Here they are written very explicitly, so that they can be translated into C almost line by line.

mfv2d.mimetic2d.apply_e10(order: int, other: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Apply the E10 matrix to the given input.

Calling this function is equivalent to left multiplying by E10.

mfv2d.mimetic2d.apply_e10_t(order: int, other: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Apply the E10 transpose matrix to the given input.

Calling this function is equivalent to left multiplying by E10 transposed.

mfv2d.mimetic2d.apply_e10_r(order: int, other: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Apply the right E10 matrix to the given input.

Calling this function is equivalent to right multiplying by E10.

mfv2d.mimetic2d.apply_e10_rt(order: int, other: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Apply the right transposed E10 matrix to the given input.

Calling this function is equivalent to right multiplying by E10 transposed.

mfv2d.mimetic2d.apply_e21(order: int, other: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Apply the E21 matrix to the given input.

Calling this function is equivalent to left multiplying by E21.

mfv2d.mimetic2d.apply_e21_t(order: int, other: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Apply the E21 transposed matrix to the given input.

Calling this function is equivalent to left multiplying by E21 transposed.

mfv2d.mimetic2d.apply_e21_r(order: int, other: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Apply the right E21 matrix to the given input.

Calling this function is equivalent to right multiplying by E21.

mfv2d.mimetic2d.apply_e21_rt(order: int, other: npt.NDArray[np.float64]) → npt.NDArray[np.float64]

Apply the right transpose E21 matrix to the given input.

Calling this function is equivalent to right multiplying by E21 transposed.

Element Sides

Many functions have to perform operations on boundaries of the elements. To make this easier to type check, the type ElementSide is introduced. This is an enum.IntEnum subtype with only four values and is used everywhere when a side of an element is an input of any function.

class mfv2d.mimetic2d.ElementSide(*values)

Enum specifying the side of an element.

SIDE_BOTTOM = 1

SIDE_RIGHT = 2

SIDE_TOP = 3

SIDE_LEFT = 4

property next: ElementSide

Next side.

property prev: ElementSide

Previous side.

To help identify what side a Line with a specific index of a Surface is from, the function find_surface_boundary_id_line() is given.

mfv2d.mimetic2d.find_surface_boundary_id_line(s: Surface, i: int) → ElementSide

Find what boundary the line with a given index is in the surface.

Caching Basis

Since only a handful of basis and integration orders are ever used in the solve and creating mfv2d._mfv2d.Basis1D and mfv2d._mfv2d.IntegrationRule1D both is not extremely cheap, the FemCache is introduced to deal with caching them. It does not deal with caching mfv2d._mfv2d.Basis2D objects, since they are just containers for two objects and so are not a significant cost to newly construct each time.

class mfv2d.mimetic2d.FemCache(order_difference: int)

Cache for integration rules and basis functions.

This type allows for caching 1D integration rules and 1D basis. The 2D basis are not cached, since the Basis2D is just a container for two 1D Basis objects, so it is probably just as cheap to create a new one as it would be to cache it.

Parameters:

order_difference (int) – Difference of orders between the integration rule and the basis in the case that it is not specified.

clean() → None

Clear all caches.

get_basis1d(order: int, int_order: int | None = None) → Basis1D

Get requested one-dimensional basis.

Parameters:

• order (int) – Order of the basis.

• int_order (int, optional) – Order of the integration rule for the basis. If it is not specified, then order + self.order_diff is used.

Returns:

One-dimensional basis.

Return type:

Basis1D

get_basis2d(order1: int, order2: int, int_order1: int | None = None, int_order2: int | None = None) → Basis2D

Get two-dimensional basis.

These are not cached, since there is zero calculations involved in their computations.

Parameters:

• order1 (int) – Order of basis in the first dimension.

• order2 (int) – Order of basis in the second dimension.

• int_order1 (int, optional) – Order of the integration rule in the first direction. If unspecified, it defaults to order1 + self.order_diff.

• int_order2 (int, optional) – Order of the integration rule in the second direction. If unspecified, it defaults to order2 + self.order_diff.

Returns:

Requested two-dimensional basis.

Return type:

Basis2D

get_integration_rule(order: int) → IntegrationRule1D

Return integration rule.

Parameters:

order (int) – Order of integration rule to use.

Returns:

Integration rule that was obtained from cache or created and cached if it was previously not there.

Return type:

IntegrationRule1D

get_mass_inverse_1d_edge(order: int) → npt.NDArray[np.float64]

Get the 1D edge mass matrix inverse.

get_mass_inverse_1d_node(order: int) → npt.NDArray[np.float64]

Get the 1D nodal mass matrix inverse.

order_diff: int

Mesh Geometry

Since mfv2d._mfv2d.Mesh is meant as a container, constructing it is repetative and cumbersome. As such, the function mesh_create() is intended to be used for to create a new mesh from position and connectivity data.

mfv2d.mimetic2d.mesh_create(order: int | Sequence[int] | npt.ArrayLike, positions: Sequence[tuple[float, float, float]] | Sequence[Sequence[float]] | Sequence[npt.ArrayLike] | npt.ArrayLike, lines: Sequence[tuple[int, int]] | Sequence[npt.ArrayLike] | Sequence[Sequence[int]] | npt.ArrayLike, surfaces: Sequence[tuple[int, ...]] | Sequence[Sequence[int]] | Sequence[npt.ArrayLike] | npt.ArrayLike) → Mesh

Create new mesh from given geometry.

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

Degree of Freedom Counts

There are some utility functions provided for computing number of elemet degrees of freedom or lagrange multipliers. This is done by compute_leaf_dof_counts().

mfv2d.mimetic2d.compute_leaf_dof_counts(order_1: int, order_2: int, ordering: UnknownOrderings) → npt.NDArray[np.uint32]

Compute number of DoFs for each element.

Parameters:

• order_1 (int) – Order of the element in the first dimension.

• order_2 (int) – Order of the element in the second dimension.

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

Returns:

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

Return type:

array of int

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)\), since it is assumed these are bilinear (hence the bilinear_interpolate() function).

mfv2d.mimetic2d.jacobian(corners: npt.NDArray[np.floating], 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.mimetic2d.bilinear_interpolate(corner_vals: npt.NDArray[np.floating], xi: npt.ArrayLike, eta: npt.ArrayLike) → npt.NDArray[np.float64]

Compute bilinear interpolation at (xi, eta) points.

The relation for the bilinear interpolation is given by equation (1), where \(u_0\) is the bottom-left corner, \(u_1\) the bottom right, \(u_2\) the top right, and \(u_3\) is the top left corner.

(1)\[u(\xi, \eta) = u_0 \cdot \frac{(1 - \xi)(1 - \eta)}{4} + u_1 \cdot \frac{(1 + \xi)(1 - \eta)}{4} + u_2 \cdot \frac{(1 + \xi)(1 + \eta)}{4} + u_3 \cdot \frac{(1 - \xi)(1 + \eta)}{4}\]

Parameters:

• corner_vals ((4,) array_like) – Array with the values in 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 interpolated values 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.mimetic2d.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.mimetic2d.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. To help with it, vtk_lagrange_ordering() is provided, since VTK does not like sensible order of unknowns.

mfv2d.mimetic2d.reconstruct(corners: npt.NDArray[np.floating], 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

mfv2d.mimetic2d.vtk_lagrange_ordering(order: int) → npt.NDArray[np.uint32]

Ordering for vtkLagrangeQuadrilateral.

VTK has an option to create cells of type LagrangeQuadrilateral. These allow for arbitrary order of interpolation with nodal basis. Due to backwards compatibility the ordering of the nodes in these is done in an unique way. As such, either the positions or ordering of the nodes must be adjusted.

This function returns the correct order which can be used for either given a specific polynomial order.

Parameters:

order (int) – Order of the element.

Returns:

Array of indices which correctly order nodes on an element of the specified order.

Return type:

array