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.

Mesh Geometry

Mesh geometry is handled by the Mesh2D type. This uses mfv2d._mfv2d.Manifold2D to store its topology and its dual topology, but also contains geometry. It is also equiped with helper methods for obtaining elements and plotting with pyvista.

class mfv2d.mimetic2d.Mesh2D(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)

Two dimensional manifold with associated geometry.

Mesh holds the primal manifold, which describes the topology of surfaces and lines that make it up. It also contains the dual mesh, which contains duals of all primal geometrical objects. The dual is useful when connectivity is needed.

Parameters:

• order (int or Sequence of int or array-like) – Orders of elements. If a single value is specified, then the same order is used for all elements. If a sequence is specified, then each element is given that order.

• positions ((N, 2) array-like) – Positions of the nodes.

• lines ((N, 2) array-like) – Lines of the mesh specified as pairs of nodes connected. These use 1-based indexing.

• surfaces ((N, 4) array-like) – Surfaces of the mesh specified in their positive orientation as 1-based indices of lines that make up the surface, with negative sign indicating reverse direction.

boundary_indices: ndarray[tuple[Any, ...], dtype[int32]]

Indices of lines that make up the boundary of the mesh. These can be useful when prescribing boundary conditions.

dual: Manifold2D

Dual topology of the mesh.

property n_elements: int

Number of (surface) elements in the mesh.

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

Orders of individual elements.

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

Array of positions for each node in the mesh.

primal: Manifold2D

Primal topology of the mesh.

surface_to_element(idx: int, /) → ElementLeaf2D

Create a 2D element from a surface with the given index.

Parameters:

idx (int) – Index of the surface.

Returns:

Leaf element with the geometry of the specified surface.

Return type:

ElementLeaf2D

Elements

Surfaces of the mesh can be converted into computational elements. These can either be leaf elements, with the type ElementLeaf2D, or nodes, which contain other nodes or leaves, with the type ElementNode2D. Both are subtypes of Element2D, which only contains the information about the parent.

class mfv2d.mimetic2d.ElementLeaf2D(parent: ElementNode2D | None, order: int, bottom_left: tuple[float, float], bottom_right: tuple[float, float], top_right: tuple[float, float], top_left: tuple[float, float])

Two dimensional square element.

This type facilitates operations related to calculations which need to be carried out on the reference element itself, such as calculation of the mass and incidence matrices, as well as the reconstruction of the solution.

Parameters:

• order_h (int) – Order of the basis functions used for the nodal basis in the first dimension.

• order_v (int) – Order of the basis functions used for the nodal basis in the second dimension.

• bottom_left ((float, float)) – Coordinates of the bottom left corner.

• bottom_right ((float, float)) – Coordinates of the bottom right corner.

• top_right ((float, float)) – Coordinates of the top right corner.

• top_left ((float, float)) – Coordinates of the top left corner.

bottom_left: tuple[float, float]

bottom_right: tuple[float, float]

divide(order_bl: int, order_br: int, order_tl: int, order_tr: int) → tuple[ElementNode2D, tuple[tuple[ElementLeaf2D, ElementLeaf2D], tuple[ElementLeaf2D, ElementLeaf2D]]]

Divide the element into four child elements of the specified order.

Parameters:

• order_bl (int) – Order of the bottom left element.

• order_br (int) – Order of the bottom right element.

• order_tl (int) – Order of the top left element.

• order_tr (int) – Order of the top right element.

Returns:

• ElementNode2D – Parent element which contains the nodes.

• (2, 2) tuple of ElementLeaf2D – Child elements of the same order as the element itself. Indexing the tuple will give bottom/top for the first axis and left/right for the second.

order: int

parent: ElementNode2D | None

top_left: tuple[float, float]

top_right: tuple[float, float]

class mfv2d.mimetic2d.ElementNode2D(parent: ElementNode2D | None, child_bl: Element2D, child_br: Element2D, child_tl: Element2D, child_tr: Element2D)

Two dimensional element that contains children.

child_bl: Element2D

child_br: Element2D

child_tl: Element2D

child_tr: Element2D

children() → tuple[Element2D, Element2D, Element2D, Element2D]

Return children of the element ordered.

parent: ElementNode2D | None

class mfv2d.mimetic2d.Element2D(parent: ElementNode2D | None)

General 2D element.

parent: ElementNode2D | None

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