Systematic implementation of higher order Whitney forms in methods based on discrete exterior calculus

We present a systematic way to implement higher order Whitney forms in numerical methods based on discrete exterior calculus. Given a simplicial mesh, we first refine the mesh into smaller simplices which can be used to define higher order Whitney forms. Cochains on this refined mesh can then be interpolated using higher order Whitney forms. Hence, when the refined mesh is used with methods based on discrete exterior calculus, the solution can be expressed as a higher order Whitney form. We present algorithms for the three required steps: refining the mesh, solving the coefficients of the interpolant, and evaluating the interpolant at a given point. With our algorithms, the order of the Whitney forms one wishes to use can be given as a parameter so that the same code covers all orders, which is a significant improvement on previous implementations. Our algorithms are applicable with all methods in which the degrees of freedom are integrals over mesh simplices — that is, when the solution is a cochain on a simplicial mesh. They can also be used when one simply wishes to approximate differential forms in finite-dimensional spaces. Numerical examples validate the generality of our algorithms.


Introduction
Partial differential equations describing field theories such as electromagnetism and elasticity often admit a natural expression in terms of differential forms. In the past decades, the role of differential forms has increased also in numerical methods and their practical implementations. The growth of their popularity has been accelerated by finite element exterior calculus [1,3] for the finite element method and by discrete exterior calculus [13,16] for finite difference kind of methods. For computations, differential forms are approximated in finite-dimensional spaces using a mesh consisting of a finite number of cells. In finite element exterior calculus these finitedimensional spaces are spanned by suitable finite elements, while discrete exterior calculus is based on cochains (also known as discrete forms) as approximations for differential forms.
Methods based on discrete exterior calculus (DEC) may go by different names, e.g. Yee-like schemes [9,10], finite integration technique [12], or generalised finite differences [6,8]. These methods enable one to distinguish the features that depend on metric from those that do not, and in a way they preserve the geometric structure of the continuous model at the discrete level. The methods can typically be made explicit, which enables a very large number of degrees of freedom. In the literature there are many examples where discrete exterior calculus has been successfully applied (see e.g. [17,[22][23][24][25][26]29]).
Although methods based on cochains have their benefits, there are also some drawbacks. When the solution is given as a cochain, it cannot be evaluated at a given point. In some situations evaluating the solution at a given point is preferable, and then one has to interpolate the cochain somehow. This raises the question: in which space does the interpolant lie? In the case of simplicial meshes, it is well known that Whitney forms [31] can be used to approximate differential forms and interpolate cochains in methods based on discrete exterior calculus. However, this only applies to lowest order Whitney forms. Although higher order Whitney forms have been defined and used elsewhere [4,5,7,11,15,21,27,28], they have not been used with cochainbased methods. Indeed, these are often considered low-order methods in that they seem to lack natural higher order generalisations.
In this paper, we provide an alternative viewpoint and show how higher order Whitney forms can be used to interpolate cochains in methods based on discrete exterior calculus. As with lowest order Whitney forms, we require a simplicial mesh to begin with. This mesh is refined into smaller simplices which have been used to define higher order Whitney forms in [28]. Cochains on this refined mesh can then be interpolated using higher order Whitney forms, so when we apply methods based on discrete exterior calculus with the refined mesh, the solution can be expressed as a higher order Whitney form.
With this approach, we can reduce the interpolation error without any modifications in the methods themselves; the only changes are in preprocessing (preparing the mesh) and postprocessing (interpolating the cochain) stages. For this reason, we do not focus on any specific method here, but instead provide a framework for interpolating cochains with higher order Whitney forms. This framework can then be applied with any method -our algorithms are applicable whenever the solution is a cochain on a simplicial mesh. They can also be used if one simply wishes to approximate differential forms in finite-dimensional spaces and might be relevant for the finite element method as well, if integrals over small simplices are chosen as degrees of freedom.
Although reducing the interpolation error alone does not lead to higher order DEC methods, it is a necessary step toward them. To obtain higher order convergence, one would also have to improve the accuracy of discrete Hodge operators. This is possible using higher order Whitney forms and the interpolation framework presented here. It would of course require changes in specific methods, and hence, we will study higher order discrete Hodge operators in a future article.
The main novelty of this paper is the systematic implementation strategy that is usable with DEC and yields Whitney forms of all orders with the same code. The idea of using small simplices to construct bases and degrees of freedom for higher order Whitney forms is not new [11,28]; however, although systematic implementations for more traditional bases and dofs exist, the approach with small simplices still lacks a systematic implementation strategy. The implementation in 3D is a delicate issue. For FEM, it has been studied in [4], but only second-and third-order 1-forms were implemented. Similarly, only second-order forms were implemented in our previous studies for DEC [18,19]. Without general algorithms, the implementation process has to be repeated separately for each order, while the workload becomes unreasonably laborious very quickly. The systematic implementation strategy of this paper and in particular Algorithms 1-3 are a significant novelty, yielding all orders with the same code. Without such algorithms, the implementation of, for example, 8th-order 1-forms in 3D would be practically impossible using small simplices (to define both basis functions and dofs).
The outline of this paper is as follows. We start with some preliminaries in Section 2. Section 3 covers lowest order Whitney forms, and in Section 4 we recall the small simplices of [28] and use them to define higher order Whitney forms. In Section 5 we give the general idea for interpolating cochains with higher order Whitney forms. The implementation of this idea is discussed in Section 6. This section contains three subsections where we present algorithms for the three required steps: refining the mesh, solving the coefficients of the interpolant, and evaluating the interpolant at a given point. Numerical examples of Section 7 validate the generality of our algorithms -the order of Whitney forms can be given as a parameter, and hence, all cases are covered by the same code.

Some preliminary concepts
In this section we recall some prerequisite concepts that are used in this paper. The discussion is brief, but readers unfamiliar with these concepts can consult the given references for more information.
We start by defining a mesh in a domain ⊂ R n . Bounded and convex pdimensional polytopes in R n are called p-cells for short. Cell complex K is a finite set of cells such that • each face of every cell in K is also in K.
• The intersection of two cells in K is either a common face of theirs or the empty set.
The set of p-cells in K is denoted by S p (K). A cell complex K is a mesh in if the union of the cells of K is . To enable this we assume that our domain is a bounded polytope in R n . A cell complex (or mesh) is simplicial if its cells are all simplices. In this case there is a unique barycentric function [31, App. II, §2] corresponding to each 0-simplex x i -this is denoted by λ i .
We will also need the concept of orientation [31, App. II, §5]. Recall that p-cell σ is oriented by orienting its plane. If τ ∈ S p+1 (K) is oriented and σ ∈ S p (K) is a face of τ , the orientation of τ induces an orientation on σ ; we say that the orientation of σ agrees with that of τ if σ is equipped with the induced orientation. For τ ∈ S p+1 (K) and σ ∈ S p (K), the incidence number d τ σ is defined as In the case of simplices, we denote by x 0 . . . x p the oriented p-simplex whose vertices are x 0 , . . . , x p and whose orientation is implied by this order of vertices. We assume the reader is familiar with exterior algebra and differential forms (see e.g. [31, I-III] and [1]). Let ω, α denote the action of p-covector ω on p-vector α. Differential p-form in a complex K [31, p. 226] is a set of smooth p-forms ω σ in the cells σ of K satisfying the following patch condition: if τ is a face of σ , then the trace ω σ | τ of ω σ equals ω τ in τ . In other words, ω σ (x), α = ω τ (x), α for all x ∈ τ and all p-vectors α in the plane of τ . This enables us to consider the set of pforms ω σ as single p-form ω such that ω(x), α is well-defined for those p-vectors α that are in the plane of the cell σ for which x ∈ σ − ∂σ . Hence, differential pforms in K can be integrated over p-cells. Denote by F p (K) the space of differential p-forms in K. Note that since the exterior derivative d commutes with trace, we have dω ∈ F p+1 (K) if ω ∈ F p (K), but the Hodge star ω is not necessarily in F n−p (K).
Formal sums σ i ∈S p (K) a i σ i of oriented p-cells with real coefficients are called p-chains of K [31, App. II, §6]. These form a vector space C p (K) for which the p-cells σ i constitute a natural basis (here σ i = 1σ i , the sum in which a j = δ ij , the Kronecker delta). The elements of the dual space C * p (K) are p-cochains of K. Following [31], we use σ i to denote also the cochain whose value is δ ij at the chain σ j . Then the p-cells σ i constitute the dual basis for C * p (K), and also cochains can be written as formal sums of cells. Negative coefficients indicate change of orientation so that −σ is the cell σ with opposite orientation. For computer implementations, chains and cochains can be considered as vectors consisting of the coefficients a i after a numbering has been chosen for the cells of K.
The boundary ∂τ of a (p + 1)-cell τ ∈ S p+1 (K) is the p-chain σ ∈S p (K) d τ σ σ . This defines the boundary map ∂ : C p+1 (K) → C p (K) for all chains by requiring it be linear. The coboundary map d : C * p (K) → C * p+1 (K) for cochains is defined by dX(c) = X(∂c). We use the same notation d as for the exterior derivative of forms. When cochains are considered as vectors, we can denote by d also the matrix with components d ij = d τ i σ j for τ i ∈ S p+1 (K) and σ j ∈ S p (K), since this is the matrix of the coboundary map. This notation is not explicitly needed in this paper, but we mention it so that dX makes sense also when cochain X is considered as a vector.
Since p-forms can be integrated over p-cells, each p-form ω yields a p-cochain whose values on chains are determined by integration of ω. Namely, the de Rham map C : F p (K) → C * p (K) is the linear map defined by where the second equality is the definition of integration on p-chains. The cochain Cω can be considered as an approximation of ω. In vector presentation, its components are the integrals of ω over the p-cells of K. We remark that Stokes' theorem implies Cdω = dCω for p-forms ω.
We invoke the inner product of R n to define norms as follows. If ω is a p-covector, denote by |ω| the norm induced by the inner product of R n [31, I, §12]. If ω is a pform, then |ω| denotes the function whose value at x is |ω(x)|. Hence, we may define the L 2 norm of the p-form ω as the L 2 norm of the function |ω|. This is denoted by |ω| L 2 ( ) . In other words, |ω| Discrete exterior calculus [13,14,16] enables a discrete presentation of boundary value problems expressed in terms of differential forms. When differential forms are approximated with cochains, the coboundary operator naturally replaces the exterior derivative. Indeed, if Cω approximates ω, then Cdω = dCω suggests that dCω is the right approximation for dω. One also has to express the Hodge star operator for cochains. There are several ways to do this, and one typically employs a dual complex so that p-cochains of K are mapped to (n−p)-cochains of the dual complex. We need not consider any specific approach to deal with the Hodge operator or any specific boundary value problem. The framework we present in this paper can be applied as long as the solution is a cochain on a simplicial mesh which is meant to approximate a differential form -that is, its coefficients correspond to the integrals of the form.

Lowest order Whitney forms
In this section, we briefly recall Whitney forms as a tool for interpolating cochains. More information can be found in [20]. Let us henceforth assume that the mesh K is simplicial.
where dλ i indicates a term omitted from the product. For each p, the Whitney map W : is the space of Whitney p-forms and denoted by W p (K).
First, recall that CWX = X for all X ∈ C * p (K). Cochains and Whitney forms are in one-to-one correspondence in the simplest possible way, the cochain X = we can interpolate X with WX, and the integrals of this interpolant match with the values of X on p-simplices of K. Further, we can approximate p-form ω ∈ F p (K) with the Whitney form WCω, and the integrals of this approximation match with those of ω on p-simplices of K: CWCω = Cω. This approximation is exact for elements of W p (K), including constant p-forms.
For computing derivatives, a useful fact is that d and W commute: for X ∈ C * p (K), we have dWX = WdX. Finally, we mention the affine invariance property of Whitney forms. Let σ = x 0 . . . x n and τ = y 0 . . . y n be two n-simplices and ϕ : σ → τ the affine map such that ϕ( ). Because of this property, computations done in a reference simplex transfer to all simplices by affine transformations and hence need be done only once. This will be useful when we consider interpolation with higher order Whitney forms.

Small simplices and higher order Whitney forms
Lowest order Whitney forms defined in the previous section include constants and are at most first-order polynomials in each simplex. There are also higher order Whitney forms, or Whitney forms of order k, which include (k − 1)th-order polynomials and are at most kth-order polynomials in each simplex. Higher order Whitney forms can be defined in different ways. We use the so-called small simplices of [28], since this approach enables us to interpolate cochains with higher order Whitney forms.
To define the small simplices, let I(n + 1, k) denote the set of multi-indices with n + 1 components that sum to k; that is, For a fixed n-simplex σ = x 0 . . . x n , each multi-index k ∈ I(n + 1, k − 1) defines a map, which we denote by k σ , from σ to itself such that the point x whose barycentric coordinates are λ i maps to the point whose barycentric coordinates are λ i +k i k . In other words, k σ is defined by For k ≥ 1, the set of kth-order small p-simplices of K is Small simplices are homothetic images of the simplices of K. See Fig. 1 for examples of small simplices.
To each kth-order small p-simplex υ corresponds a kth-order Whitney p-form w(υ) ∈ F p (K), as given in the following definition. Henceforth, when σ = x 0 . . . x n is a fixed n-simplex, we denote by λ k σ the function n i=0 (λ i ) k i . The space of kth-order Whitney p-forms is the span of all such forms: Since the higher order Whitney forms in the definition above are products of barycentric functions and lowest order Whitney forms, it is immediate that they too are affine invariant and elements of F p (K). However, it is equally evident that there is no similar correspondence between Whitney forms and cochains of C * p (K) in the higher order case. Namely, W p (K) is already isomorphic to C * p (K), and increasing the order increases the dimension of the space, so there is no way for W p k (K) to be isomorphic to C * p (K). However, it is only natural that higher order approximations require more degrees of freedom, and the spanning higher order forms in Definition 4.1 were given corresponding to small simplices of K. This suggests that instead of C * p (K) we should concentrate on cochains over the small simplices of K -the next section makes this idea precise.

Interpolating with higher order Whitney forms
As can be seen from Fig. 1, the small simplices do not pave , and hence, they do not form a subdivision of K. However, we can always refine the mesh K such that the refinement contains the small simplices as cells. In this section, we consider how to interpolate cochains of this refined mesh with higher order Whitney forms.
Choose k and let K k denote a refinement of K into kth-order small simplices; more precisely, K k can be any mesh in that contains the small simplices of order k as cells. K k is allowed to have also other cells. For instance, between the small simplices of a tetrahedron in three dimensions, there are holes that are either octahedra or inverted tetrahedra (see Fig. 1). The refinement K k is not unique; the holes can be accepted as cells as such, or one may further divide them into simplices. The only requirements are that K k is a mesh in (i.e. satisfies the definition of mesh given in Section 2) and contains the small simplices of order k as cells. We can then consider chains C p (K k ) and cochains C * p (K k ) of K k and the de Rham map C k : F p (K) → C * p (K k ) (which is well-defined in F p (K k ), but we restrict the domain to the subset F p (K)).
To interpolate with higher order Whitney forms, we are looking for some kind of interpolating map V : C * p (K k ) → F p (K) akin to the Whitney map W. Since the coefficients a i of the cochain X = σ i ∈S p (K k ) a i σ i can be considered as integrals of some differential form that X is meant to approximate, we would like to find VX ∈ W p k (K) such that σ i VX = a i ∀σ i ∈ S p (K k ). In other words, we would like to have C k VX = X for all X ∈ C * p (K k ), preserving the property of W. However, it is evident that this is generally not possible, since the number of cells in K k is greater than the dimension of W p k (K). We will therefore have to relax the condition σ i VX = a i ∀σ i ∈ S p (K k ) somehow.
Since the spanning higher order forms in Definition 4.1 were given corresponding to small simplices of K, the first relaxation that comes to mind is to require σ i VX = a i not for all σ i ∈ S p (K k ) but for all σ i ∈ S p k (K) -that is, for those cells that are small simplices. However, the spanning forms in Definition 4.1 are not linearly independent, so the dimension of W p k (K) is lower than the number of small simplices in S p (K k ). Hence, we must relax the condition even further.
A simple and feasible way to deal with this linear dependency is to choose a subset S p k (K) of S p k (K) such the forms corresponding to small simplices inŜ ). This is how our earlier requirement has been relaxed. To summarise, VX has the correct integrals on the chosen subset of small simplices.
To make the above precise, we have to specify some details -namely, how to choose the subsetŜ p k (K) such that the corresponding Whitney forms constitute a basis for W p k (K). In the following we give the general idea in n-dimensions; details specific to three dimensions are considered in the next section.
Let σ n denote a generic n-simplex considered as a cell complex, and letW p k (σ n ) denote the subspace of p-forms in W p k (σ n ) with zero trace on the boundary of σ n . We rely on the decomposition [2, Theorem 7.3] where we have extended elements in W p k (σ q ) to elements of W p k (σ n ) using the barycentric extension (see [2, (7.1)]; simply consider all barycentric functions in σ n instead of σ q ), which will henceforth be applied implicitly when appropriate. This implies the global decomposition Hence, it suffices to choose the subsetŜ p k (σ q ) in a generic q-simplex σ q , considering only those small simplices that are not contained in the boundary of σ q . The same choice can then be applied throughout the mesh, yielding the required map V : C * p (K k ) → F p (K). Let us therefore consider q-simplex σ q and letS p k (σ q ) denote the set of those small simplices in S p k (σ q ) that are not contained in the boundary of σ q . The dimension ofW p k (σ q ) compares with the cardinality ofS p k (σ q ) as follows: 3) The first line has been proved in [1], and the second line is proved as Lemma A.3 in Appendix A. These counts determine how many small simplices one must omit before their cardinality matches the dimension of the Whitney forms in each simplex. The omitted simplices can be chosen in many ways, as long as the remaining forms span the same space. For this, let us recall the linear relations between higher order Whitney forms. When σ is a p-face of a (p + 1)-simplex τ , let τ − σ denote the node of τ that is opposite to σ . Then, for each τ ∈ S p+1 (σ q ) we have (see e.g. [19] for a proof) This shows that second-order p-forms are linearly dependent, and multiplying both sides by λ k σ q with k ∈ I(q + 1, k − 2) gives relations for kth-order p-forms. These relations can be used to ensure that we only omit redundant small simplices and the remaining p-forms still span the same space.
When p = 0 or p = q, (5.3) says dim(W p k (σ q )) = #S p k (σ q ). Thus none of the small simplices will be omitted in that case. When 0 < p < q, we multiply (5.4) applied to each (This requirement on k is set to obtain relations specifically for W p k (σ q ) instead of W p k (σ q ).) By Lemma A.2 in Appendix A, the number of such k is p + k − 1 q , and hence, we get q + 1 p + 2 Notice that the binomial coefficient p + k − 1 q involved in (5.3) appears here too; these formulas have a nice geometric interpretation which will be helpful in implementations. All of this will be clarified in practice in the next section. We conclude this section with the approximation property of higher order Whitney forms [20] and some remarks.
be the linear map obtained with a choice of kth-order small simplices as explained above, and let ω be a smooth p-form in . There exists a constant C ω,k such that whenever h > 0, C > 0, and K is a simplicial mesh in such that diam(σ ) ≤ h and (σ ) ≥ C for all simplices σ of K.

Remark 5.2
To vary the order of Whitney forms in different subdomains of , one has to divide adjacent elements into small simplices of different orders. This results in a division that is not a mesh in the sense of discrete exterior calculus. The interfaces of such subdomains require special treatment and changes in specific methods.

Remark 5.3 d and V
commute for 0-cochains, i.e. dVX = VdX for X ∈ C * 0 (K k ). For p > 0, this property cannot be achieved without giving up on other properties. In general dVX is given by VdC k VX, so the derivative can still be computed after one first corrects the values of X on omitted small simplices.

Systematic implementation strategy
Let us next turn our attention to the implementation of the ideas presented in the previous section. For simplicity and due to practical interests, we will henceforth restrict to the three-dimensional case n = 3. Although more work would be required in the case n > 3, it should be visible how the same algorithms generalise to any dimension. In contrast, we will not restrict to any specific order k. Our algorithms enable the implementation of kth-order 0-, 1-, 2-, and 3-forms in R 3 such that the order k can be given as a parameter and the same code covers all orders.
Let us first consider kth-order small simplices and the choice of the subsetŝ S p k (σ q ) in three dimensions. As explained, we will work in a generic q-simplex σ q and consider only those small simplices that are not contained in the boundary of σ q . Using (5.3) with 0 ≤ p ≤ q ≤ 3, we find that there are three types of small simplices that require attention: This means we have to omit k 2 small 1-simplices in σ 2 .
. This means we have to omit 3 k 3 small 1-simplices in σ 3 .
Our strategy to omit redundant small simplices was inspired by the example with k = 3 given in [11, p. 31]. We make use of the following geometrical observation, illustrated in Fig. 2 (recall also Fig. 1). In σ 2 , the holes that are not small simplices are inverted triangles and in correspondence to elements of I(3, k − 2). Their cardinality is k 2 , which is also the number of relations we get forW 1 k (σ 2 ). These relations correspond to the inverted triangles, and we omit one small 1-simplex for each. In σ 3 , there are two kind of holes: k 3 inverted tetrahedra corresponding to elements of I(4, k − 3) and k+1 3 octahedra corresponding to elements of I(4, k − 2). We have 4 k 3 relations forW 1 k (σ 3 ); 4 per inverted tetrahedron (one for each of its 2-faces). These relations allow us to omit any 3 of the edges of each inverted tetrahedron to obtain a basis forW 1 k (σ 3 ). ForW 2 k (σ 3 ), we have one relation per octahedron, and omitting one 2-face for each octahedron yields a basis forW 2 k (σ 3 ). As stated, it suffices to make these choices only once in a generic q-simplex, and the same choice can then be applied throughout the mesh. For other small simplices than the three types covered above, none of the small simplices are omitted. However, we remark that in practice one should exclude duplicate entries from the setsS 0 k (σ q ). Namely, the same small 0-simplex may appear multiple times in (4.1). To handle this, we can choose to always label the small 0-simplices as the images of the first vertex of σ q .
Let us denote the chosen subsets ofS p k (σ q ) byŜ p k (σ q ) and choose a numbering for their elements. As explained in the previous section, these subsets determine a map V : C * p (K k ) → F p (K) which enables us to interpolate cochains of K k with kth-order Whitney forms. The implementation process can be divided into three steps: 1. Given a simplicial mesh K in , we form a refinement K k containing the kthorder small simplices as cells. 2. Given a cochain X ∈ C * p (K k ), we solve the coefficients of the interpolant VX in the chosen basis. 3. Given the coefficients of the interpolant VX, we show how to evaluate it at a given point x ∈ .
These steps are discussed in detail in the following three subsections.

Refining the mesh
Suppose we have a three-dimensional simplicial mesh K in . In this subsection, we discuss how the mesh can be refined systematically to obtain a refinement K k that contains the kth-order small simplices as cells. We make the following assumptions for the mesh as data structure: • The p-simplices of K have been indexed, and the positions of the 0-simplices (which are points in ) are contained in a list. • The vertices of each p-simplex σ ∈ S p (K) can be accessed in a definite order, and this order determines the orientation of σ . • Given a list of vertices x 0 , . . . , x p , we have means of finding the p-simplex with these vertices (if one exists).
These can be achieved when, in addition to the position list, we store the indices of the (p − 1)-faces and the parent (p + 1)-simplices of each p-simplex σ . When refining the mesh, one might wish to store the indices of the small simplices that correspond to the elements of the setsŜ p k (σ q ) in each q-simplex σ q ∈ S q (K). They will be needed when interpolating cochains, and having saved the indices in memory, one need not repeatedly find them for each interpolation. While possible, this is not necessary if new indices are allocated in increasing order, as the indices of the small simplices will then be implied from those of the big simplices. In this case finding them is a very small task, while storing the indices can use a lot of memory. Nevertheless, if desired, the index of each small simplex can be stored right after it is added to K k , and hence, we need not make this explicit in the rest of this subsection.
The strategy for refining the mesh is as follows. We start with an empty mesh K k and make it the desired refinement of K by going through the simplices of K in the order of increasing dimension and adding the corresponding small simplices into K k . With this order, we only have to consider the small simplices that are not in the boundary of the simplices, since those in the boundary have already been covered along with some lower-dimensional simplex. For 2-and 3-simplices, we will also fill the holes that would otherwise be left between the small simplices. Multi-indices are needed to label both small simplices and the holes, and hence, we assume access to the sets I(l, m) stored in some data structure.
The first step is straightforward: we copy the 0-simplices of K into K k . Second, the small 0-and 1-simplices of each 1-simplex of K are added into K k . Third, we go through the 2-simplices of K and add the corresponding small 0-, 1-, and 2-simplices into K, also filling the inverted triangles in between. Fourth, we add the small 0-, 1-, 2-, and 3-simplices and the holes of each 3-simplex of K into K k . The holes that are octahedra are the only cells in K k that are not simplices; however, it is possible to divide them into four tetrahedra. In discrete exterior calculus, it is sometimes desirable that cells are well-centered. In this case one might wish to divide each octahedron into four tetrahedra by adding a 1-simplex of smallest possible length. If this is not uniquely determined, the octahedron can accordingly be divided into two pyramids or kept as it is. The potential division has no effect on the algorithms of this paper.
When adding the small simplices of σ q = x 0 . . . x q ∈ S q (K) into K k , we first add the small 0-simplices that are in the interior of σ q . These are obtained as images of x 0 through k σ with k ∈ I(q, k − 1) such that k i = 0 for i > 0. Each small 0-simplex is covered exactly once, so we do not have to check if another 0-simplex has already been added in the same position. For adding other small simplices (and the holes), we need means of finding the indices of the small 0-simplices corresponding to their vertices (which have already been added into K k ). It is useful to have a function that returns the index of the small 0-simplex k σ (x i ) for any k, σ = x 0 . . . x q , and i ∈ {0, . . . , q}. The indices of the 0-simplices can then be used to add higher-dimensional cells (or find if they exist); these are added in the order of increasing dimension.
The approach is simple on paper and also relatively easy to implement. Algorithm 1 summarises our strategy for refining the mesh.
number of possibilities we have to precompute, and then the matrix B(p, σ q , σ r ) will always be one of these. For example, consider B(1, σ 3 , σ 2 ). There are four 2faces of σ 3 , and there are six permutations for the vertices for σ 2 . Hence, there are 24 possibilities for the matrix B(p, σ q , σ r ). In 3D, the number of possibilities never exceeds 24 (see Table 1).
Having precomputed these matrices B(p, σ q , σ r ), we do not have to integrate anything when solving the coefficients. When σ r is a face of σ q and we wish to take the Whitney forms corresponding to small simplices of σ r into account, we use the matrix B(p, σ q , σ r ) (the one appropriate for σ q and σ r ) and the coefficients we have solved in σ r earlier. Matrix-vector multiplication yields the integrals over elements ofŜ p k (σ q ). The integrals of higher order Whitney forms over small simplices that are needed for the matrices A(p, q) and B(p, σ q , σ r ) can be computed analytically; we elaborate on this in Appendix B. Another option is to use at least kth-order quadrature formulas for numerical integration. The matrices are formed only once, and hence, we can precompute them in as high numerical precision as desired.
Algorithm 2 summarises our strategy for solving the coefficients of the interpolant.

Evaluating the interpolant at a given point
Suppose we have solved the coefficients c[υ i ] of the interpolant in the chosen basis for W p k (K) and we wish to evaluate the interpolant at a given point x ∈ . In this subsection, we discuss how to compute the value of υ i ∈Ŝ p k (K) c[υ i ]w(υ i ) at x. Note that this is not an entirely trivial task, for we have used the decomposition (5.2) when forming a basis and hence have to accumulate contributions ofW p k (σ q ) for different σ q ∈ S q (K), p ≤ q ≤ 3.
We assume it is possible to search for a 3-simplex σ 3 ∈ S 3 (K) that contains x (or one is known in advance). Notice that w(υ) is zero in σ 3 if υ is not contained in σ 3 , and hence, it suffices to consider the basis forms corresponding to small simplices of σ 3 . If x happens to be in multiple 3-simplices, their intersection is some q-simplex σ q ∈ S q (K) for q < 3. In this case we may either compute the value in σ q (considering only the small simplices of σ q ) or choose one of them and compute the value there; the trace on σ q agrees with the value computed in σ q .
Supposing we have found a 3-simplex σ 3 containing x, the next step is to compute the barycentric coordinates of x and the values of the lowest order Whitney forms corresponding to p-faces of σ 3 . These can then be used with the coefficients c[υ i ] to compute the value of the interpolant. Indeed, since the basis forms are products of barycentric functions and lowest order Whitney forms, we can write the interpolant as where the second sum is over p-faces σ p i of σ 3 and d i is a linear combination of products of barycentric functions. To compute d i , recall the decomposition (5.1) and accumulate the contributions from all faces of σ 3 using the coefficients c[υ i ]. Taking orientations into account, our strategy for evaluating the interpolant at a given point is formulated in Algorithm 3.

Numerical examples
To show that our algorithms can be implemented in practice, we provide numerical examples with higher order Whitney forms. In all of the test cases, our domain is the rhombic dodecahedron with vertices (±1, ±1, ±1), (±2, 0, 0), (0, ±2, 0), and (0, 0, ±2). We build a simplicial mesh K in and its refinement K k . For a p-form ω, the cochain C k ω is obtained by integrating ω using high-order quadrature formulas for p-simplices. Then we can approximate ω with the kth-order Whitney form VC k ω. We compute the L 2 norm (again using quadrature formulas) of the error |VC k ω − ω|. The experiments are performed for k ∈ {1, . . . , 12} using several test functions ω. All test functions have L 2 norm close to one so that the use of absolute error is appropriate. To study convergence properly before running out of machine accuracy, computations are done in quadruple precision.
First, we confirm that our algorithms work as expected using polynomial test functions on a coarse mesh. The mesh has 24 tetrahedra and the maximum edge length is 2.0. Our test functions ω p,j are (where the label j ∈ {1, 2, 3}) The results are displayed in Table 2. As expected, the approximation becomes exact as soon as the test function is in the space of kth-order Whitney forms.
We also study the convergence of the approximations with respect to the maximum edge length h of the initial mesh K. Our test function is the 1-form ω(x, y, z) = 1 4 sin(2y) cos(2z)e x 2 /4 dx + sin(2z) cos(2x)e y 2 /4 dy + sin(2x) cos(2y)e z 2 /4 dz .   We approximate ω using four meshes that have 24,192,1536, and 12288 tetrahedra with maximum edge lengths 2.0, 1.0, 0.5, and 0.25 respectively. The results are displayed in Table 3 and illustrated in Fig. 3 (in log-log scale). We conclude that higher Fig. 3 Illustration of the results displayed in Table 3 = With Lemmas B.1 and B.2 we can compute averages of barycentric products over small simplices. The integrals of higher order Whitney forms over small simplices are then easily obtained using the following result.

Proposition B.3
Let σ q be a q-simplex, τ ∈ S p (σ q ) a p-face of σ q , and ω a smooth 0-form. For any p-simplex υ ⊂ σ q , we have where 1 |υ| υ ω is the average of ω over υ, x is any point in υ, and vect(υ) is the p-vector of υ.