Mimetic interpolation

In mathematics, mimetic interpolation is a method for interpolating differential forms. In contrast to other interpolation methods, which estimate a field at a location given its values on neighboring points, mimetic interpolation estimates the field's [math]\displaystyle{ k }[/math]-form given the field's projection on neighboring grid elements. The grid elements can be grid points as well as cell edges or faces, depending on [math]\displaystyle{ k = 0, 1, 2, \cdots }[/math].

Mimetic interpolation is particularly relevant in the context of vector and pseudo-vector fields as the method conserves line integrals and fluxes, respectively.

Interpolation of integrated forms

Let [math]\displaystyle{ \omega^k }[/math] be a differential [math]\displaystyle{ k }[/math]-form, then mimetic interpolation is the linear combination

[math]\displaystyle{ \bar{\omega}^k = \sum_i \left(\int_{M_i} \omega^k \right) \phi_i^k }[/math]

where [math]\displaystyle{ \bar{\omega}^k }[/math] is the interpolation of [math]\displaystyle{ \omega^k }[/math], and the coefficients [math]\displaystyle{ \int_{M_i} \omega^k }[/math] represent the strengths of the field on grid element [math]\displaystyle{ M_i }[/math]. Depending on [math]\displaystyle{ k }[/math], [math]\displaystyle{ M_i }[/math] can be a node ([math]\displaystyle{ k = 0 }[/math]), a cell edge ([math]\displaystyle{ k=1 }[/math]), a cell face ([math]\displaystyle{ k=2 }[/math]) or a cell volume ([math]\displaystyle{ k=3 }[/math]). In the above, the [math]\displaystyle{ \phi_i^k }[/math] are the interpolating [math]\displaystyle{ k }[/math]-forms, which are centered on [math]\displaystyle{ M_i }[/math] and decay away from [math]\displaystyle{ M_i }[/math] in a way similar to the tent functions. Examples of [math]\displaystyle{ \phi_i^k }[/math] are the Whitney forms for simplicial meshes in [math]\displaystyle{ n }[/math] dimensions.

An important advantage of mimetic interpolation over other interpolation methods is that the field strengths [math]\displaystyle{ \int_{M_i} \omega^k }[/math] are scalars and thus coordinate system invariant.

Interpolating forms

In many cases it is desirable for the interpolating forms [math]\displaystyle{ \phi_i^k }[/math] to pick the field's strength on particular grid elements without interference from other [math]\displaystyle{ \phi_j^k }[/math]. This allows one to assign field values to specific grid elements, which can then be interpolated in-between. A common case is linear interpolation for which the interpolating functions ([math]\displaystyle{ 0 }[/math]-forms) are zero on all nodes except on one, where the interpolating function is one. A similar construct can be applied to mimetic interpolation

[math]\displaystyle{ \int_{M_j} \phi_i^k = \delta_{ij}. }[/math]

That is, the integral of [math]\displaystyle{ \phi_i^k }[/math] is zero on all cell elements, except on [math]\displaystyle{ M_i }[/math] where the integral returns one. For [math]\displaystyle{ k=0 }[/math] this amounts to [math]\displaystyle{ \phi_i^0(x_j) = \delta_{ij} }[/math] where [math]\displaystyle{ x_j }[/math] is a grid point. For [math]\displaystyle{ k = 1 }[/math] the integral is over edges and hence the integral [math]\displaystyle{ \int_{M_j} \phi_i^1 }[/math] is zero expect on edge [math]\displaystyle{ M_i }[/math]. For [math]\displaystyle{ k=2 }[/math] the integral is over faces and for [math]\displaystyle{ k = 3 }[/math] over cell volumes.

Conservation properties

Mimetic interpolation respects the properties of differential forms. In particular, Stokes' theorem

[math]\displaystyle{ \int_M \overline{d \omega^k} = \int_{\partial M} \bar{\omega}^k }[/math]

is satisfied with [math]\displaystyle{ \overline{d \omega^k} }[/math] denoting the interpolation of [math]\displaystyle{ d \omega^k }[/math]. Here, [math]\displaystyle{ d }[/math] is the exterior derivative, [math]\displaystyle{ M }[/math] is any manifold of dimensionality [math]\displaystyle{ k }[/math] and [math]\displaystyle{ \partial M }[/math] is the boundary of [math]\displaystyle{ M }[/math]. This confers to mimetic interpolation conservation properties, which are not generally shared by other interpolation methods.

Commutativity between the interpolation operator and the exterior derivative

De Rham complex. Top: the spaces of differential forms in three dimensions. Bottom: the corresponding discretized versions of the differential forms obtained after interpolation. The commutativity condition ensures that the dashed and dash-dotted paths give the same result for each de Rham cell.

To be mimetic, the interpolation must satisfy

[math]\displaystyle{ I_{k+1} ( d \omega^k ) = d (I_k \omega^k) }[/math]

where [math]\displaystyle{ I_k }[/math] is the interpolation operator of a [math]\displaystyle{ k }[/math]-form, i.e. [math]\displaystyle{ \bar{\omega}^k = I_k \omega^k }[/math]. In other words, the interpolation operators and the exterior derivatives commute. Note that different interpolation methods are required for each type of form ([math]\displaystyle{ k = 0, 1, \cdots }[/math]), [math]\displaystyle{ I_{k+1} \neq I_k }[/math]. The above equation is all that is needed to satisfy Stokes' theorem for the interpolated form

[math]\displaystyle{ \int_M \overline{d \omega^k } = \int_M I_{k+1} ( d \omega^k ) = \int_M d (I_k \omega^k) = \int_{\partial M} I_k \omega^k = \int_{\partial M} \bar{\omega}^k. }[/math]

Other calculus properties derive from the commutativity between interpolation and [math]\displaystyle{ d }[/math]. For instance, [math]\displaystyle{ d^2 \omega = 0 }[/math],

[math]\displaystyle{ \int_M I_{k+2} (d^2 \omega^k) = \int_M d (I_{k+1} d \omega^k) = \int_{\partial M} I_{k+1} d \omega^k = \int_{\partial M} d (I_k \omega^k) = 0. }[/math]

The last step gives zero since [math]\displaystyle{ \int_{\partial M} d(\cdot) = 0 }[/math] when integrated over the boundary [math]\displaystyle{ \partial M }[/math].

Projection

The interpolated [math]\displaystyle{ \bar{\omega}^k }[/math] is often projected onto a target, [math]\displaystyle{ k }[/math]-dimensional, oriented manifolds [math]\displaystyle{ T }[/math],[math]\displaystyle{ \int_T \bar{\omega}^k = \sum_i \left( \int_{M_i} \omega^k \right) \left(\int_T \phi_i^k \right). }[/math]For [math]\displaystyle{ k=0 }[/math] the target is a point, for [math]\displaystyle{ k=1 }[/math] it is a line, for [math]\displaystyle{ k=2 }[/math] an area, etc.

Applications

Many physical fields can be represented as [math]\displaystyle{ k }[/math]-forms. When discretizing fields in numerical modeling, each type of [math]\displaystyle{ k }[/math]-form often acquires its own staggering in accordance with numerical stability requirements, e.g. the need to prevent the checkerboard instability. This led to the development of the exterior finite element and discrete exterior calculus methods, both of which rely on a field discretization that are compatible with the field type.

The table below lists some examples of physical fields, their type, their corresponding form and interpolation method, as well as software that can be leveraged to interpolate, remap or regrid the field:

field example	field type	k-form equivalent	target	staggering	Interpolation method (example)	example of software
temperature	scalar	0-form	point	nodal	bilinear, trilinear	ESMF
electric field	vector	1-form	line	edge	edge	MINT
magnetic field	pseudo-vector	2-form	area	face	face	MINT
density	pseudo-scalar	3-form	volume	cell	area weighted, conservative	SCRIP, ESMF

Example

Indexing of nodes and edges for a quadrilateral cell, as used in the example. The edges are chosen to point in the east and north directions.

Consider quadrilateral cells in two dimensions with their node indexed [math]\displaystyle{ i = 0, 1, 2, 3 }[/math] in the counterclockwise direction. Further, let [math]\displaystyle{ \xi_1 }[/math] and [math]\displaystyle{ \xi_2 }[/math] be the parametric coordinates of each cell ([math]\displaystyle{ 0 \leq \xi_1, \xi_2 \leq 1 }[/math]). Then[math]\displaystyle{ \begin{align} \phi_0^0 & = (1-\xi_1) (1-\xi_2) \\[0.6ex] \phi_1^0 & = \xi_1 (1-\xi_2) \\[0.6ex] \phi_2^0 & = \xi_1 \xi_2 \\[0.6ex] \phi_3^0 & = (1-\xi_1) \xi_2. \end{align} }[/math]

are the bilinear interpolating forms of [math]\displaystyle{ I_0 }[/math] in the unit square ([math]\displaystyle{ \xi_1, \xi_2 }[/math]). The corresponding [math]\displaystyle{ I_1 }[/math] edge interpolating forms are

[math]\displaystyle{ \begin{align} \phi_0^1 & = (1-\xi_2) d\xi_1 \\[0.6ex] \phi_1^1 & = \xi_1 d\xi_2 \\[0.6ex] \phi_2^1 & = \xi_2 d\xi_1 \\[0.6ex] \phi_3^1 & = (1-\xi_1) d\xi_2, \end{align} }[/math]

The four vector fields which are dual to the four edge interpolating forms attached to a quadrilateral cell. The vector fields are strongest on their supporting edge and decrease to zero towards the opposite edge. Note how the vector fields bend to enforce perpendicularity with respect to edges that are adjacent to the supporting edge.

were we assumed the edges to be indexed in counterclockwise direction and with the edges pointing to the east and north. At lowest order, there is only one interpolating form for [math]\displaystyle{ I_2 }[/math],

[math]\displaystyle{ \phi_0^2 = d\xi_1 \wedge d\xi_2, }[/math]

where [math]\displaystyle{ \wedge }[/math] is the wedge product.

We can verify that the above interpolating forms satisfy the mimetic conditions [math]\displaystyle{ I_1 d\omega^0 = d(I_0 \omega^0) }[/math] and [math]\displaystyle{ d(I_1 \omega^1) = I_2 d \omega^1 }[/math]. Take for instance,

[math]\displaystyle{ \begin{align} d(I_0 \omega^0) & = d \left(f_0 (1-\xi_1) (1-\xi_2)+f_1 \xi_1 (1-\xi_2) + f_2 \xi_1 \xi_2 + f_3 (1-\xi_1) \xi_2 \right) \\[0.6ex] & = (f_1-f_0)(1-\xi_2)d\xi_1 + (f_2-f_1) \xi_1 d\xi_2 + (f_2-f_3) \xi_2 d\xi_1 + (f_3-f_0) (1-\xi_1)d\xi_2 \\[0.6ex] & = (f_1-f_0)\phi_0^1 + (f_2-f_1)\phi_1^1 + (f_2-f_3)\phi_2^1 + (f_3-f_0)\phi_3^1 \\[0.6ex] & = I_1 d\omega^0 \end{align} }[/math]where [math]\displaystyle{ f_0 = \omega^0(0,0) }[/math], [math]\displaystyle{ f_1 = \omega^0(1,0) }[/math], [math]\displaystyle{ f_2 = \omega^0(1,1) }[/math] and [math]\displaystyle{ f_3 = \omega^0(0,1) }[/math] are the field values evaluated at the corners of the quadrilateral in the unit square space. Likewise, we have

[math]\displaystyle{ \begin{align} d(I_1 \omega^1) & = d( g_0(1-\xi_2)d\xi_1 + g_1 \xi_1 d\xi_2 + g_2 \xi_2 d\xi_1 + g_3 (1-\xi_1)d\xi_2) \\[0.6ex] & = (g_0 + g_1 - g_2 - g_3) d\xi_1 \wedge d\xi_2 \\[0.6ex] & = I_2 d \omega^1 \end{align} }[/math] with [math]\displaystyle{ g_i = \int_{M_i} \omega^1 }[/math], [math]\displaystyle{ i \in \{0, 1, 2, 3\} }[/math], being the 1-form [math]\displaystyle{ \omega^1 }[/math]projected onto edge [math]\displaystyle{ i }[/math]. Note that [math]\displaystyle{ g_i }[/math] is also known as the pullback. If [math]\displaystyle{ F: \mathbb{R} \rightarrow \mathbb{R}^2 }[/math] is the map that parametrizes edge [math]\displaystyle{ i }[/math], [math]\displaystyle{ x = F_i(t) }[/math], [math]\displaystyle{ 0 \leq t \leq 1 }[/math], then [math]\displaystyle{ g_i = \int {F_i}^* \omega^1 }[/math]where the integration is performed in [math]\displaystyle{ t }[/math] space. Consider for instance edge [math]\displaystyle{ i=2 }[/math] , then [math]\displaystyle{ F_2(t) = x_3 + t(x_2 - x_3) }[/math] with [math]\displaystyle{ x_2 }[/math] and [math]\displaystyle{ x_3 }[/math] denoting the start and points. For a general 1-form [math]\displaystyle{ \omega^1 = a(x) d\xi_1 + b(x) d \xi_2 }[/math], one gets [math]\displaystyle{ g_2 = \int_0^1 \left( a(t) \frac{\partial \xi_1}{\partial t} + b(t) \frac{\partial \xi_2}{\partial t} \right) dt }[/math].