Construction and application of algebraic dual polynomial representations for finite element methods on quadrilateral and hexahedral meshes

Given a sequence of finite element spaces which form a de Rham sequence


Introduction
The differential operators grad, curl and div play a fundamental role in the representation of physical field laws.In terms of differential geometry these operators are either represented by a topological, metric-free operation called the exterior derivative or by the metric-dependent operation called the codifferential or coderivative denoted by δ or d ⋆ , see [1].The codifferential is the Hilbert adjoint of the exterior derivative and the metric dependence of the codifferential is a direct consequence of the definition of inner product, which involves the metric.Although the distinction between exterior derivative and codifferential is most clearly seen in terms of differential geometry, it can also be introduced in terms of vector calculus as weak differential operators, see for instance [2,Appendix A].
The construction of the exterior derivative in a discrete setting is relatively straightforward since it follows from the continuous definition.The challenge lies in the construction of the sequence of discrete function spaces such that, together with the exterior derivative, they constitute an exact sequence (a de Rham sequence), see for example [3].
The subject of constructing discrete representations for the codifferential operator is a more challenging one.This operator contains all metric and material properties and its discretization typically establishes the distinguishing properties of each numerical method.Different aspects of the codifferential operator can be taken into account when constructing its discrete counterpart.Apart from its metric-dependence, the discrete codifferential is not just the restriction of the continuous operator to an appropriate finite dimensional subspace.In weak formulations this generally corresponds to using integration by parts to convert the codifferential into an exterior derivative.This approach stems directly from the fact that finite element formulations typically use one de Rham sequence (a primal sequence of spaces).By employing integration by parts, the discrete codifferential can be constructed using only one single sequence of spaces, see for example [3].
In the finite difference/finite volume community the definition of the codifferential operator follows a different route, relying on the construction of a dual grid, see [4].For the evaluation of the codifferential, the scalar or vector field is converted to the dual grid, there the exterior derivative is computed after which the result is mapped back to the primal grid.The mappings which convert the variables on the primal grid to the variables on the dual grid are discrete representations of the Hodge-⋆ operator and therefore are called Hodge matrices.This approach, instead of defining the discrete codifferential as the dual of the exterior derivative, directly approximates the explicit representation of the codifferential where d is the exterir derivative, d ⋆ the codifferential operator, d denotes the dimension of the ambient space, k refers to the k-form for which the derivative is taken, and ⋆ denotes the Hodge-⋆ operator see [1,5,6].
In [7] both approaches are compared for a spectral element method.On a single grid, integration by parts is used to convert the codifferential to an exterior derivative, while in the dual grid approach dual spectral element basis functions are constructed on a dual grid.
As mentioned before, the construction of discrete codifferential operators is challenging, especially for high order spectral element methods.The underlying reason has to do with the fact that the codifferential is non-local.The higher the approximating degree of the method, the denser the matrix representation of the codifferential and the higher the condition number of the resulting system of equations becomes.These two properties are particularly impactful to any numerical method.This is the focus of this work: to propose a high order discrete codifferential operator with improved sparsity and condition number.
In this paper we use integration by parts and interpret the mass or Gram matrices as Hodge matrices.These Hodge matrices correspond to discrete Hodge-⋆ operators that convert the primal representation to a dual representation.By doing so, the integration by parts formula becomes metric-free and the codifferential is interpreted as the derivative of dual variables.This leads to much sparser system matrices-especially for higher order methods.We set up a de Rham sequence for primal spaces and for each space in this complex we construct a dual space.The sequence of these dual spaces also forms a de Rham complex.The construction of the dual de Rham sequence explicitly requires that the de Rham sequence should also be satisfied on the boundary.The construction of a dual basis used in this work is similar to the inverse Gram constructions described in [8,9].
Earlier dual bases have been used in isogeometric methods for projection of B-splines [10,11].In finite elements they have been used for mortar methods [8,12].In [13] these ideas have been combined to form isogeometric mortar methods.For other implementations in isogeometric methods see [14,15].Different approaches for the construction of dual bases have also been discussed in [9], and for construction of dual splines in [16].
As applications, we will demonstrate the use of a dual basis on a constrained minimization problem of the Poisson equation in 3D.We also show the discrete well-posedness of this problem which is expressed in terms of degrees of freedom only.It will be shown that the use of an algebraic dual basis results in a very sparse matrix where two of the sub-matrices consist of 1, −1 and 0 only and do not change with the shape and size of the elements of the mesh as long as the mesh topology remains the same.This observation is relevant for the incompressible flow equation where we encounter a similar div-grad pair.These techniques may also be valuable in electromagnetism to represent the involution constraint div B = 0 in a way that is very sparse and independent of the shape and size of the mesh.
In the second application we will solve for the pair of Dirichlet-Neumann problems discussed in [17] and prove their equivalence in the discrete sense.This is in general not trivial.It is shown that the duality relation continues to hold point-wise in these finite-dimensional approximations, on arbitrary grids, through the use of algebraic dual polynomials.
In the third application, we will solve for a grad-div eigenvalue problem on affine and non-affine meshes using dual degrees of freedom and will show optimal convergence rates for the eigenvalues.
The remainder of this paper is structured as follows: The construction of dual polynomial spaces in the one dimensional case is presented in Section 2. In this section it is also shown how nodal sampling and edge sampling from polynomial spaces extend to Sobolev spaces.The treatment of the multiple elements case and the derivative of a dual representation will be given in this section.In Sections 3 and 4 these constructions are extended to two dimensions and three dimensions, respectively.In Section 5 a dual polynomial representation is used for the mixed formulation of the Poisson equation in the three-dimensional case with multiple elements.In Section 6 equivalence of the Dirichlet-Neumann problems [17] is proved and demonstrated by a particular example.In Section 7 we address a grad-div eigenvalue problem and show optimal convergence rates on affine and non-affine grids.Finally, in Section 8 conclusions are drawn and future work is discussed.

Construction of 1D dual finite elements
We will use the definition of finite element spaces in terms of the triplet (K , P, N ) by Ciarlet, [18], see also Ern and Guermond, [19, §1.2] and Brenner and Scott,[20,§3.1].Definition 1.A finite element consists of the triplet (K , P, N ) where i K is a compact, connected, Lipschitz subset of R d with non-empty interior; ii P is a (finite dimensional) linear vector space with domain K .Usually, P is a polynomial vector space of dimension d P ; iii N is a set of linear functionals {N i }, i = 1, . . ., d P , acting on elements of P, such that the linear map, The linear functionals {N i } are called the local degrees of freedom.The following Proposition taken from [19] defines the basis functions: Proposition 1.There exists a basis {Ψ 1 , . . ., Ψ dP } in P such that , where L N (ξ ) is the Legendre polynomial of degree N and L ′ N (ξ ) its derivative.These nodes are referred to as the Gauss-Lobatto-Legendre (GLL) points, [21].Let P be the space of polynomials of degree N defined on the interval K .For any p h ∈ P define the degrees of freedom by Because polynomials are continuous, (1) is well-defined.The superscript 0 in N 0 i indicates that we sample the polynomial p h in points.The basis which satisfies the Kronecker-delta property from Proposition 1 is given by the set of Lagrange polynomials through the GLL-points , i = 0, 1, . . ., N .
Remark 1.Note that the degrees of freedom are linear functionals on P. The nodal sampling of functions in P is essentially the Dirac delta distribution which is well defined when the space P consists of continuous functions, see [22,Example 2.10.2].Extension of this functional to Sobolev spaces in this way is in general not possible.The extension to Sobolev spaces will be given in Definition 4.
Example 2. Let K and ξ i be defined as in Example 1.Let Q be the space of polynomials of degree (N −1).For this example we choose to define the degrees of freedom as For polynomials, the integral in ( 2) is well-defined.The superscript 1 in N 1 i expresses the fact that the degrees of freedom are associated to line segments [ξ i−1 , ξ i ] (geometrical dimension d = 1).The basis functions, e j (ξ ), which must satisfy the Kronecker-delta property from Proposition 1 need to satisfy Lemma 1.The basis functions e j (ξ ) on the GLL-grid defined in Example 2 are given by e j (ξ ) = − where h k (ξ ) are the Lagrange polynomials through GLL-points defined in Example 1. Proof.

∫ ξ
where we repeatedly use the Kronecker-delta property of the Lagrange polynomials.If the Lagrange polynomials h k (ξ ) are polynomials of degree N, then dh k (ξ )/dξ is a polynomial of degree (N − 1).It is easy to show that {e i } N i=1 forms a basis for Q. □ Corollary 1. From (3) it follows that dh j dξ (ξ ) = e j (ξ ) − e j+1 (ξ ) .

So if p h
∈ P is expanded in terms of Lagrange polynomials as then its derivative is given by where we have used the fact that e 0 (ξ ) = e N+1 (ξ ) = 0. Let E 1,0 be the N × (N + 1) matrix then we can write (4) as Taking the derivative of a nodal expansion changes the nodal degrees of freedom discussed in Example 1 to the integral degrees of freedom of the derivative, as discussed in Example 2. The matrix E 1,0 is called the incidence matrix, which converts the nodal degrees of freedom of a function p to the integral degrees of freedom of its derivative dp dx .

Construction of dual basis
Consider the finite element constructed in Example 1.Any element p h ∈ P can be represented as where N 0 i (p h ) are the nodal degrees of freedom and h i (ξ ) are the associated basis functions.To simplify the notation, we will write this as where . . . .
Here, and in the remainder of the paper, basis functions will always be expressed as a row vector and the degrees of freedom as a column vector.Let p h , q h ∈ P be both represented as in (6), then the L 2 -inner product is given by where M (0) denotes the mass matrix (or Gram matrix) associated with the nodal basis functions Definition 2. Let N 0 (p h ) be the degrees of freedom for p h ∈ P, then the dual degrees of freedom, Ñ 1 (p h ) for p h ∈ P are defined by Remark 2. In Definition 2 the superscript 1 on Ñ 1 corresponds to lines (geometric dimension 1) that are geometric duals of nodes (geometric dimension 0) for the one-dimensional domain K = [− Proof.If the corollary is true, we have where we made use of the linearity of the operator N 0 .This shows that Ψ 1 = Ψ 0 ( M (0) ) −1 forms a basis for the dual degrees of freedom.□

Corollary 3. For any element p h
∈ P we have

Corollary 4.
The mass matrix M(1) is the inverse of the mass matrix M (0) .

Proof.
M(1) := where In Fig. 1 the Lagrange polynomials through the GLL-points and the associated dual polynomials are presented for N = 4.
Analogous to the construction of the dual nodal polynomials, we can also construct the dual polynomials to the edge functions.Let an element p h ∈ Q be represented as In the simplified notation this can be written as . . .We can write the L 2 -inner product for two functions p h , q h ∈ Q expanded in this way as with M (1) the mass matrix (or Gram matrix) associated with the edge polynomials Definition 3. Let N 1 (p h ) be the degrees of freedom for p h ∈ Q, then the associated dual degrees of freedom Ñ 0 (p h ) are defined as ) .
Here again we follow Remark 2 to denote the superscript 0 on Ñ 0 corresponding to nodes (geometric dimension 0) that are geometric dual of 1 for K = [−1, 1].In general, for the d-dimensional case the dual degrees of freedom of N 1 will be denoted by Ñ d−1 , see also Section 3.
Following Corollary 2, the dual edge functions are then given by In Fig. 2 the edge polynomials e i (ξ ) and their dual polynomials ẽi (ξ ) are shown for N = 4.
Lemma 2. Let Ψ k (ξ ) and Ψ 1−k (ξ ) for k = 0, 1 be the primal and the dual basis as defined above, then these bases are bi-orthogonal with respect to each other where I is the (N + 1) × (N + 1) identity matrix for k = 0 and the N × N identity matrix for k = 1.
Proof.Using the definition of the dual basis, where in the second term we have used the fact that mass matrix M (k) is symmetric.□ In Remark 1 it was stated that nodal sampling of a function is only possible in the space P of continuous functions.In a Sobolev space the elements consist of equivalence classes of functions that satisfy an integral equation and in this case nodal sampling may not be defined.

Lemma 3. Let p h
∈ P, then the nodal degrees of freedom are given by where in the last equality we used Lemma 2. □ Lemma 3 allows us to extend nodal sampling to square integrable functions.Definition 4. For f ∈ L 2 (K ) we define the nodal degrees of freedom by Using now the fact that Ψ 1 (ξ ) = Ψ 0 (ξ ) ( M (0) ) −1 this 'nodal sampling' can be written as which is just the L 2 -projection of f on the space P. Analogously we have These definitions are used, for e.g., in the calculation of degrees of freedom of the right hand side term and the boundary conditions in Sections 5 and 6.

Transformation rules for 1D function spaces
So far we have discussed the construction of polynomial basis on a reference domain K = [−1, 1].In case of a more general interval I = [a, b] we have to introduce mappings.For any x ∈ I, in case of a linear mapping Φ, we have and its Jacobian J = dx dξ = (b−a) 2 .In this case we have that the nodal basis and the edge basis transform as The associated mass matrices are then given by

Transformation rules for dual function spaces
The construction of dual function spaces for I = [a, b] follows exactly the same procedure as in Section 2.1.Let Ψ k (x), for k = 0, 1, be the transformed basis functions, and N k (p h ) for k = 0, 1, be the associated degrees of freedom, then we have that

Multi-element case
Let I = [a, b] with a, b ∈ R and a < b.Let the domain I be partitioned in K el elements bounded by the points In the multi-element case we require that polynomials q h , expanded by the Lagrange polynomials in each element, are in the space G(I) ⊂ H 1 (I), with G(I k ) = P.The global basis functions, Ψ 0 i (x) are related to the element basis functions by where k = 1, . . ., K el and j = 0, . . ., N.Here we use a simple linear transformation between the canonical element and the physical elements of the mesh, as presented in Section 2.2.On the standard domain [−1, 1] we will use the variable ξ , while in the global, multi-element case we will denote the independent variable by x.Just as in the single element case, the dual degrees of freedom are obtained by premultiplying the global degrees of freedom with the global (assembled) mass matrix.The corresponding dual basis functions are constructed by post-multiplying the row vector of global basis functions by the inverse of the global, assembled mass matrix.An example with I = [−1, 1], K el = 5 and N = 1 is shown in Fig. 3.
The global mass matrix which converts the primal degrees of freedom to the dual degrees of freedom used in the construction of Fig. 3 is given by , and its (approximate) inverse is given by where the underlined values are the nodal values plotted for the dual basis in Fig. 3. Remark 3. Explicit construction of dual basis functions requires the multiplication with the inverse of the global mass matrix.For large, multi-dimensional problems this is unfeasible.But the explicit construction of dual polynomials is in general not required to construct the system matrices of a given problem.It is sufficient to use the properties of the dual representation, such as the fact that it forms a bi-orthogonal basis to the primal basis, see Lemma 2. It is in the post processing step when we need to reconstruct dual representations that we solve a linear system of equations to convert dual degrees of freedom to primal degrees of freedom and then use primal basis functions for reconstruction, whenever necessary.
We will refer to the space spanned by the dual polynomials as G(I), which explicitly shows that the functions are expanded in terms of dual polynomials, despite the fact that G(I) = G(I).
In the multi-element case the representation in terms of edge functions is discontinuous between elements.Its dual representation, which is a linear combination of primal edge functions, will therefore also be discontinuous between elements.The global representation in terms of edge functions, which consists of discontinuous piecewise polynomials of degree N − 1 when a linear map is used, will therefore form elements of the subspace S (I) ⊂ L 2 (I) only, with S(I k ) = Q.The global basis functions, Ψ 1 i (x), are related to the local basis functions through where i = 1, . . ., N • K el , k = 1, . . ., K el and j = 1, . . ., N. In Fig. 4 a global discontinuous representation in terms of edge functions is shown together with its global dual representation for I = [−1, 1], K el = 5 and N − 1 = 0.

Differentiation of dual variables
Since the functions in G(I) are continuous, we can consider the restriction of these functions to the boundary.We will denote this space as G(∂I).For the one-dimensional case this space consists of the function's values in the points x = a and x = b.The dual representation of a function in a point is in this case equal to the primal representation in a point.

Definition 5. Consider the multi-element case and let q h
∈ G (I) be expanded in Lagrange polynomials and φ h in dual edge polynomials in each element.Since φ h / ∈ H 1 (I), we define its weak derivative: For given φh ∈ G (∂I) we define d dx : where, ( φh (a), φh (b)) ⊺ ∈ G(∂I).
To clarify the notation we will denote the degrees of freedom on the boundary by B 0 instead of N 0 .Therefore, we have Since d dx (φ h , φh ) ∈ G(I) we can expand it, using the global basis functions dual to (12) as where just as in (9) Ψ 1 (x) is a row vector of global basis functions and Ñ 1 ( d dx (φ h , φh ) ) is a column vector of global degrees of freedom.The function q h is expanded in terms of the global basis functions as With these global representations the integral on the left hand side in ( 14) can be written as ) , because the two dual bases Ψ 0 (x) and Ψ 1 (x) are bi-orthogonal by construction, see Lemma 2. Note that we can use this property without actually constructing the dual basis functions.So the integral on the left hand side of ( 14) reduces to the vector product of the degrees of freedom of the two representations.
The first integral on the right hand side of ( 14) evaluates to because the derivative of a piecewise Lagrange representation q is expressed in terms of piecewise edge polynomials, Corollary 1, and φ h is expanded in terms of dual edge polynomials.Both bases are bi-orthogonal therefore the explicit dependence on the basis functions cancels from this integral and the integral can be expressed in terms of degrees of freedom and a topological incidence matrix, only.Finally, the boundary terms reduce to Since ( 14) needs to hold for all q h ∈ G(I) -and therefore for all expansion coefficients N 0 (q h ) -this reduces to where the inclusion matrix N maps the degrees of freedom on the boundary B 0 (q h ) to the global degrees of freedom.For K el = 5 and N = 1, the incidence matrix is given by Note that the matrices E 1,0 and N are independent of the mesh size, h, and the shape of the elements and only depend on the topology of the mesh.The flux is considered positive from left to right.On the left side of the interval, the outward unit normal points to the left and therefore the entry in N is −1, while on the right boundary the flux and the outward normal point in the same direction and therefore the entry in the N matrix is +1.
This makes the derivative of the dual representations also a topological operation in contrast to the codifferential applied to the primal variables which depends explicitly on the metric.Note also that the matrices E 1,0 and N do not depend on the polynomial degree N, so even for high order methods the derivative of the dual variables is expressed in very sparse matrices, which only contain the non-zero entries −1 and 1, and the dual degrees of freedom.
As an example of the derivative of dual representations, we take a piecewise constant function which approximates − cos(2π x) and take its derivative as described above for K el = 15 in Fig. 5, and K el = 100 in Fig. 6.In this case the boundary values φh are set to the value of the function − cos(2π x) at the boundary of domain, B0 ( φh Remark 4.Although the Hodge star operator in differential geometry is well-defined, there is some freedom numerically as can be seen by the ever growing number of papers on discrete Hodge operators and Hodge matrices; primarily in the finite difference/finite volume literature.In the finite element community this is less common, but we can interpret the dual representation as a finite dimensional Hodge dual of the primal representation.In finite element formulations one generally employs inner-products for the weak formulation.An inner-product can also be written as a combination of a (topological) wedge product and a Hodge star operator using where ω and η are two differential k-forms, see, for instance, [3, §4, pg.320] or [6, pg. 18].The dual representation is the finite element analogue of the Hodge star operator applied to the differential form η. That the mass matrix, which converts primal to dual representations, acts as Hodge operator in a finite element setting was already observed by Tarhasaari et al. [23].In this paper the mass matrices M (0) and M (1) (and later the mass matrices along the boundary, M 0 b and M 1 b ) play the role of Hodge operator.By contracting the mass matrix into a new variable, called here the dual variable, the inner-product is converted to a metric-free wedge product.

Discrete de Rham sequence
Using the function spaces defined in Sections 2.3 and 2.4, we can write the primal de Rham sequence as , and the dual de Rham sequence as

Two-dimensional dual spaces
In order to address more challenging problems, it is important to consider in detail the finite element spaces in 2D.
For d = dim Ω = 2, we have two sets of function spaces that obey the de Rham cohomology [3,24] H (curl; Ω) and , where curl of a scalar field ψ is the vector field (∂ψ /∂y, −∂ψ/∂x) ⊺ , the div applied to a vector field u = (u, v) ⊺ is the scalar field ∂u/∂x+∂v/∂y, gradient of a scalar field φ is the vector field (∂φ/∂x, ∂φ/∂y) ⊺ and the rot for a vector field u = (u, v) ⊺ is the scalar field given by (∂v/∂x − ∂u/∂y).We will introduce the three finite element spaces, , and the corresponding dual spaces C (Ω), D (Ω), S (Ω), respectively.We will first present the construction of finite element spaces on the reference element K = [−1, 1] 2 , and then on an arbitrary element Ω k ⊂ R 2 .This will be followed by construction of finite element spaces in case of multiple elements, and then the derivation of differential operators on the dual representations.

The function space C
], i = 0, . . ., N, be the GLL-points, and P denote the space of polynomials of degree N defined on the interval [−1, 1], see Example 1.Consider now the polynomial tensor product space C (K ) := P ⊗ P. Given the set x of 2D nodes x k defined as x := {x i(N+1)+j = (ξ i , η j ) | i, j = 0, . . ., N}, we can introduce for any p h ∈ C (K ) the degrees of freedom as The basis which satisfies the Kronecker-delta property from Proposition 1 is given by the Lagrange (or nodal) , i, j = 0, . . ., N, where h i are the 1D nodal interpolants introduced in Example 1.A visual representation of these basis functions for N = 2 is presented in Fig. 7a.

The dual finite element
The construction of the dual basis functions follows the ideas presented in Section 2.1.Here we outline the direct application to the 2D case of constructing the dual basis of the space C (K ).
Definition 6.The degrees of freedom of the dual element are given by where ) .
Since the dual basis functions ε(2) j need to satisfy the Kronecker-delta property by Corollary 2 we have that the dual basis functions can be expressed in terms of the primal basis functions as A visualization of the dual basis functions in the reference domain for N = 2 is presented in Fig. 7b.

The function space D(K
. ., N, be the GLL-points, and P and Q denote the space of polynomials of degree N and N − 1, respectively, defined on the interval [−1, 1], as in Examples 1 and 2. Consider now the polynomial tensor product spaces D ξ := P ⊗ Q and D η := Q ⊗ P. We introduce for any polynomial vector field where e ξ , and e η are the unit vectors in the ξ -and η-directions, respectively.In a polynomial vector space these integrals are well-defined. It is possible to show, see [7,[25][26][27], that the basis functions which satisfy the Kronecker-delta property from Proposition 1 are the 2D edge polynomials, ϵ (1)  k , k = 1, . . ., 2N(N + 1), defined as where h i are the 1D nodal interpolants introduced in Example 1, and e j are the 1D edge interpolants introduced in Example 2. A visual representation of these basis functions for N = 2 is presented in Fig. 8a.
Let ω h ∈ C (K ) be represented as ) , using Corollary 1, we have where E 1,0 is the incidence matrix for the edge and nodal degrees of freedom of a 2D element. 1   If R(curl; C (K )) denotes the range of the curl operator applied to elements from C (K ), this implies that R (curl; C (K )) ⊂ D(K ).This is a necessary requirement for C (K ) and D(K ) to form part a finite dimensional de Rham sequence.

Dual finite element
The construction of the dual basis functions of the space D(K ) is done in the same way as for the dual basis functions of the space C(K ).In this case, the dual basis functions are expressed in terms of the primal basis functions as A visual representation of these basis functions for N = 2 is presented in Fig. 8b.
Definition 7. The dual degrees of freedom for p h ∈ D (K ) are given by . ., N, be the GLL-points, and Q represent the space of polynomials of degree N − 1, as in Example 2. Consider now the polynomial tensor product space S (K ) := Q ⊗ Q.The degrees of freedom for this finite element can be introduced for any polynomial p h ∈ S (K ) as These integrals are well-defined in a polynomial space.It is possible to demonstrate, see [7,[25][26][27], that the basis functions which satisfy the Kronecker-delta property from Proposition 1 are the surface polynomials, ϵ (2)  k , k = 1, . . ., N 2 , defined as where, as before, e j are the 1D edge interpolants introduced in Example 2. For ease of notation we will write these basis functions as ) .
A visual representation of these basis functions for N = 2 is presented in Fig. 9a.
An element q h ∈ D(K ) can be represented as a linear combination of the basis (17) as ) . ( If we take the divergence of this vector field and use (4) repeatedly, we have So, we see that the divergence modifies the degrees of freedom (the expansion coefficients) and changes the basis functions from basis functions in D(K ) to basis functions for S(K ).We can write this as 1 The incidence matrix E 1,0 is not the same as used in (5) which is for degrees of freedom of a 1D element.where the incidence matrix E 2,1 is a sparse matrix which only contains the non-zero entries, −1 and 1, that can be obtained from (20).
Application of (4) shows that R(div; D(K )) ⊆ S(K ), which is required for the spaces D(K ) and S(K ) to be part of the finite dimensional de Rham sequence.
Since we have that div curl ω h ≡ 0 for all ω h ∈ C (K ), we have, by combining ( 18) and ( 21) that ≡ 0, which is a well-known property of incidence matrices in mimetic methods.

Dual finite element
The construction of dual basis functions of the space S(K ) follows the same steps as performed for the spaces C(K ) and D(K ).Definition 8.The dual degrees of freedom for p h ∈ S (K ) are given by The associated dual basis functions are expressed in terms of the primal basis functions as A visual representation of these basis functions for N = 2 is presented in Fig. 9b.

Transformation rules for 2D function spaces
In the previous sections we have introduced the construction of function spaces on a reference domain K = [−1, 1] 2 .In this section we will extend this construction to an arbitrary Ω k ⊂ R 2 .Let Φ k be the diffeomorphism between the canonical domain K and the arbitrary domain Ω k , and J be its Jacobian tensor such that In this work, we construct the diffeomorphism using transfinite mapping [28].The mapping is performed element-wise such that the global map is continuous.

Transformation rules for spaces
k is the pullback operator.We can then reverse the relation to obtain

Transformation rules for spaces D (Ω k )
Given a vector field ū ∈ D(K ), the transformed vector field u ∈ D(Ω k ) on Ω k is given by It is possible to compute the inverse of this transformation, resulting in The inverse relation can be computed in a similar fashion yielding

Transformation rules for dual function spaces
The construction of dual spaces for an arbitrary element Ω ⊂ R 2 follows the same procedure as that for the canonical domain K .
Let Ψ k (x) for k = 0, 1, 2 be the transformed basis functions in space C (Ω), D (Ω) and S (Ω), respectively, and let the associated mass matrix be given by Then we have the transformed dual basis functions and the dual degrees of freedom given by

Multi-element case
Let a two-dimensional, bounded domain Ω with Lipschitz continuous boundary be partitioned in K el non-overlapping conforming elements.In the multi-element case, the global dual representation is constructed similarly to the onedimensional case described in Section 2. For the global piecewise polynomials in D(Ω) ⊂ H(div; Ω), the degrees of freedom at the boundary, which represent integrated fluxes, are shared by adjacent elements.The normal component of these vector fields is continuous between adjacent elements, while the tangential component is discontinuous.The dual representation has different degrees of freedom and different basis functions, but being a linear combination of basis functions in D(Ω), these functions also have continuous normal component and discontinuous tangential components between neighbouring elements.
Finally, the piecewise polynomials in S(Ω) ⊂ L 2 (Ω) are discontinuous between elements and therefore so are the elements of S(Ω).

Vector operations on dual variables
Just as in the one-dimensional case we will define weak derivatives for which we use the dual representation.
Note that the standard gradient cannot be applied to elements of S(Ω) because these elements lack the required smoothness.We see that in (23) the integral on the left hand side is a bilinear form on D(Ω) × D(Ω), the first integral on the right hand side is a bilinear form on S(Ω) × S(Ω) and the boundary integral on the right is a bilinear form on D(∂Ω) × D(∂Ω).
Just as in the one-dimensional case we use standard finite element assembly of local basis functions and degrees of freedom to form global basis functions and degrees of freedom, which we will denote by Ψ 1 (x) and N 1 , respectively.An element q h ∈ D(Ω) can be expanded as Since grad maps into D(Ω) we can expand grad (s h , ŝh ) as This shows that the integral on the left hand side of ( 23) evaluates to If the expansion of q h is given by (24), div q h is expanded as div q h (x) = Ψ 2 (x)E 2,1 N 1 (q h ) , (26) and The first integral on the right hand side of (23) gives where once again, due to the bi-orthogonality of the primal and dual bases, the integral can be expressed in terms of the degrees of freedom only (no integration points or integration weights are required).
The restriction of q h to the boundary ∂Ω, is given by In this case N 1 is the 2D inclusion matrix which maps degrees of freedom on the boundary, D(∂Ω), to the global degrees of freedom in D(Ω), where the subscript 1 on the matrix denotes that degrees of freedom are associated to edges (geometric dimension 1).Just as its 1D analogue, see for instance (16), this sparse matrix only has −1 and 1 as non-zero entries.See Example 3 on how to construct N 1 .The basis functions Ψ 1 b (x) = Ψ 1 (x)N 1 form a basis for the trace space D (∂Ω) and the degrees of freedom N ⊺ 1 N 1 (q h ) the associated degrees of freedom.The mass matrix of this trace space is given by Using this mass matrix, we can set up the dual basis functions in the trace space given by Ψ 0 . With these dual basis functions, we can express ŝh ∈ D(∂Ω) as where, using Definition 4 we have With the expansions ( 28) and ( 30) we can write the boundary integral in (23) as Collecting all integrals in (23) we have This equation needs to be satisfied for all q h , i.e.N 1 (q h ), so we need to have Example 3. In Fig. 10 a small mesh is shown and the degrees of freedom N 1 i for elements in D(Ω) and the boundary degrees of freedom B 1 j for functions in D(∂Ω).The matrix N 1 , which plays a role in (34), in this case is given by .
N 1 is a dim(D(Ω)) × dim(D(∂Ω)) matrix.On the eastern and northern boundary the outward unit normal points to the positive x-and y-direction, respectively, and therefore the boundary degrees of freedom B 1 j equal the corresponding degrees of freedom along the boundary, N 1 i .While for the western and southern boundary the outward unit normal points in the negative x-and y-direction, respectively, and therefore the boundary degree of freedom B 1 j is minus the degree of freedom N 1 i .
Fig. 11.Global nodal degrees of freedom N 0 i for scalar fields in C (Ω) and the boundary degrees of freedom B 0 j for functions in C (∂Ω).

The rot acting on D(Ω) × C(∂Ω)
The space of continuous functions C (Ω) can be restricted to the boundary ∂Ω, which we will refer to as C (∂Ω).Using the mass matrix in this restricted function space, we can set up the dual representation for this one dimensional space as discussed in Section 2.3.This dual representation will be denoted by C(∂Ω).
As we have seen in Section 3.5 the vector fields in D(Ω) only have continuity of normal component over element boundaries; the tangential component will be discontinuous.So D (Ω) lacks the smoothness to apply the conventional rot operator and therefore we will define the rot as, rot : Here we see that the integral on the left is a bilinear form on C(Ω) × C (Ω), the first integral on the right hand side is a bilinear form on D(Ω) × D(Ω), while the boundary integral is a bilinear form over C(∂Ω) × C (∂Ω).
Expanding all variables in their appropriate basis functions reveals that where N 0 is the sparse inclusion matrix containing only the values −1, 0 and 1 that maps the boundary nodal degrees of freedom to the global nodes.The subscript 0 denotes the fact that this matrix is acting on the nodes (geometric dimension 0).
Example 4. In Fig. 11 the small mesh is shown again but this time with the nodal degrees of freedom N 0 i for elements in C (Ω) and the boundary degrees of freedom B 0 j for functions in C (∂Ω).The matrix N 0 , which plays a role in (36), in this case is given by 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 N 0 is a dim(C (Ω)) × dim(C (∂Ω)) matrix.Just like the edges at the outer boundary in Example 3 had an orientation given by the direction of the outward unit normal, all points in the domain have an orientation given by the sense of rotation around the points which we have chosen to be anti-clockwise.Since we endow the boundary degrees of freedom with the same orientation, the matrix N 0 only contains the non-zero value 1.
Lemma 4. For the inclusion matrices N k−1 , N k and the incidence matrix E k,k−1 , 0 < k < d, we have the identity (37) -dimensional degrees of freedom in the primal representation to global (k − 1)-dimensional degrees of freedom in the primal representation.It is an identity for the degrees of freedom on the boundary and sets the internal degrees of freedom to zero.Consequently, the matrix I − N k−1 N ⊺ k−1 sets the degrees of freedom along the boundary to zero and leaves the internal degrees of freedom unchanged.Then E k,k−1 maps the (k − 1)-dimensional degrees of freedom to the k-dimensional degrees of freedom.The k-dimensional degrees of freedom ) N k−1 (φ h ) will be zero along the boundary, since these degrees of freedom are linear combination of (k − 1)-dimensional degrees of freedom along the boundary which are zero.The internal k-dimensional degrees of freedom will, in general, be non-zero.Pre-multiplication with N ⊺ k restricts the k-dimensional degrees of freedom to the boundary and therefore which means that Taking the transpose of these matrices on both sides of the equality sign gives the desired result.□ A more constructive proof of this result will be given in Appendix B.

Extended derivatives of dual variables-the dual de Rham sequence
Straightforward calculation reveals that while div curl g h ≡ 0 for g h ∈ C (Ω), rot ( grad (s h , ŝh ), dh ) will in general not be identically zero for s h ∈ S(Ω), ŝh ∈ D(∂Ω) and dh ∈ C(∂Ω).Because for all c h ∈ C (Ω) we have because the div curl ≡ 0. Eq. (38) shows that the range of grad does not necessarily map into the null space of rot.As a consequence the dual spaces and their vector operations do not form a de Rham sequence.We can remedy this by requiring that Note that the integral of ŝh and (curl c h • n) is a bilinear form over D(∂Ω) × D(∂Ω), while the integral of dh and c h is a bilinear form over C(∂Ω) × C (∂Ω).Since the integrals involve only mutual dual representations, they can be evaluated without reference to the basis functions and only in terms of the degrees of freedom.This gives This implies that we need to have Since Remark 5.The above expression is synonymous to the operation of 1D dual gradient on the elements of ∂Ω which is a 1D manifold.This means that the trace spaces on ∂Ω also satisfy a 1D de Rham sequence.
Definition 11.We now extend the definition of grad in (23) with the condition (41) to a map GRAD : ) −1 are the dual basis functions defined over the boundary ∂Ω.
Definition 12.We extend the rot operator in (35) to a map ROT : Using Definitions 11 and 12 we have Here we used the fact that E 2,1 E 1,0 ≡ 0 and Lemma 4 that Therefore the R( GRAD; S(Ω) × D(∂Ω)) ⊂ K( ROT; D(Ω) × C(∂Ω)), where K( ROT; D(Ω) × C(∂Ω)) denotes here the null space of the ROT operator applied to the space D(Ω) × C(∂Ω).So in order for dual operators to form a de Rham sequence, not only the derivative of the dual variables needs to be defined as was done in (34), but also the derivative of the trace variables, (41).
For dh ∈ D (Ω) × C (∂Ω), we have Corollary 5. Using the norms in Definition 13 and the degrees of freedom for grad(s h , ŝh ) in (34) we have ) .
Using the norms in Definition 13 and the degrees of freedom for rot(d h , dh ) in (36) we have ( ) .

Discrete de Rham sequence
Using the function spaces defined in Section 3.1-Section 3.6, the de Rham sequence for primal spaces is given as and de Rham sequence for the dual spaces is given as

Three-dimensional dual spaces
Similar to Section 3 we can define primal and dual spaces for three dimensional problems.For Ω ⊂ R 3 the de Rham sequence is given by .
We will first define the four finite element spaces for and the corresponding dual spaces G(K ), C(K ), D(K ), S(K ).We will then extend this construction to an arbitrary Ω k ⊂ R 3 .This will be followed by an explanation on the treatment of multiple element case and derivation of differential operators on the dual spaces.
Let ξ i , η j , ζ k , i, j, k = 0, . . ., N, be the GLL-points.Let P be the space of polynomials of degree N as introduced in Example 1 and Q be the space of polynomials of degree N − 1 as introduced in Example 2.

The function space G(K ) ⊂ H 1 (K )
Let G(K ) := P ⊗ P ⊗ P be the tensor product space.Any element p h ∈ G(K ) can be represented by where N 0 (p h ) are the degrees of freedom defined at set of 3D nodes and Ψ 0 (ξ , η, ζ ) are the Lagrange (or nodal) basis through x k given by ) , where Definition 15.The dual degrees of freedom for p h ∈ G(K ) are given by where and the associated dual basis, following Corollary 2, is given by

The function space C (K ) ⊂ H (curl; K )
Consider the polynomial tensor spaces given by, C ξ := Q ⊗ P ⊗ P, C η := P ⊗ Q ⊗ P and C ζ := P ⊗ P ⊗ Q.We define the finite element space in 3D as where the degrees of freedom defined at edges are given by , and the basis that satisfies the Kronecker delta property, in the sense of integrals over the edges, are given by where,

, N
Let p h be any element from G(K ), then following Corollary 1 grad p h expanded in the basis functions of C (K ), where the incidence matrix E 1,0 is for degrees of freedom of a 3D element. 2 This implies that R (grad; G(K )) ⊂ C (K ).

Definition 16. The dual degrees of freedom for c
where and the dual basis are given by, Corollary 2,

The function space D(K ) ⊂ H (div; K )
We define the finite element space as D(K ∈ D(K ) can be expressed as where, the degrees of freedom defined at the surfaces are given by and the basis function that satisfy the Kronecker delta property are given by where, Let c h be any element from C (K ), then the curl operation on c h , following Corollary 1, gives expanded in the basis functions of D(K ).Here the incidence matrix E 2,1 is defined for degrees of freedom of a 3D element. 3  This implies that R (curl; C (K )) ⊂ D(K ).

Definition 17. The dual degrees of freedom and dual basis functions for p
where

The function space S(K )
In this case, we define the finite element space in 3D as S(K 2 The incidence matrix E 1,0 is not the same as used earlier in (5) for degrees of freedom of a 1D element or in (18) for degrees of freedom of a 2D element. 3The incidence matrix E 2,1 is not the same as used earlier in (21) for degrees of freedom of a 2D element.
where the degrees of freedom evaluated over a volume are given by and the basis functions that satisfy the Kronecker-delta property are given by where, Let q h ∈ D(K ), then the div operation on q h gives div q h (ξ , η, ζ expressed in terms of the basis functions of S(K ).This implies that R (div; D(K )) = S(K ).
Definition 18.The dual degrees of freedom and the dual basis for f h ∈ S (K ) are given by, where So far we have introduced finite element spaces on the primal and the dual complex, and the differential operators on the primal complex for a reference element K = [−1, 1] 3 .In the next sections we will present the transformation rules for spaces for an arbitrary Ω k ⊂ R 3 .

Transformation rules for 3D function spaces
In this section we introduce the function spaces for an arbitrary Ω k ⊂ R 3 .The transformations are performed elementby-element provided that the global transformation is continuous.Let Φ k be the diffeomorphism and J be the Jacobian tensor such that,

Transformation rules for spaces
k is the pullback operator.We can then reverse the relation to obtain

Transformation rules for spaces C (Ω k )
Given a vector field v ∈ C (K ), the transformed vector field v ∈ C (Ω k ) on Ω k is given by the following expression Following a similar approach as before, the inverse of this transformation is given by 4.5.3.Transformation rules for spaces D (Ω k ) Given a vector field ū ∈ D(K ), the transformed vector field u ∈ D(Ω k ) on Ω k is given by the following expression It is possible to compute the inverse of this transformation, resulting in

Transformation rules for spaces S (Ω k )
Given a function ḡ ∈ S(K ) on K , we define g ∈ S(Ω k ) on Ω k by The inverse relation can be computed in a similar fashion yielding

Transformation rules for dual function spaces
The construction of dual function spaces for an arbitrary Ω k ⊂ R 3 follows the same steps as that for the reference element.
Let Ψ k (x) for k = 0, 1, 2, 3 be the transformed basis functions for the spaces G (Ω), C (Ω), D (Ω) and S (Ω), respectively.Let the associated mass matrices be given by Then the transformed dual basis and the corresponding dual degrees of freedom are given by As mentioned earlier, we do not explicitly need to construct these basis functions.In practice we simply use the property that the dual basis are bi-orthogonal to the primal basis.To post-process the dual degrees of freedom, they are first converted to primal degrees of freedom by solving a linear system of equations and then reconstructed using primal basis functions.

Finite element spaces for multiple elements in 3D
In this section we introduce the finite element spaces for multiple elements in the 3D case.We follow the same procedure as for the multiple elements in 1D and 2D cases.
The global basis functions in G (Ω) ⊂ H 1 (Ω) are in C 0 (Ω).The nodal degrees of freedom at the boundaries are shared between the elements.Since the dual representation is a linear combination of primal basis functions, G (Ω) also contains continuous functions.
The global basis functions in C (Ω) ⊂ H (curl; Ω) have continuous tangential components and discontinuous normal component between the neighbouring elements.As the dual representations in C (Ω) are linear combinations of basis functions in C (Ω), these functions also have continuous tangential component and discontinuous normal component between the neighbouring elements.
Global representation of vector fields in D (Ω) ⊂ H (div; Ω) have a continuous normal component and discontinuous tangential components between the neighbouring elements.As the dual polynomials in D (Ω) are linear combinations of basis functions in D (Ω) they also have continuous normal component and discontinuous tangential components between the neighbouring elements.
Lastly, piecewise polynomials in S (Ω) ⊂ L 2 (Ω) are discontinuous between elements and therefore the dual representations in S (Ω) are also discontinuous between the elements.Using Definitions 20 and 22, we have ) .
Where we used that E 3,2 E 2,1 ≡ 0 and Lemma 4 that in which case ≡ 0. This proves that the dual curl of the dual gradient of scalar field is zero.

The divergence acting on C (Ω) × G (∂Ω)
Let the restriction of the space G (Ω) to the domain boundary, ∂Ω, be denoted as G (∂Ω), and the corresponding dual space as G (∂Ω).
Definition 23.We define the divergence operator on dual space as div : Since div maps into G (Ω), the left hand side term is a bilinear form in G (Ω) × G (Ω).The first term on the right hand side is a bilinear form on C (Ω) × C (Ω), and the boundary integral term is a bilinear form on G (∂Ω) × G (∂Ω).
Using assembled degrees of freedom, we can write (63) as div where, N 0 is a sparse inclusion matrix only containing the non-zero values −1 and 1, which maps the degrees of freedom of C (∂Ω) to degrees of freedom of C (Ω).
Definition 24.We extend the divergence applied to the dual representation as DIV : Using Definitions 22 and 24 we have ≡ 0, and using Lemma 4 we have that This proves that the dual divergence of the dual curl of a vector field is zero.

Discrete de Rham sequence
Using the function spaces and differential operators defined in Section 4.1-Section 4.7, we can define the primal de Rham's sequence as , and the dual de Rham sequence as

Mixed formulation of the Poisson equation
In this section we want to assess the use of dual spaces for the mixed formulation of the Poisson problem.We present an application of 3D finite element spaces to a constrained minimization problem of the Poisson equation.We will compare the results from the two formulations: (1) with primal spaces only, and (2) with primal and dual spaces.In this application it will be shown that the use of dual spaces can give much sparser systems with a lower condition number without compromising on the accuracy of the solution.Let Ω ⊂ R 3 , then for φ ∈ L 2 (Ω) and q ∈ H(div; Ω) we define the functional for prescribed functions f ∈ L 2 (Ω) and φ ∈ H 1/2 (∂Ω).The optimality conditions for this functional are given by This corresponds to a Poisson equation for φ with Dirichlet boundary condition φ = φ along the boundary.We will consider two different discretizations for this problem.For the first approximation we choose (q h , φ h ) ∈ D(Ω) × S(Ω), we will call this primal-primal formulation, while in the second case we approximate the solution as (q h , φ h ) ∈ D(Ω) × S(Ω), we will call this primal-dual formulation.
If we use this in the variational formulation (72), we get ( where the degrees of freedom of f and the boundary term is evaluated using Definition 4. In (73) the incidence matrix E 3,2 is a sparse topological matrix.All metric properties are contained in the mass matrices M (2) and M (3) .For high order methods, these matrices are dense matrices with large full blocks that destroy the sparsity of the incidence matrix with which they are multiplied.We refer to this formulation as the primal-primal formulation, because both q h and φ h are expanded in primal basis functions.If the mesh is deformed, all sub-matrices in (73) will change because the mass matrices M (k) will change and need to be recomputed.
In this case the discrete system will be ( where the degrees of freedom of f and the prescribed boundary conditions φ are the same as evaluated earlier in (73).
We see that if we expand φ h in terms of dual polynomials, the discrete divergence and gradient blocks in (74) are sparse and no longer depend on the metric of the mesh geometry.It is only the mass matrix M (2) that needs to be constructed for each element separately, unless all elements have the same size, shape and polynomial degree.
Remark 6.We can immediately convert (73) to (74).The mass matrix M (3) in the second row of (73) can be cancelled on both sides of the equation, while the mass matrix M (3) in the first row can be contracted with the degrees of freedom ), but these new unknowns are just the dual degrees of freedom Ñ 0 (φ h ) according to (55).

The discrete inf-sup condition
At the continuous level one establishes well-posedness by showing that the divergence operator, div, from H (div; Ω) into L 2 (Ω) is surjective in which case the Closed Range Theorem states that the Hilbert adjoint of the divergence operator is bounding and therefore injective, [30].At the continuous level surjectivity of the divergence operator is proven through the auxiliary problem: For an arbitrary p ∈ L 2 (Ω) find ψ ∈ H 1 0 (Ω) such that The Lax-Milgram lemma ensures uniqueness of ψ.If we set u p = grad ψ then we have div u p = p and since p ∈ L 2 (Ω) was arbitrary, this shows surjectivity of the divergence operator.
At the finite dimensional level the discrete inf-sup condition is derived in exactly the same way.Let p h be any element from S(Ω), then we need to prove that there is a u h ∈ D(Ω) which is mapped by the divergence operator onto p h .Just as in the continuous setting, we use an auxiliary problem: Find Using the dual basis functions this translates into ) , which has to hold for all vectors Ñ 0 (φ h ) and therefore we have Using ( 52) and (57) we define Ñ 0 (ψ h ) ∈ D(Ω) in which case (75) can be written as or This equation states that for all vectors Ñ 0 (p h ) there exists a vector N 2 (u h ) which is mapped by the discrete divergence, E 3,2 to degrees of freedom in S(Ω) and then mapped by the mass matrix M (3) onto S(Ω).

Table 1
Condition number for the primal-primal formulation and the primal-dual formulation for

Test case
In the following test case we compare the primal-primal formulation with the primal-dual formulation.Using the primal-dual formulation we obtain a sparser algebraic formulation and lower condition number without compromising on quantitative results.Consider the domain (taken from [31]) shown in Fig. 12.The deformed mesh coordinates (x, y, z) ∈ Ω are obtained by transforming the reference coordinates (ξ , η, ζ ) ∈ [−1, 1] 3 with the mapping .
We compare both the formulations with a manufactured solution φ ex = sin(2π x) sin(2π y) sin(2π z) which gives and we use Dirichlet boundary conditions over the entire domain given by φ = φ ex | ∂Ω .
In Fig. 13 we compare the sparsity pattern of primal-primal and primal-dual formulation for the algebraic formulations (73) and (74) for the given test case.In the top-left and the top-right figure we show the sparsity plot for a single element with N = 3.The 8586 non-zero entries of the primal-dual formulation are much less than the 14094 non-zero entries of the primal-primal formulation.In the bottom-left and bottom-right we show similar comparison but now for the multiple element case with total number of elements K el = 2 × 2 × 2 and N = 3.The non-zero entries in the primal-dual formulation 70632 are much less than those in the primal-primal formulation 114696.
In Table 1 we list the condition number of the algebraic system of the two formulations for K el = 1 × 1 × 1 with polynomial degree N = 2, 4, 8.We observe that the condition number for the primal-primal formulation is significantly higher than that of the primal-dual formulation.Also the rate of growth of condition number for increasing N is higher for the primal-primal formulation.In this sense the use of dual polynomials can also be interpreted as a form of inverse type mixed preconditioning [32] ( ) ( ) ( ) .
In the top-left plot of Fig. 14 we show the L 2 -error in the constraint (div q h − f ) and the interpolation error in the right hand side term , for primal-dual formulation.On the x-axis we have the element length (assuming no In the plot at the bottom of Fig. 14 we see the convergence in L 2 -error of φ h .The results from the primal-primal formulation and the primal-dual formulation overlap.We see optimal rate of convergence of errors, of O(N).
In terms of accuracy Fig. 14 shows that the results from primal-primal formulation and primal-dual formulation are equal up to machine precision, and the rate of convergence are optimal for both the formulations.
For the next test problem we choose a pair of div-grad problems that are dual to each other in a continuous setting.We will show that this duality continues to hold true at the discrete level by using the primal-dual representations for these problems.

The Dirichlet-Neumann problems
The mathematical theory of finite elements often makes use of the equivalence of dual problems.In general this equivalence no longer holds at the discrete level.In this section we want to show that this equivalence continues to hold at the finite dimensional level when dual representations are employed.In the proof of Lemma 2.2 in [17], for instance, use is made of an equivalence between a Dirichlet and a Neumann problem.We start with the two problems: Given φ ∈ H 1/2 (∂Ω) 2. The Neumann problem: Find q ∈ H(div; Ω) such that { div q = φ on ∂Ω −grad (div q) + q = 0 in Ω . ( If φ solves the Dirichlet problem (76), then q solves the Neumann problem (77), if and only if φ = div q.Furthermore, it follows that [17] ∥ φ∥ H 1/2 (∂Ω) = ∥φ∥ H 1 (Ω) = ∥q∥ H(div;Ω) .
The finite dimensional problem is to find suitable finite dimensional function spaces, D(Ω) ⊂ H(div; Ω) and S (Ω) ⊂
In Fig. 17a div q h is plotted for the deformed grid with c = 0.3, for N = 8.On the deformed mesh we expect the solution to be less accurate than on the orthogonal mesh, but φ h computed on the same mesh is graphically still identical to Fig. 17a.The difference between div q h and φ h is shown in Fig. 17b.This confirms that for this test case the discrete equivalence (87) holds.
In order to corroborate that the norms ∥ φh ∥ H ( grad;Ω ) and ∥q h ∥ H(div;Ω) are identical for this problem, Table 2 lists these norms on three different meshes, the orthogonal mesh, c = 0.0, the slightly deformed mesh, c = 0.15 and the highly deformed mesh, c = 0.3.This table shows that on all mesh configurations and for all polynomial degrees we have ∥ φh ∥ H ( grad;Ω ) = ∥q h ∥ H(div;Ω) .All the three mesh configurations show convergence to a limiting value of 2.35561.

Eigenvalue problems for the vector Laplace operator
In this section we consider the eigenvalue problem of the vector Laplace operator given by: Let Ω ⊂ R 2 and find λ ∈ R and a non-vanishing u such that −grad div u = λu in Ω div u = 0 on ∂Ω .

Table 2
Norms ∥ φh ∥ H ( grad;Ω ) and ∥q h ∥ H(div;Ω) on three different meshes as a function of the polynomial degree This eigenvalue problem is essentially the same as the eigenvalue problem discussed in [33] apart from the fact that here Neumann boundary conditions are applied.For a more thorough treatment of eigenvalue problems, we refer to [34].
The corresponding weak formulation reads: Find λ ∈ R and u ∈ H(div; Ω) such that For the finite dimensional approximation of this problem, we choose u h ∈ D(Ω) which then leads to the generalized matrix eigenvalue problem Alternatively, we could rephrase this problem in mixed formulation as: Find λ ∈ R and (u, The main difference between (96) and (98) is that (96) gives a large number of eigenvalues λ = 0, because for any ψ ∈ H(curl; Ω), u = curl ψ is an eigensolution associated to λ = 0.The mixed formulation (98) has no eigenvalue λ = 0.
Note that the mixed formulation is essentially an eigenvalue problem for div u and not u.In the finite dimensional setting (98) we take u h Fig. 18.The three meshes on domain [0, π] 2 with 8 × 8 elements.Left: orthogonal mesh; Middle: curved mesh; Right: non-affine mesh.

Orthogonal
Curved Non-affine If we eliminate the degrees of freedom N 1 (u h ) (using Schur complement method) from (99) and rephrase the resulting eigenvalue problem in terms of the dual degrees of freedom for p h by using Ñ 0 (p h ) = M (2) N 2 (p h ), we obtain the dual formulation: Find λ ∈ R and p h ∈ Sh (Ω) such that The eigenvalues calculated for (97) and (100) are identical up to 4 decimal places, therefore we only present results for (100).We use three different mesh configurations (1) orthogonal mesh; (2) curved mesh; (3) non-affine mesh, as shown in Fig. 18.
In Fig. 19 we show the convergence plots for the first five eigenvalues of (100) for Ω = [0 , π] 2 given by λ = 2, 5, 5, 8, 10.On the x-axis we have the element length h = π/K e , where K e is the number of elements in one direction.
On the y-axis we have the absolute error of the eigenvalues.In Fig. 19 along the horizontal rows we have the convergence plots with varying mesh: orthogonal, curved and non-affine, and in that order.Along the vertical columns we have the convergence plots with varying N = 1, 3, 5, and in that order.The numerical eigenvalues are given in Appendix A. The rate of convergence for the eigenvalues with varying meshes and polynomial degrees N are given in Table 3.
For an order N element the optimal convergence rate is given by 2N, [34].For all the three mesh configurations we observe optimal order convergence rates as can be seen in Table 3.

Conclusions
In this paper a dual polynomial space is constructed.The duality pairing between variables from a primal and a dual representation reduces to the vector product between the primal and dual degrees of freedom.The grad, curl and div operators applied to the dual (and primal) representation are topological relations between the dual degrees (and primal degrees) of freedom.These topological relations do not depend on the metric.Furthermore, the addition of differential operations in the boundary ensures that the dual function spaces form a de Rham sequence.The first example where we show the advantage of using dual representation concerns the mixed formulation of the Poisson problem.We derive the inf-sup stability condition in terms of degrees of freedom only.We show optimal convergence for multi-element test case on a curved 3D domain.When a primal-dual formulation is used, two sub-matrices in the mixed formulation become very sparse and the condition number for system matrix is reduced, even though very high order methods are used and these two sub-matrices do not change when the mesh is deformed.The second example shows the equivalence of a Dirichlet-Neumann pair of equations at the discrete level that is not so trivial using primal spaces only.This equivalence Fig. 19.Convergence results for the first five eigenvalues.In the first column we have the results for the orthogonal mesh, in the second column for the curved mesh, and in the third column for the non-affine mesh.In the first row we have the results for N = 1, in the second row for N = 3, and in the third row for N = 5. is proven and illustrated with a test case.And in the last example, we solve the grad-div eigenvalue problem in terms of primal and dual representations, and show optimal convergence on both affine and non-affine meshes.
We have also seen that the use of dual representations allows us to work directly with the degrees of freedom, without explicitly referring to the basis functions.It suffices to make use of its properties.This allows for a direct comparison with staggered finite volume methods.For example, in (74), E 3,2 acts directly on the degrees of freedom of q h and E 3,2 ⊺ acts on the dual degrees of freedom for φ h .
To deal with computationally demanding problems, in future, this work will be extended to the framework of domain decomposition methods.The compatible construction of the discrete trace space presented in this work allows one to set up a fully compatible hybrid finite element formulation.The advantage of a hybrid formulation is that the dual basis functions and degrees of freedom can be defined using element mass matrices only, instead of the global mass matrices presented in this paper.For preliminary work in this direction, see [35].
Furthermore, in this paper the construction of dual polynomial spaces is based on multiplication with mass matrices (or the inverse of mass matrices).These matrices change with change in shape and size of the element.As an improvement, we will present a construction of dual spaces where the mass matrix is also independent of the shape and the size of the element by using wedge product instead of the inner product, see [7, §6.1].To prove the identity we first need to identify some properties of the different matrices.The first point to note is that by changing the numbering of the degrees of freedom we simply permute the rows and columns of the matrices.Naturally, this change in the numbering of the degrees of freedom does not impact the identity.For this reason we are free to choose any numbering that best suites the proof, and the result will hold for any other numbering.Therefore, in this proof we choose the following numbering.contains the contributions that only relate to the interior degrees of freedom.What is important to note here is the zero block that states that the interior degrees of freedom do not contribute to the boundary degrees of freedom, this is essential for the proof.
If we now use (B.2) we can easily see that

Fig. 1 .
Fig. 1.The nodal Lagrange polynomial basis functions and the associated dual polynomials for N = 4.

Fig. 2 .
Fig. 2. The edge polynomial basis functions and the associated dual polynomials for N = 4.

Fig. 3 .
Fig. 3. Global Lagrange basis function for K el = 5 and N = 1 (left) and the corresponding dual basis function (right).

Fig. 4 .
Fig. 4. Global edge basis function for K el = 5 and N = 1 (left) and the corresponding dual function (right).

Fig. 5 .
Fig. 5. Application of the derivative applied to dual variables for K el = 15 and N = 1.

Fig. 6 .
Fig. 6.Application of the derivative applied to dual variables for K el = 100 and N = 1.

3 .
Now we consider the global degrees of freedom, and the global dual degrees of freedom are obtained by pre-multiplying the global primal degrees of freedom with the global mass matrix.The corresponding global dual basis functions are constructed by post-multiplying the global primal basis functions by the inverse of the global mass matrix.Once again, these operations are far too expensive to perform explicitly, but in practice we do not need to construct these dual degrees of freedom and basis functions explicitly and merely use their properties.Global representations of piecewise polynomials in C (Ω) ⊂ H(curl; Ω) are in C 0 (Ω); the nodal degrees of freedom along an element boundary are shared by neighbouring elements.Since the dual representation is a linear combination of the primal basis functions C(Ω) also contains continuous functions.

Fig. 10 .
Fig. 10.Global edge degrees of freedom N 1 i for vector fields in D(Ω) and the boundary degrees of freedom B 1 j for functions in D(∂Ω).

Fig. 12 .
Fig. 12. Test domain for number of elements K el = 3 × 3 × 3 and polynomial order N = 5.The bold lines show the element domains and the thin lines show the mesh generated by GLL quadrature points.

Fig. 14 .
Fig. 14.Top left: L 2 -error in constraint (div q h− f ), and the interpolation error in the right hand side term (f h − f ).Top right: H (div)-error in flux q h .Bottom: L 2 -error in φ h .
For each set of degrees of freedom (associated to geometrical objects of dimension k, with 0 ≤ k ≤ d) we first number the degrees of freedom on the boundary and then the ones on the interior of the domain.Therefore, for each set of degrees of freedom associated to a k-dimensional geometric object, the boundary degrees of freedom are numbered from 1 to d b k and the ones in the interior from d b k + 1 to d k , where d b k is the number of degrees of freedom on the boundary and d k is the total number of degrees of freedom.If we make this choice of numbering, then N k is a d k × d b k matrix with the following formN k = [ D k is a d b k × d bk diagonal matrix filled with only -1 and 1.In the same way the incidence matrixE k,k−1 is a d k × d k−1 of the incidence matrix that only relate boundary degrees of freedom, E k,k−1 bi contains the contribution of the boundary degrees of freedom to the interior, and finally E k,k−1 ii

4 )⊺ k− 1 (
I b k is the d b k × d b k identity matrix.Combining (B.2) and (B.3) we get ( E k,k−1 ) ⊺ N k = together with (B.5), it is straight forward to see that N k−1 N 1, 1].In general, for R d the dual degrees of freedom of N 0 are denoted by Ñ d .
Corollary 2. The dual basis functions are given by 4, 8.