An arbitrary-order discrete de Rham complex on polyhedral meshes: Exactness, Poincar\'e inequalities, and consistency

In this paper we present a novel arbitrary-order discrete de Rham (DDR) complex on general polyhedral meshes based on the decomposition of polynomial spaces into ranges of vector calculus operators and complements linked to the spaces in the Koszul complex. The DDR complex is fully discrete, meaning that both the spaces and discrete calculus operators are replaced by discrete counterparts, and satisfies suitable exactness properties depending on the topology of the domain. In conjunction with bespoke discrete counterparts of $L^2$-products, it can be used to design schemes for partial differential equations that benefit from the exactness of the sequence but, unlike classical (e.g., Raviart--Thomas--N\'ed\'elec) finite elements, are nonconforming. We prove a complete panel of results for the analysis of such schemes: exactness properties, uniform Poincar\'e inequalities, as well as primal and adjoint consistency. We also show how this DDR complex enables the design of a numerical scheme for a magnetostatics problem, and use the aforementioned results to prove stability and optimal error estimates for this scheme.


Introduction
The design of stable and convergent schemes for the numerical approximation of certain classes of partial differential equations (PDEs) requires to reproduce, at the discrete level, the underlying geometric, topological, and algebraic structures. This leads to the notion of compatibility, which can be achieved either in a conforming or non-conforming setting. Relevant  (1.1) where Ω denotes the operator that maps a real value to a constant function over Ω, H 1 (Ω) the space of scalar-valued functions over Ω that are square integrable along with their gradient, H(curl; Ω) (resp. H(div; Ω)) the space of vector-valued functions over Ω that are square integrable along with their curl (resp. divergence). In order to serve as a basis for the numerical approximation of PDEs, discrete counterparts of this sequence of spaces and operators should enjoy the following key properties: (P1) Complex and exactness properties. For the sequence to form a complex, the image of each discrete vector calculus operator should be contained in the kernel of the next one. Moreover, the following exactness properties should be reproduced at the discrete level: Im Ω = Ker grad (since Ω is connected); Im grad = Ker curl if the first Betti number of Ω is zero; Im curl = Ker div if the second Betti number of Ω is zero; Im div = L 2 (Ω) (since we are in dimension three). (P2) Uniform Poincaré inequalities. Whenever a function from a space in the sequence lies in some orthogonal complement of the kernel of the vector calculus operator defined on this space, its (discrete) L 2 -norm should be controlled by the (discrete) L 2 -norm of the operator up to a multiplicative constant independent of the mesh size. (P3) Primal and adjoint consistency. The discrete vector calculus operators should satisfy appropriate commutation properties with the interpolators and their continuous counterparts. Additionally, these operators along with the corresponding (scalar or vector) potentials should approximate smooth fields with sufficient accuracy. Finally, whenever a formal integration by parts is used in the weak formulation of the problem at hand, the vector calculus operators should also enjoy suitable adjoint consistency properties. The notion of adjoint consistency accounts for the failure, in non-conforming settings, to verify global integration by parts formulas exactly.
In the context of Finite Element (FE) approximations, discrete counterparts of the de Rham complex are obtained replacing each space in the sequence with a finite-dimensional subspace. These subspaces are built upon a formulation of magnetostatics, and in [9], which contains a detailed study of the interpolation and stability properties of the low-order VEM spaces. In order to derive an actual discretisation scheme starting from the sequence of virtual spaces, a variational crime involving projections is required. A different approach is pursued in [31,29], where a discrete de Rham (DDR) complex is presented, based on decompositions of full polynomial spaces into the range of vector calculus operators and their L 2 -orthogonal complements. This complex involves discrete spaces and operators that appear, through discrete L 2products, in the formulation of discretisation methods. The analysis in [31,29] focuses on a subset of properties (P1)-(P2) involved in the stability analysis of numerical schemes: local exactness ([31, Theorems 4.1 and 5.1]), global complex property, discrete counterparts of Im grad = Ker curl for domains that do not enclose voids and Im div = L 2 (Ω) ([29, Theorem 3]), as well as Poincaré inequalities for the divergence and the curl ( [29,Theorems 18 and 20,respectively]). This approach completely avoids, both in the construction and in the analysis, the use of (virtual or piecewise polynomial) functions with global regularity, and is closer in spirit to Mimetic Finite Differences and Mixed Hybrid Finite Volume methods.
Regarding consistency properties (P3) for polytopal methods, and starting from low-order methods, results for Compatible Discrete Operator approximations of the Poisson problem based on nodal unknowns can be found in [12]; see in particular the proof of Theorem 3.3 therein, which contains an adjoint consistency result for a gradient reconstructed from vertex values. In the same framework, an adjoint consistency estimate for a discrete curl constructed from edge values can be found in [13,Lemma 2.3]. A rather complete set of consistency results for Mimetic Finite Difference operators can be found in [8], where they appear as intermediate steps in the error analyses of Chapters 5-7. A notable exception is provided by the adjoint consistency of the curl operator, which is not needed in the error estimate of [8,Theorem 7.3] since the authors consider an approximation of the current density based on the knowledge of a vector potential.
Moving to consistency properties for arbitrary-order polytopal methods, error estimates that involve the adjoint consistency of a gradient and the consistency of the corresponding potential have been recently derived in [14] in the framework of the H 1 -conforming Virtual Element method. The same method is considered in [27,Section 3.2], where a different analysis is proposed based on the third Strang lemma. The estimate of the consistency error in [27,Theorem 19] involves, in particular, the adjoint consistency of a discrete gradient reconstructed as the gradient of a scalar polynomial rather than a vector-valued polynomial. We note, in passing, that the concept of adjoint consistency for (discrete) gradients is directly related to the notion of limitconformity in the Gradient Discretisation Method [35], a generic framework which encompasses several polytopal methods. Primal and dual consistency estimates for a discrete divergence and the corresponding vector potential similar (but not identical) to the ones considered here have been established in [32] in the framework of Mixed High-Order methods. Note that these methods (the H 1 -conforming Virtual Element method and the Mixed High-Order method) do not lead to a discrete de Rham complex. In the framework of arbitraryorder compatible discretisations, on the other hand, primal consistency results for the curl appear as intermediate results in [4], where an error analysis for a Virtual Element approximation of magnetostatics is carried out assuming interpolation estimates for three-dimensional vector valued virtual spaces; see Remark 4.4 therein. However, [4] does not establish any adjoint consistency property of the discrete curl (the formulation of magnetostatics considered in this reference does not require it).
Content of the paper. We present a new DDR sequence based, contrary to [31,29], on explicit complements of the ranges of vector calculus operators inspired by the ones used in [4]; these complements are easier to implement, and enable a complete proof of the full set of properties (P1)-(P3). To the best of our knowledge, this is the first time that such a complete panel of results is available for an arbitrary-order polyhedral method compatible with the de Rham complex. The complements considered here are linked to the spaces appearing in the Koszul complex (see, e.g., [2,Chapter 7]) and enjoy two key properties on general polyhedral meshes: they are hierarchical (see Remark 1 below) and their traces on polyhedral faces or edges lie in appropriate polynomial spaces (cf. Proposition 8). These properties make it possible to prove discrete integration by parts formulas for the discrete potentials (see Remarks 8,17,18, and 19 below) which, in turn, are essential to the proof of the adjoint consistency properties.
The key ingredients to establish primal consistency are the polynomial consistency of discrete vector calculus operators along with the corresponding potentials, and their boundedness when applied to the interpolates of smooth functions. The proofs of adjoint consistency, on the other hand, rely on operator-specific techniques and are all grounded in the above-mentioned discrete integration by parts formulas for the corresponding potential reconstructions. Specifically, the key point for the adjoint consistency of the gradient are estimates for local H 1 -like seminorms of the scalar potentials. The adjoint consistency of curl requires, on the other hand, the construction of liftings of the discrete face potentials that satisfy an orthogonality and a boundedness condition. These reconstructions are inspired by the minimal reconstruction operators of [8,Chapter 3], with a key novelty provided by a curl correction which ensures the well-posedness of the reconstruction inside mesh elements and relies on fine results from [25,3].
In order to showcase the theoretical results derived here, we carry out a full convergence analysis for a DDR approximation of magnetostatics. This is, to the best of our knowledge, the first full theoretical result of this kind for arbitrary-order polytopal methods.
The key innovation of the DDR complex presented here, compared to the one in [31,29], precisely lies in the fact that it enables all mathematical results required to prove error estimates for schemes built from this sequence. The only analytical results available in [31,29] are Poincaré inequalities for the curl and the divergence and, as a matter of fact, it seems that the sequence in these references does not satisfy the critical discrete integration by parts formulas mentioned above, and is therefore not amenable to an adjoint consistency analysis.
The rest of the paper is organised as follows. In Section 2 we establish the general setting. Section 3 contains the definition of the DDR sequence along with key intermediate results for the discrete vector calculus operators (including the commutation property in (P3)) and the proof of (P1). In Section 4, we introduce tools for the design and analysis of schemes based on the DDR sequence: polynomial potential reconstructions and L 2 -products on the discrete spaces. Discrete Poincaré inequalities corresponding to (P2) are covered in Section 5. Section 6 contains the statement and proofs of the primal and adjoint consistency results corresponding to (P3). The application of the theoretical tools to the error analysis of a DDR approximation of magnetostatics is considered in Section 7, where numerical evidence supporting the error estimates is also provided. The paper is completed by three appendices. Appendix A contains results on local polynomial spaces including those on the traces of the trimmed spaces constructed from the Koszul complements. Appendix B contains an in-depth and novel study of the div-curl problems defining the curl liftings on polytopal elements: well-posedness, orthogonality and boundedness properties. Finally, Appendix C details the conventions of notation adopted throughout the paper, and lists the main spaces and operators of the DDR complex.

Domain and mesh
For any (measurable) set ⊂ R 3 , we denote by ℎ ≔ sup{| − | : , ∈ } its diameter and by | | its Hausdorff measure. We consider meshes M ℎ ≔ T ℎ ∪ F ℎ ∪E ℎ ∪V ℎ , where: T ℎ is a finite collection of open disjoint polyhedral elements such that Ω = ∈T ℎ and ℎ = max ∈T ℎ ℎ > 0; F ℎ is a finite collection of open planar faces; E ℎ is the set collecting the open polygonal edges (line segments) of the faces; V ℎ is the set collecting the edge endpoints. It is assumed, in what follows, that (T ℎ , F ℎ ) matches the conditions in [28,Definition 1.4], so that the faces form a partition of the mesh skeleton ∈T ℎ . We additionally assume that the polytopes in T ℎ ∪ F ℎ are simply connected and have connected Lipschitz-continuous boundaries. This notion of mesh is related to that of cellular (or CW) complex from algebraic topology; see, e.g., [43,Chapter 7].
The set collecting the mesh faces that lie on the boundary of a mesh element ∈ T ℎ is denoted by F . For any mesh element or face ∈ T ℎ ∪ F ℎ , we denote, respectively, by E and V the set of edges and vertices of .
Throughout the paper, unless otherwise specified, we write in place of ≤ with depending only on Ω, the mesh regularity parameter of [28,Definition 1.9], and the considered polynomial degree. We note that this mesh regularity parameter ∈ (0, 1) is bounded away from 0 when a shape-regular matching simplicial submesh of T ℎ exists such that each ∈ T ℎ is partitioned into simplices of size uniformly comparable to . We also use ≃ as a shorthand for " and ".

Orientation of mesh entities and vector calculus operators on faces
For any face ∈ F ℎ , an orientation is set by prescribing a unit normal vector and, for any mesh element ∈ T ℎ sharing , we denote by ∈ {−1, 1} the orientation of relative to , that is, = 1 if points out of , −1 otherwise. With this choice, is the unit vector normal to that points out of . For any edge ∈ E ℎ , an orientation is set by prescribing the unit tangent vector . Denoting by ∈ F ℎ a face such that ∈ E , its boundary is oriented counter-clockwise with respect to , and we denote by ∈ {−1, 1} the (opposite of the) orientation of relative to that : = 1 if points on in the opposite orientation to , = −1 otherwise. We also denote by the unit vector normal to lying in the plane of such that ( , ) forms a system of right-handed coordinates in the plane of , so that the system of coordinates ( , , ) is right-handed in R 3 . It can be checked that is the normal to , in the plane where lies, pointing out of . For any mesh face ∈ F ℎ , we denote by grad and div the tangent gradient and divergence operators acting on smooth enough functions. Moreover, for any : → R and : → R 2 smooth enough, we define the two-dimensional vector and scalar curl operators such that rot ≔ − /2 (grad ) and rot = div ( − /2 ), where − /2 is the rotation of angle − 2 in the oriented tangent space to .

Lebesgue and Sobolev spaces
Let be a measurable subset of R 3 . We denote by L 2 ( ) the Lebesgue space spanned by functions that are square-integrable over . When is a subset of an -dimensional variety, we will use the boldface notation L 2 ( ) ≔ L 2 ( ) for the space of vector-valued fields over with square-integrable components.
Given an integer and ∈ {Ω} ∪ T ℎ ∪ F ℎ , H ( ) will denote the Sobolev space spanned by square-integrable functions whose partial derivatives of order up to are also square-integrable. Denoting again by the dimension of , we An arbitrary-order discrete de Rham complex on polyhedral meshes 9

Polynomial spaces and decompositions
For a given integer ℓ ≥ 0, P ℓ denotes the space of -variate polynomials of total degree ≤ ℓ, with the convention that P ℓ 0 = R for any ℓ and that P −1 ≔ {0} for any . For any ∈ T ℎ ∪ F ℎ ∪ E ℎ , we denote by P ℓ ( ) the space spanned by the restriction to of the functions in P ℓ 3 . Denoting by 1 ≤ ≤ 3 the dimension of , P ℓ ( ) is isomorphic to P ℓ (see [28,Proposition 1.23]). In what follows, with a little abuse of notation, both spaces are denoted by P ℓ ( ). We additionally denote by ℓ P, the corresponding L 2 -orthogonal projector and let P 0,ℓ ( ) denote the subspace of P ℓ ( ) made of polynomials with zero average over . For the sake of brevity, we also introduce the boldface notations P ℓ ( ) ≔ P ℓ ( ) 3 for all ∈ T ℎ and P ℓ ( ) ≔ P ℓ ( ) 2 for all ∈ F ℎ .
Let again an integer ℓ ≥ 1 be given, and denote by ⊂ E ℎ a collection of edges such that ≔ ∈ forms a connected set. We denote by P ℓ c ( ) ≔ ∈ C 0 ( ) : ( ) | ∈ P ℓ ( ) for all ∈ the space of functions over whose restriction to each edge ∈ is a polynomial of total degree ≤ ℓ and that are continuous at the edges endpoints; these endpoints are collected in the set V ⊂ V ℎ . Denoting by the coordinates vector of a vertex ∈ V ℎ , it can be easily checked that the following mapping is an isomorphism: (2.2) For all ∈ T ℎ ∪ F ℎ , denote by a point inside such that contains a ball centered at of radius ℎ , where is the mesh regularity parameter in [28,Definition 1.9]. For any mesh face ∈ F ℎ and any integer ℓ ≥ 0, we define the following relevant subspaces of P ℓ ( ): (where ⊥ is a shorthand for the rotated vector − /2 ) so that These decompositions of P ℓ ( ) (as well as those of P ℓ ( ) in (2.6) below) result from [2,Corollary 7.4]. Notice that the direct sums in the above expression are not L 2 -orthogonal in general. The L 2 -orthogonal projectors on the spaces (2.3) are, with obvious notation, ℓ G, , c,ℓ G, , ℓ R, , and c,ℓ R, . Similarly, for any mesh element ∈ T ℎ and any integer ℓ ≥ 0 we introduce the following subspaces of P ℓ ( ): Also in this case, the direct sums above are not L 2 -orthogonal in general. The L 2 -orthogonal projectors on the spaces (2.5) are ℓ G, , c,ℓ G, , ℓ R, , and c,ℓ R, .

Remark 1 (Hierarchical complements)
Unlike the L 2 -orthogonal complements considered in [31], the Koszul complements in (2.4) and (2.6) satisfy, for all ∈ T ℎ ∪ F ℎ and all ℓ ≥ 1, Remark 2 (Vector calculus isomorphisms on local polynomial spaces) For any polygon , polyhedron , and polynomial degree ℓ ≥ 0, a consequence of the polynomial exactness [2,Corollary 7.3] is that the following mappings are isomorphisms: An estimate of the norms of the inverses of these differential isomorphisms is provided in Lemma 9 in Appendix A.

Remark 3 (Composition of
Using the definition of the L 2 -orthogonal projectors, and denoting by ℓ P, the L 2 -orthogonal projector on P ℓ ( ), it holds In what follows, we will need the local Nédélec and Raviart-Thomas spaces: These spaces sit between P ℓ−1 ( ) and P ℓ ( ) and are therefore referred to as trimmed in the FE literature. Notice that we have selected the index in (2.12) so as to reflect the maximum polynomial degrees of functions in each space and, as a result, it is shifted by +1 with respect to [31,29].

Recovery operator
As mentioned above, the direct sums in (2.4) and (2.6) are not L 2 -orthogonal.
The following lemma however shows that, for any of these decompositions, a given polynomial can be recovered from its orthogonal projections on each space in the sum.
Lemma 1 (Recovery operator) Let be a Euclidean space, be a subspace of , and c be a complement (not necessarily orthogonal) of in . Let and c be, respectively, the orthogonal projections on and c . Then, the mappings Id − c : → and Id − c : → are isomorphisms. We can therefore define the recovery operator ℜ , c (·, ·) : × c → such that This operator satisfies the following properties: Proof Let us denote by · the norm in . To prove that Id− c is invertible, we show that the mapping c has a norm < 1, which implies The space being finite dimensional, it suffices to see that, for any ∈ with = 1, we have ( c ) < 1. Since is an orthogonal projector, by Pythagoras' theorem we have where the second inequality is a consequence of the fact that c is an orthogonal projection. This concludes the proof that Id − c is an isomorphism. The invertibility of Id − c is obtained similarly, exchanging the roles of and c .
Let us prove the first relation in (2.14). The second follows using the same arguments. We expand (Id − c ) −1 in (2.13) using the series (2.16) (and similarly for (Id − c ) −1 ) to write We have ≥0 ( c ) c = ≥1 ( c ) (we have used = to introduce the pre-factor ) and the operator acting on above therefore reduces to , and returns since ∈ . As for the operator acting on , using again = shows that it is equal to 0. This concludes the proof of the first relation in (2.14).
Fix now ∈ and set ≔ −ℜ , c ( , c ). Applying (2.14) to = and = c shows that = c = 0. Since = ⊕ c , we can write = + with ∈ and c ∈ c , and the definition of the orthogonal projectors on and c therefore yields, with (·, ·) the scalar product on , Hence, = 0 and (2.15) is established. ⊓ ⊔ The following lemma shows that the norm of the recovery operator for the decompositions (2.4) and (2.6) is equivalent to the sum of the norms of its arguments, uniformly in ℎ. In other words, it states that the decompositions are not just algebraic but also topological (uniformly in ℎ). Since the recovery operator will mostly be of interest to us for these pairs of spaces, to alleviate the notations from here on we will write where · denotes the norm induced by · L 2 ( ) on the space of endomorphisms of P ℓ ( ). As a result, Remark 4 (Recovery operator and L 2 -orthogonal complements) When working with L 2 -orthogonal complements to G ℓ ( ) and R ℓ ( ), instead of the Koszul complements in (2.3) and (2.5), the recovery operator is trivial since it consists in the sum of its two arguments (its topological property (2.19) is also obvious). As mentioned in the introduction, however, the Koszul complements enable proofs of commutation and consistency properties that do not seem straightforward with orthogonal complements; the trade-off lies in having to deal with a less trivial recovery operator (although it remains a purely theoretical tool, see Remark 11), whose topological properties are more complex to establish.
Proof 1. Proof of (2.18). We estimate ℓ G, for an element ∈ T ℎ , the other cases being identical. The linear mapping R 3 ∋ ↦ → ℎ −1 ( − ) ∈ R 3 maps onto a polyhedron of diameter 1, transports the spaces P ℓ ( ), G ℓ ( ) and G c,ℓ ( ) on their equivalent over , and simply scales the L 2 -norm of functions. As a consequence, , and we only have to estimate the latter quantity.

Discrete de Rham complex
We define a discrete counterpart of the de Rham complex (1.1). Throughout the rest of this section, we fix an integer ≥ 0 corresponding to the polynomial degree of the discrete sequence. The rules used in the notations are detailed in Appendix C, and the main DDR-related notations are summarised in Table 4.

Discrete spaces
The DDR spaces are spanned by vectors of polynomials whose components, each attached to a mesh entity, are selected in order to: 1) enable the reconstruction of consistent local discrete vector calculus operators and (scalar or vector) potentials in full polynomial spaces of total degree ≤ (or ≤ + 1 for the potentials associated with the gradient); 2) give rise to exact local sequences on mesh elements and faces.
Specifically, the discrete counterparts of H 1 (Ω), H(curl; Ω), H(div; Ω) and L 2 (Ω) are respectively defined as follows: R, ∈ R −1 ( ) and c R, ∈ R c, ( ) for all ∈ T ℎ , R, ∈ R −1 ( ) and c R, ∈ R c, ( ) for all ∈ F ℎ , and ∈ P ( ) for all ∈ E ℎ , G, ∈ G −1 ( ) and c G, ∈ G c, ( ) for all ∈ T ℎ , and ∈ P ( ) for all ∈ F ℎ , and Remark 5 (Component of grad,ℎ on the mesh edge skeleton) By the isomorphism (2.2) with = E ℎ , we can replace the space P +1 c (E ℎ ) in the definition of grad,ℎ by the Cartesian product space ∈E ℎ P −1 ( ) × R V ℎ . This product space is easier to manipulate in practical implementations of the DDR complex.  Table 1: Polynomial components attached to each mesh vertex ∈ V ℎ , edge ∈ E ℎ , face ∈ F ℎ , and element ∈ T ℎ for each of the DDR spaces. The space P +1 c (E ℎ ) in the definition of grad,ℎ has been replaced by  Remark 6 (Components of curl,ℎ and div,ℎ ) For each mesh element or face ∈ T ℎ ∪ F ℎ , the pair of components ( R, , c R, ) of a vector in curl,ℎ defines an element in RT ( ). Similarly, for any ∈ T ℎ , each pair of element components ( G, , c G, ) of a vector in div,ℎ defines an element in N ( ). In the exposition, we prefer to distinguish these components as they play very different roles in the construction.
The polynomial components attached to mesh vertices, edges, faces, and elements for each of the DDR spaces are summarised in Table 1 (notice that we have accounted for Remark 5 for grad,ℎ ). An inspection of Table 1 reveals that its diagonal contains full polynomial spaces on the mesh entities of dimension corresponding to the index of the space in the sequence (with the convection that P ( ) ≔ R for any vertex ∈ V ℎ ). The components collected in the upper triangular portion of the table are non-zero only for ≥ 1, and encode additional information required for the reconstruction of high-order discrete vector calculus operators and potentials. In particular, the complements R c, ( ), R c, ( ), and G c, ( ) complete the information contained, respectively, in the face curl, element curl and tangential trace, and element divergence to construct the corresponding face or element vector potentials; see In what follows, given • ∈ {grad, curl, div} and a mesh entity of dimension greater than or equal to the index of •,ℎ , we denote by •, the restriction of this space to , i.e., •, contains the polynomial components attached to and to all the mesh entities that lie on its boundary.

Remark 7 (Comparison with Raviart-Thomas-Nédélec finite elements)
When is a tetrahedron or a hexahedron, the local spaces in the DDR sequence can be compared to classical (Raviart-Thomas-Nédélec) FE spaces. The number of degrees of freedom in each case for polynomial degrees ∈ {0, 1, 2} (the most commonly used) is reported in Table 2. For ≥ 1, the DDR construction leads to slightly larger spaces on tetrahedra and to significantly smaller spaces on hexahedra. The number of degrees of freedom for the DDR spaces could be further reduced adapting the serendipity techniques of Virtual Elements [6]; this topic is left for a future work.
For codes aiming at general meshes, the implementation of the DDR spaces requires the local (element-by-element) computation of discrete vector operators and potentials, which is an additional cost with respect to traditional FE codes. It should be noticed, however, that: 1) these computations are an embarrassingly parallel task that scales with the number of mesh elements, and are therefore asymptotically less expensive than the resolution of the algebraic systems (see, e.g., Figure 3); 2) this cost can be substantially reduced when dealing with meshes composed of a finite number of element shapes using standard reference element techniques; 3) it possible to combine the FE and DDR approaches on a given mesh (using the former on elements of standard shape and the latter on elements of more general shape, possibly resulting from local mesh refinement).

Interpolators
In the following, for all ℎ ∈ grad,ℎ , we set ≔ ( E ℎ ) | ∈ P +1 ( ). (3.4) The interpolators on the DDR spaces are defined collecting component-wise L 2 -projections. Specifically grad,ℎ : C 0 (Ω) → grad,ℎ is such that, for all ∈ where t, ≔ × ( | × ) denotes the tangent trace of over . Finally, div,ℎ : The restriction of the above interpolators to a mesh entity of dimension larger than or equal to the index of the corresponding space in the sequence (see Table 1) is denoted replacing the subscript ℎ by . Finally, we let P,ℎ : L 2 (Ω) → P (T ℎ ) denote the global L 2 -orthogonal projector such that, for all ∈ L 2 (Ω), ( P,ℎ ) | = P, | for all ∈ T ℎ .

Discrete vector calculus operators
We define in this section the discrete vector calculus operators that appear in the DDR sequence, obtained collecting the L 2 -orthogonal projections of local discrete operators mapping on full polynomial spaces. In what follows, the operators that only appear in the discrete sequence (3.37) through projections are denoted in sans serif font, while those appearing verbatim (without projection) in the sequence are in standard font.

Gradient
The discrete counterpart of the gradient operator in the DDR sequence maps on curl,ℎ , and therefore requires to define local gradients on mesh edges, faces, and elements. For any ∈ E ℎ , the edge gradient : grad, → P ( ) is defined as: For all ∈ grad, = P +1 ( ), where the derivative is taken along according to the orientation of . For any ∈ F ℎ , the face gradient G : grad, → P ( ) is such that, for all = ( , E ) ∈ grad, and all ∈ P ( ), The existence and uniqueness of G in P ( ) follow from the Riesz representation theorem applied to this space equipped with the usual L 2 -product. Similar considerations hold for the other discrete vector calculus operators defined below, and will not be repeated.
Remark 8 (Validity of (3.10)) The relation (3.10) holds, in fact, for any ∈ R ( ) ⊕ R c, +2 ( ). To check it, take ∈ R ( ) and notice that the left-hand side vanishes owing to div = 0, while the right-hand side vanishes owing to the definition (3.9) of G and again div = 0. This means, in particular, that (3.10) holds for any For all ∈ T ℎ , the element gradient G : grad, → P ( ) is defined such that, for all ∈ grad, and all ∈ P ( ), where we have performed an integration by parts on the first term in the right-hand side to pass to the second line.

Lemma 3 (Consistency properties)
The edge, face, and element gradients, and scalar trace satisfy the following consistency properties: Proof Let us prove (3.12). Take ∈ H 1 ( ). For all ∈ P ( ), denoting by 1 and 2 the coordinates of the vertices 1 and 2 of , oriented so that points from 1 to 2 , we have where we have used an integration by parts in the first line, obtained the second equality applying the definition of grad, ∈ P +1 ( ) (which satisfies ( grad, ) ( ) = ( ) for all ∈ V and −1 P, ( grad, ) = −1 P, ) together with ′ ∈ P −1 ( ), and used another integration by parts to conclude. This proves that ( grad, ) ′ = P, ( ′ ).
The equality (3.15) follows from (3.10) written for ∈ R c, ( ) (this choice is made possible by (2.7)) after replacing the full face gradient G by its definition (3.9), simplifying the boundary terms, and invoking again the isomorphism property (2.9), this time with ℓ = .
Finally, (3.16) can be established from (3.14) following the ideas in [31, The following proposition contains a stronger version of [31, Eq. (5.16)], with test function taken in the Nédélec space N +1 ( ) instead of P ( ).
Proposition 1 (Link between element and face gradients) For all ∈ T ℎ and all ( , ) ∈ grad, × N +1 ( ), Proof Writing (3.11) with = curl ∈ P ( ) and recalling the relation div curl = 0, we have The global discrete gradient ℎ : grad,ℎ → curl,ℎ is obtained collecting the projections of each local gradient on the space attached to the corresponding mesh entity: For all (3.18)

Remark 9 (Practical implementation)
In schemes based on the DDR sequence, the discrete gradient (3.18) only appears as an argument of the discrete L 2product on curl,ℎ (see (4.15) below), that is, composed with the scalar trace and potential reconstruction on this space. Thus, leveraging (4.29) below, one never has to implement ℎ , as only the full element gradients (G ) ∈T ℎ are required. Similar considerations hold for the discrete curl defined by (3.31) below (see also [29,Remark 7] on this matter). Notice that this strategy differs from the one often pursued in the context of Virtual Elements, which consists in directly taking the appropriate components of ℎ as degrees of freedom. This difference is linked to the fact that the present construction embeds what could be interpreted in Virtual Element terms as an enhancement, enabling us to reduce the degree of certain internal polynomial components.

Curl
We next consider the DDR counterpart of the curl operator, which maps on div,ℎ and therefore has components at mesh faces and inside mesh elements. For all ∈ F ℎ , the face curl Reasoning as in [31,Proposition 4.3], we get curl, Proposition 2 (Local complex property) Let ∈ F ℎ and denote by : → curl, the restriction to of the global gradient ℎ defined by (3.18). Then, it holds Remark 10 (Two-dimensional complex) The relations (3.13) and (3.21) show that the following two-dimensional sequence forms a complex: Having assumed simply connected, adapting the arguments of [31, Theorem 4.1], one can additionally prove that this complex is exact, that is, Ker = grad, R, Im = Ker , and Im = P ( ).

Proof (Proposition 2)
Let ∈ grad, . Using the definition (3.19) of and where the suppression of −1 R, in the second line is possible since rot ∈ R −1 ( ), the third line is obtained using the definitions (3.9) of G with = rot (additionally noticing that div (rot ) = 0) and (3.8) of , while the conclusion is obtained reasoning as in [ where t, R, ∈ R ( ) is defined, using the isomorphism property (2.8) with Remark 11 (Validity of (3.23)) Observing that both sides of (3.23) vanish when ∈ P 0 ( ), it is inferred that this relation holds in fact for any ∈ P +1 ( ). We also notice that, since R, (by virtue of (3.22) and (2.14)), t, R, can be replaced by t, in the left-hand side of (3.23).
The actual computation of t, does not require the implementation of the recovery operator in the right-hand side of (3.22), but rather hinges on the solution of the following equation: For all ( , Indeed, the test functions of the form ( , 0) with = c R, . These two conditions combined yield (3.22). Similar considerations hold for the three-dimensional potential reconstructions defined in Sections 4.2 and 4.3 below.
Proposition 3 (Properties of the tangential trace) It holds and thus, using (2.14) and Remark 11, we obtain This proves the first relation in (3.24). The second relation is a straightforward consequence of (3.22) and (2.14).

Proof of (3.26). Let
∈ grad, . For all ∈ P +1 ( ), it holds where the first equality follows recalling that = ′ on , integrating by parts on each edge, noting that ( ) ′ | = − rot · (see [31,Eq. (4.20)]), and cancelling out the vertex values that appear twice with opposite sign, while the conclusion is obtained recalling the definition (3.9) of G and observing that div (rot where we have used the inclusion (3.21) in the cancellation, while the conclusion follows from (3.27). This implies t, R, Plugging the above results into (3.22) with = and using the recovery formula The following polynomial consistency property is proved as in [31,Lemma 5.2] (recall the shift of exponent in the notation of the Nédélec space with respect to this reference): Proposition 4 (Link between element and face curls) For all ( , ) ∈ curl, × P +1 ( ), it holds Proof For any ∈ P +1 ( ), writing (3.28) for = grad ∈ P ( ) and using the fact that curl(grad ) = 0 and that (grad ) | × = rot ( | ) for all ∈ F (see [31,Eq. (3.6)]), we infer that · rot ( | ).
Using Remark 11, we arrive at ∫ .

33)
Proposition 5 (Local exactness property) It holds, for all ∈ T ℎ , where denotes the restriction to of the global curl ℎ defined by (3.31) Proof Let us start by proving that = 0 for all ∈ curl, , that is, where we have used grad ∈ G −1 ( ) to introduce the projector −1 G, . Hence, using the definition (3.32) of , we have, for all ∈ P ( ), Since is arbitrary in P ( ), this shows that = 0.
Let us now prove the inclusion Ker( ) ⊂ Im . We fix an element ∈ div, such that = 0 and prove the existence of ∈ curl, such that where the conclusion follows from the relation (3.35) linking volume and face curls. Since grad spans G −1 ( ) as spans P ( ), this proves that By the isomorphism (2.10), this condition defines R, uniquely. ⊓ ⊔

Remark 12 (Variations)
In the spirit of [7, Section 9], one could consider alternatives of the DDR sequence (3.37) obtained varying certain couples of polynomial degrees in such a way as to preserve the exactness properties. Thus one could, e.g., replace R −1 ( ) with R ( ) in the definition (3.2) of curl,ℎ and, correspondingly, G c, ( ) with G c, +1 ( ) in the definition (3.3) of div,ℎ . With these changes, the results of Proposition 5 (and, in particular, (3.36)) remain valid. Assessing the impact such and similar changes on the consistency is, however, more delicate. These developments are left for a future work.

Commutation properties
Lemma 4 (Local commutation properties) It holds, for all ∈ T ℎ , grad, Remark 13 (Global commutation properties) Global commutation properties can be readily inferred from the local ones stated in Lemma 4 when interpolating functions that have sufficient global regularity.

Remark 14 (Role of commutation properties in the design of robust methods)
The commutation properties of Lemma 4 play a key role in the design of discretisation methods robust with respect to the variations of physical parameters. See, e.g., [30] concerning a DDR method for the Reissner-Mindlin plate bending problem robust with respect to plate thickness.

Proof (Lemma 4)
We start by noticing that all the interpolates defined in (3.38)-(3.40) are well-defined under the assumed regularities.
1. Proof of (3.38). By (3.12) it holds, for all ∈ E , Let now ∈ F . Writing the definition (3.9) of G with = grad, | and where we have removed the projectors using their definition in the second equality and we have integrated by parts to conclude. Recalling the definition (2.12) of RT ( ), we can first let span R −1 ( ) to infer and then R c, ( ) to infer c, R, grad is similar: we write the definition (3.11) of G for = grad, and ∈ RT ( ), use property (A.4) along with (3.15) to replace the trace | in each face integral, remove the projectors using their definitions, and integrate by parts. This concludes the proof of (3.38). (3.39). For all ∈ F , by (3.20)

Proof of
where we have used rot t, = (curl ) | · , see [31,Eq. (3.7)]. Writing the definition (3.28) for where we have removed −1 R, using its definition and, recalling (A.5), we have introduced the L 2 -orthogonal projector RT, on RT ( ) in the boundary integral. By (2.12) together with (2.15) written with the choices ( , , c ) = (RT ( ), R −1 ( ), R c, ( )) and (3.24), Plugging this relation into (3.41), we infer where we have used again (A.5) to remove the projector in the boundary term and we have integrated by parts to conclude. Letting curl ), thus concluding the proof of (3.39).

Complex and exactness properties
The properties collected in the following theorem show that the sequence (3.37) forms a (cochain) complex. (3.42). From the consistency properties (3.12), (3.13) and (3.16) of the full gradients and the definition (3.18) of ℎ , it is readily inferred that ℎ grad,ℎ = 0 for all ∈ R, hence grad,ℎ R ⊂ Ker ℎ .

Theorem 1 (Complex property) It holds
To prove converse inclusion Ker ℎ ⊂ grad,ℎ R, let ℎ ∈ grad,ℎ be such that ℎ ℎ = 0. By the definitions (3.18) of ℎ and (3.8) of , this means that ′ = 0 for all ∈ E ℎ , that is, ( E ℎ ) | is constant over . Since Ω has only one connected component, accounting for the single-valuedness of E ℎ at vertices, we thus infer the existence of ∈ R such that E ℎ = . Let now ∈ F ℎ and ∈ R c, ( ). We have c, R, G = 0, and thus where the second equality comes from the definition (3.9) of G , and the conclusion is obtained accounting for the fact that E = and integrating by parts. Since is generic in R c, ( ), recalling the isomorphism (2.9) this implies . As, for all ∈ F ℎ , the previous results give = ( , E ) = grad, , we also have +1 = by (3.14). Similarly, let ∈ T ℎ and ∈ R c, ( ). Writing the definition (3.11) of G for ∈ R c, ( ), and accounting for c, for all ∈ T ℎ , which concludes the proof that ℎ = grad,ℎ .
2. Proof of (3.43). The inclusion (3.43) follows from the local property: . We next notice that it holds, for all ∈ G c, ( ), where we have used the definition (3.28) of C and the property (3.26) of the tangential trace reconstruction in the first equality, the fact that curl ∈ R −1 ( ) to cancel the projector, and the link (3.17) between volume and face gradients to conclude. This shows that c, G, C = 0 and concludes the proof of (3.46).
3. Proof of (3.44). Immediate consequence of (3.34) after observing that and are the restrictions of ℎ and ℎ to , respectively.
4. Proof of (3.45). The inclusion Im ℎ ⊂ P (T ℎ ) is an obvious consequence of the definition (3.33) of the global divergence. To prove the converse inclusion, let ℎ ∈ P (T ℎ ). Since the continuous divergence operator div : The exactness properties of the DDR sequence, depending on the topology of the domain, are collected in the following theorem.

Remark 16 (Meaning of vanishing Betti numbers)
In broad terms, the condition 1 = 0 means that Ω does not have any tunnel, while 2 = 0 means that Ω does not enclose any void. A typical example of Ω that has 1 ≠ 0 is (the interior of) a torus, and an example of Ω with 2 ≠ 0 is a domain enclosed between two concentric spheres. We start by constructing a function E ℎ ∈ P +1 c (E ℎ ) such that = ( E ℎ ) ′ for all ∈ E ℎ . Let 0 , ∈ V ℎ be two distinct mesh vertices of coordinates 0 and , respectively, and denote by E ⊂ E ℎ a set of edges that form a connected path from 0 to (such a path always exists since Ω is connected). By the fundamental theorem of calculus, there is a unique function E ∈ P +1 c (E ) such that E ( 0 ) = 0 and ( E ) ′ | = for all ∈ E , the derivative being taken in the direction of ( E is obtained integrating, in the direction defined on each edge by , the functions ( ) ∈E ). We want to show that the value E ( ) taken at is independent of the choice of the path . To this end, denote by another path from 0 to formed by the edges in E , and denote by − the same path but with reversed orientation. We assume, for the moment, that E and E are disjoint. By similar considerations as before, there exists a unique E ∈ P +1 c (E ) such that E ( 0 ) = 0 and ( E ) ′ | = for all ∈ E . Since 1 = 0 (i.e., there is no "tunnel" crossing Ω), the path ≔ − ′ formed by the edges in E ≔ E ∪ E is a 1-boundary, i.e., there is a set of faces F ⊂ F ℎ giving rise to a connected surface ≔ ∈F such that = . We fix an orientation for and, for all ∈ F , we denote by ∈ {−1, 1} the orientation of relative to . For all ∈ E , there is a unique face ∈ F such that ∈ E , and we let ≔ denote the orientation of relative to . Since = 0 for all ∈ F , it holds where the second equality is obtained from (3.19) with identically equal to 1, while the conclusion follows observing that all the edges that are interior to appear exactly twice in the sum, with opposite signs. Thus, reasoning as in [31,Proposition 4.2], there exists E ∈ P +1 c (E ) such that ( E ) ′ = for all ∈ E , which we can be uniquely identified by additionally prescribing that E ( 0 ) = 0. Under this condition, by uniqueness we infer ( E ) | = ( E ) | for all ∈ E and ( E ) | = ( E ) | for all ∈ E . Since E is continuous at the vertices of , this shows that E ( ) = E ( ). This argument can be extended to paths and such that E ∩ E ≠ ∅, the only difference being that one should reason, in this case, on each connected component of the manifold (corresponding to a "loop" inside the path = − ). Repeating this reasoning for each vertex ∈ V ℎ and all possible paths connecting 0 and , we conclude that there exists a unique E ℎ ∈ P +1 c (E ℎ ) such that E ℎ ( 0 ) = 0 and, recalling the notation (3.4), Let now ∈ F ℎ . We look for a ∈ P −1 ( ) such that ≔ ( , E ) ∈ grad, satisfies = . Plugging = 0 into (3.19), we infer, for all ∈ P ( ), ∫ where the second equality is a consequence of (3.50), while the conclusion follows from (3.27). This shows that R, = −1 R, G ) for all ∈ P −1 ( ).
Let us now enforce c R, = c, where we have used the definition (3.9) of G in the second equality. Recalling the isomorphism (2.9), the above condition defines the sought ∈ P −1 ( ) uniquely.
Writing the definition (3.28) of C with ∈ G c, ( ), using ℎ ∈ Ker ℎ to see that c, G, C = 0, invoking (3.26) to write t, = G and using the link (3.17) between element and face gradients with = , we see that, for any ∈ T ℎ , R, = −1 R, G with ≔ , ( ) ∈F , E for all ∈ P −1 ( ). Proceeding then as for above, we can select ∈ P −1 ( ) to additionally have c R, = c, R, G . This concludes the proof of (3.49). The definitions of the element gradient G and curl C required us to introduce discrete scalar and tangential traces on the mesh faces. In this section, for each ∈ T ℎ and • ∈ {grad, curl, div}, we define discrete potential reconstructions inside acting on the discrete space •, . These potentials have consistency properties, and enable the design of discrete L 2 -inner products on DDR spaces that are also consistent.
Using the polynomial consistency properties G grad, = grad and +1 grad, | = | , valid for all ∈ P +1 ( ) (see (3.16) and (3.14), respectively), the following polynomial consistency of +1 grad, is obtained: Moreover, applying (4.1) to ∈ R c, ( ) (see Remark 1), using the definition (3.11) of G with = , and recalling that div : R c, ( ) → P −1 ( ) is onto, we obtain Recalling the definition (2.17) of ℜ R, (·, ·), the vector potential reconstruction curl, : curl, → P ( ) is such that, for all ∈ curl, , curl, where curl, R, ∈ R ( ) is defined, using the fact that curl : G c, +1 ( ) → R ( ) is an isomorphism (see (2.10)), by ∫ curl, R, Remark 18 (Discrete integration by parts formula for curl, ) Formula (4.5) can be extended to test functions in the Nédélec space N +1 ( ) defined by (2.12). To check it, simply notice that both sides vanish whenever ∈ G ( ) (use curl grad = 0 and the definition (3.28) of C ). Since R, curl, = curl, R, (see (4.4) and (2.14)), we infer that ∫ curl, Apply (4.6) to = curl, with ∈ P ( ), use the consistency properties t, curl, = P, t, = t, and C curl, = curl (see (3.25) and (3.29), respectively), and integrate by parts. Since curl : N +1 ( ) → R ( ) is onto (due to the isomorphism property (2.10)), we obtain the relation . As a consequence, Using similar arguments as in the proof of Proposition 3, we also have Recalling the definition (2.17) of ℜ G, (·, ·), the vector potential reconstruction div, : div, → P ( ) is such that, for all ∈ div, , div, where div, G, ∈ G ( ) is defined by ∫ div, G, Remark 19 (Discrete integration by parts formula for div, ) Observing that div, G, = G, div, (use (2.14)) and that (4.10) holds for any ∈ P +1 ( ) (as can be proved taking constant in and observing that both sides of this equation vanish due to the definition (3.32) of ), we infer ∫ div, Writing this relation for = div, with ∈ RT +1 ( ), observing that div, = P, (div ) = div by (3.40) and P, ( | · ) = | · for all ∈ F by (A.4), and integrating by parts the right-hand side of the resulting expression, we infer G, div, div, by definition of div, , div, and (2.14), we deduce that div, div, = P, ∀ ∈ RT +1 ( ). We now define discrete L 2 -inner products on the DDR spaces. These products are all constructed in a similar way: by assembling local contributions composed of a consistent term based on the potential reconstruction and a stabilisation term that provides a control over the polynomial components on the lower dimensional geometrical objects. Specifically, each L 2 -product (·, ·) grad,ℎ : grad,ℎ × grad,ℎ → R, (·, ·) curl,ℎ : curl,ℎ × curl,ℎ → R, and (·, ·) div,ℎ : div,ℎ × div,ℎ → R is the sum over ∈ T ℎ of its local counterpart defined by: with symmetric, positive semidefinite stabilisation bilinear forms s •, , • ∈ {grad, curl, div} defined as follows: where we recall that the index t, denotes the tangential trace on , and These local stabilisation bilinear forms s •, are polynomially consistent, i.e., they vanish whenever one of their arguments is the interpolate of a polynomial of total degree ≤ + 1 for • = grad, or ≤ for • ∈ {curl, div}. The consistency properties on interpolates of smooth functions of the potential reconstructions and stabilisation forms proved in Section 6 below make these discrete L 2 -products natural candidates for use in the discretisation of PDEs in weak formulation; see the application of DDR to the magnetostatic problem in Section 7.

Component L 2 -norms, bounds, and equivalence properties
The analysis of the stability and consistency properties of the DDR sequence is facilitated by the introduction of L 2 -like norms naturally associated with the choices of polynomial components in the DDR spaces. Specifically we set, for all ℎ ∈ grad,ℎ ,
The next proposition follows from (2.19) and Lemma 9 in Appendix A, in a similar way as in the proof of [29,Proposition 13].
Proposition 6 (Boundedness of local potentials) It holds, for all ∈ T ℎ and all ∈ F , We next establish the equivalence of the norms corresponding to the discrete L 2 -products and the component norms. Proof We only prove the result for • = div, the other cases being similar (see also [29,Proposition 14] for • = curl in the case of orthogonal complements, instead of the Koszul complements (2.3), (2.5)). Let ∈ T ℎ and ∈ div, . By definition (4.16) of the L 2 -product on div, , we have 2 div, = div, where the first inequality follows from a triangle inequality together with the discrete trace inequality ℎ div, Lemma 1.32]), while the conclusion is a consequence of (4.24) together with the definition of ||| ||| div, . This proves in (4.25).
To prove the converse inequality, we start from (4.13) to write where the first inequality follows from the L 2 -boundedness of the orthogonal projectors −1 G, and c, G, together with a triangle inequality, and the conclusion is obtained invoking the same discrete trace inequality as before together with the definition of div, . This proves in (4.25). ⊓ ⊔ Lemma 6 (Boundedness of local interpolators) It holds, for all ∈ T ℎ , ||| grad, ||| grad,

Remark 21 (Boundedness in other norms)
The boundedness of grad, and curl, could easily be stated using norms in larger spaces (typically, 0 ( ) for grad, , and the same spaces on which usual Nédélec interpolators are defined for curl, -see [10, Section 2.5.3]). However, the role of Lemma 6 is to enable primal consistency estimates (Theorem 6); since these estimates require higher regularity on the solutions, the bounds (4.26) and (4.27) stated in non-minimal norms are sufficient to our purpose.

Links between discrete vector potentials and vector calculus operators
In the next proposition, we show that the element gradient and curl can be recovered applying the suitable potential reconstruction to the corresponding discrete vector calculus operator, in a similar way as in (3.26) for the tangential face reconstruction and face gradients.

Proposition 7 (Link between discrete vector potentials and vector calculus operators) For all ∈ T ℎ , it holds
where the conclusion follows from the link between element and face gradients established in Proposition 1. By the isomorphism (2.10) with ℓ = + 1 and since G c, +1 ( ) ⊂ N +1 ( ), this establishes (4.31) and concludes the proof of (4.29).
As before, it therefore remains to analyse the projections on G ( ). Apply (4.11) to = and a generic ∈ P +1 ( ), and use the inclusion Im ⊂ Ker (see Proposition 5) to get ∫ div, where the conclusion is obtained applying the link between element and face curls of Proposition 4. This yields G, div, = G, C , proving (4.30). ⊓ ⊔

Corollary 1 (Bounds on discrete gradients and curl) For all
∈ F ℎ , it holds For all ∈ T ℎ , it holds ∀ ∈ grad, , (4.33) Proof The definitions of |||·||| curl, , |||·||| curl, , and show that the edge gradient contributions in the left-hand sides of (4.32) and (4.33) are bounded by the corresponding right-hand sides. To bound the face and element gradient contributions in the left-hand sides of (4.32) and (4.33), simply apply (4.23) to = and use (3.26) along with (4.29). The estimate (4.34) is established in a similar way, using (4.30). ⊓ ⊔ An arbitrary-order discrete de Rham complex on polyhedral meshes 41

Poincaré inequalities
In this section we state and prove Poincaré-type inequalities for the operators in the DDR sequence. Notice that we consider here the complex without boundary conditions, but one could alternatively consider the complex with (homogeneous) boundary conditions, for which similar inequalities are expected to hold. The details are left for a future work.

Discrete Poincaré inequalities
where, for any mesh element ∈ T ℎ , face ∈ F , and edge ∈ E of vertices 1 and 2 , is the tetrahedron of vertices , , 1 , and 2 . This corresponds to the construction on the dual barycentric mesh of [12, Section 4.1].
We also notice, in passing, that condition (5.1) is not needed when considering the subspace of grad,ℎ with homogeneous boundary conditions. For the sake of completeness, we state in what follows Poincaré inequalities for the curl and the divergence that are easy consequences of the results of [29].

Proof of the discrete Poincaré inequality for the gradient
We first prove a preliminary result, which will also be useful to establish adjoint consistency properties for the discrete gradient operator in Section 6.4.
Lemma 7 (Estimates on local H 1 -seminorms of potentials) For all ∈ F ℎ and all ∈ grad, , it holds For all ∈ T ℎ and all ∈ grad, , it holds (5.6) Proof 1. Proof of (5.5). Let ∈ grad, and define , ∈ R as the average of E over . Introducing grad, , (see (3.14)), using ℎ ≃ ℎ and card(E ) 1, and invoking a discrete trace inequality on +1 − , , we have Since E is continuous, recalling that = ( E ℎ ) | for all ∈ E and using a Poincaré-Wirtinger inequality along followed by the definition (4.21) of |||·||| curl, yields We now turn to the second term in (5.7). Using the isomorphism property (2.9), we select ∈ R c, +2 ( ) such that div = +1 − grad, , . By Lemma 9 in Appendix A, we have The discrete trace inequality of [28,Lemma 1.32] and the consistency property (3.13) of G then yield Hence, applying the definition (3.10) of +1 to − grad, , ∈ grad, , taking above as a test function, and using Cauchy-Schwarz inequalities, we obtain +1 − grad, Simplifying and recalling (4.32) and (5.8), we infer +1 − , ℎ ||| ||| curl, which, plugged together with (5.8) into (5.7), gives the following estimate on the second term in the left-hand side of (5.5): Integrating by parts the definition (3.10) of +1 applied to a generic ∈ P ( ) (see Remark 8) , using Cauchy-Schwarz inequalities, (4.32), a discrete trace inequality, and (5.9) then yields the bound on the first term in the left-hand side of (5.5).
Finally, for 3 , we apply the definition (4.1) of +1  where the second inequality follows from (4.33) and a triangle inequality to write Using discrete trace inequalities and the estimates (5.11) and (5.12) on 1 and 2 , (5.15) leads to ||| 2 curl, . Plugging this bound together with the estimates on 1 and 2 into (5.10) concludes the proof of the bound on the second term in the right-hand side of (5.6). To bound the first term in the left-hand side of (5.6), we proceed as for grad where the conclusion is a consequence of (5.6) followed by the definition (4.21) of the |||·||| curl,ℎ -norm. Let ∈ T ℎ . By definition (4.17) of s grad, we have where the first inequality follows writing in the norms and using triangle and discrete trace inequalities, while the conclusion is obtained invoking (5.5), (5.6), ℎ ≃ ℎ ≃ ℎ and the definition of ||| ||| curl, . Using ℎ 1, summing (5.17) over ∈ T ℎ , and adding the resulting estimate to (5.16) we infer that ℎ grad,ℎ ||| ℎ ℎ ||| curl,ℎ . The Poincaré inequality (5.2) then follows from the norm equivalence (4.25). ⊓ ⊔ 6 Consistency results

Primal consistency
In this section we state consistency results for the discrete potentials, vector calculus operators, stabilisation bilinear forms, and discrete L 2 -products. Because of the nature of the interpolator on curl, (which requires higher regularity of functions), we introduce the following notation: For ∈ T ℎ and ∈ H max( +1,2) ( ), The corresponding global broken seminorm |·| H ( +1,2) ( T ℎ ) is such that, for all The proofs of the following theorems are postponed to Section 6.3.
Theorem 6 (Consistency of the potential reconstructions) It holds, for all ∈ T ℎ , +1 grad, grad, Theorem 7 (Primal consistency of the discrete vector calculus operators) It holds, for all ∈ T ℎ , Theorem 8 (Consistency of stabilisation forms) For all ∈ T ℎ , the stabilisation forms defined by (4.17)-(4.19) satisfy the following consistency properties: The following corollary is a straightforward consequence of Theorems 6 and 8, and of the definitions (4.14)-(4.16) of the discrete L 2 -products. Its proof is therefore omitted.

Adjoint consistency
Whenever a (formal) integration by parts is used to write the weak formulation of a PDE problem underpinning its discretisation, a form of adjoint consistency is required in the convergence analysis. We state here the adjoint consistency of the operators in the DDR sequence (3.37). Since this sequence does not incorporate boundary conditions, the corresponding adjoint consistency will be based on essential (homogeneous) boundary conditions. The regularity requirements will be expressed in terms of the broken Sobolev spaces and norms such that, for any ℓ ≥ 1, The corresponding seminorms for vector-valued functions are denoted using boldface letters, as usual. We denote in what follows by H 1 0 (Ω), H 0 (div; Ω), and H 0 (curl; Ω) the subspaces of H 1 (Ω), H(div; Ω), and H(curl; Ω) spanned by functions whose trace, normal trace, and tangential trace vanish on the boundary Ω of Ω, respectively.

Proof of the adjoint consistency for the gradient
Proof (Theorem 9.) It holds, by definition (4.15) of the local discrete L 2product in curl,ℎ and (4.29), (6.20) Using Remark 17, we have, for all Subtracting this quantity from (6.20), we obtaiñ grad, where is introduced into the boundary term by single-valuedness of the discrete trace, and using | · = 0 whenever ⊂ Ω. Integrating by parts the third term in the right-hand side of the above expression, we obtaiñ grad, ) .

(6.21)
We set = P, and use (6.3) and the approximation properties of P, stated in [28,Theorem 1.45] to see that curl, curl, − P, L 2 ( ) + − P, Using Cauchy-Schwarz inequalities on the integrals and on the stabilisation bilinear form in (6.21), the bound (4.33) together with the norm equivalence An arbitrary-order discrete de Rham complex on polyhedral meshes 53 (4.25), and the consistency property (6.9) of the stabilisation term, we arrive at The conclusion follows from the estimate (5.6), and Cauchy-Schwarz inequalities on the sums. ⊓ ⊔

Proof of the adjoint consistency for the curl
The proof of the adjoint consistency for the curl hinges on liftings defined as solutions of local problems. For any ∈ F ℎ , the face lifting curl, : curl, → H(rot; ) ∩ H(div; ) is such that, for all ∈ curl, , curl, Let now ∈ T ℎ . The curl correction : curl, → H(curl; ) ∩ H(div; ) is such that, for all ∈ curl, , div = − div C in , (6.24a) The curl correction lifts the difference between the face curl and the normal component of the element curl C as a function defined over . Its role is to ensure the well-posedness of the problem defining the element lifting curl, : curl, → H(curl; ) ∩ H(div; ) such that, for all ∈ curl, , In Appendix B we prove that these lifting operators are well-defined, and that they satisfy the following two key properties: -Orthogonality of the face lifting: For all ∈ F ℎ , Proof By the mesh regularity assumption, there is a simplex ⊂ whose inradius is ℎ . Following the arguments in the proof of [28, Lemma 1.25], we infer the norm equivalence Let us take as the Nédélec interpolant in N +1 ( ) of ; can be uniquely extended as an element of N +1 ( ). By the arguments in the proof of [40, Theorem 3.14 and Corollary 3.17], and since ⊂ , it holds We then write, introducing +1 P, and using triangle inequalities, where we have used the approximation property of +1 P, together with the norm equivalence (6.30) in the second line, and concluded by introducing and invoking (6.31) to write +1 P, This concludes the proof of (6.28). The proof of (6.29) is done in a similar way, introducing curl( +1

P,
) and using the approximation property curl − curl( +1 P, Proof (Theorem 10) For all ∈ T ℎ , select ∈ N +1 ( ) given by Lemma 8. Using (4.16) to expand (·, ·) div,ℎ together with (4.30), and recalling (4.6), we see that it holds, for all ℎ ∈ curl,ℎ , . (6.32) Using Cauchy-Schwarz and triangle inequalities, it is readily inferred for the first term where the conclusion follows using the approximation properties (6.4) and (6.28) to bound the first factor, and (4.34) along with the norm equivalence (4.25) to bound the second. For 2 , combining the consistency property (6.10) of s div, with discrete Cauchy-Schwarz inequalities and the definition of the · div,ℎ -norm readily gives For 3 , Cauchy-Schwarz inequalities, the approximation property (6.29), and the definition of the norm · curl,ℎ yield Let us now consider the last term in the right-hand side of (6.32). Since ( ) | × ∈ RT +1 ( ) as a consequence of (A.5) with ℓ = + 1, by (6.26) we can replace t, by curl, in the boundary integral. Using the fact that both curl, and the (rotated) tangential trace of are continuous across interfaces, along with the fact that 1 + 2 = 0 for all ∈ F ℎ between two elements 1 , 2 , and | × = 0 for all ⊂ Ω, we then have where the conclusion follows recalling that curl, = ( curl, ) t, for all ∈ T ℎ and all ∈ F (see (6.25c)), and integrating by parts. Using Cauchy-Schwarz inequalities, it is inferred The approximation properties (6.28)-(6.29) of along with the boundedness (6.27) of curl, yield Plugging (6.33)-(6.36) into (6.32), (6.16) follows. ⊓ ⊔ 6.6 Proof of the adjoint consistency for the divergence Proof (Theorem 11) Combining the definition (6.17) of the adjoint consistency error for the divergence with (4.11) summed over ∈ T ℎ , we infer that it holds, for all ( , ℎ ) as in the theorem and all ℎ ∈ P +1 (T ℎ ) with ≔ ( ℎ ) | for all ∈ T ℎ , where the cancellation of P, is justified by its definition along with ∈ P ( ), while the insertion of into the boundary integral is possible thanks to its single-valuedness at interfaces along with the fact that it vanishes on Ω.
Taking absolute values and using Cauchy-Schwarz inequalities in the righthand side along with ℎ ≃ ℎ for all ∈ T ℎ and all ∈ F , we infer (6.37) Taking ℎ such that = +1 P, | for all ∈ T ℎ and using the approximation properties of the L 2 -orthogonal projector [28,Theorem 1.45], it is inferred that the first factor in the right-hand side of (6.37) is ℎ +1 | | H +2 ( T ℎ ) . Moving to the second factor, we use, for all ∈ T ℎ , [32,Lemma 8] followed by the local seminorm equivalence (4.25) to write ℎ L 2 ( ) ||| ||| div, div, . The same norm equivalence and the definition of the · div, -norm also yields div, The second factor in the right-hand side of (6.37) is therefore ℎ div,ℎ , and the proof is complete. ⊓ ⊔

Convergence analysis for a DDR discretisation of magnetostatics
We analyse in this section the DDR approximation of the following magnetostatics model, in which the unknowns are the magnetic field ∈ H(curl; Ω) and the vector potential ∈ H(div; Ω): The free current belongs to curl H(curl; Ω) and we assume, for the sake of simplicity, that the magnetic permeability is piecewise-constant on the considered meshes, with : Ω → [ − , + ] for some constant numbers 0 < − ≤ + .

Scheme
As shown in [29], a scheme based on the discrete de Rham tools can be written by replacing, in the weak formulation of (7.1), the continuous L 2 -products by discrete ones built on the local products. Denote by the constant value of over ∈ T ℎ and define the bilinear forms ℎ : curl,ℎ × curl,ℎ → R, ℎ , ℎ ∈ curl,ℎ and all ℎ , ℎ ∈ div,ℎ , The discrete problem then reads: Find ℎ ∈ curl,ℎ and ℎ ∈ div,ℎ such that The equations of this problem can be recast in the standard variational form , where A ℎ : ( curl,ℎ × div,ℎ ) 2 → R and L ℎ : curl,ℎ × div,ℎ → R are the bilinear and linear forms, respectively, such that
Proof As shown in the proof of [29,Theorem 10], the exactness of the rightmost part of the sequence (3.37), which holds owing to (3.48) and (3.45), and the Poincaré inequalities for ℎ and ℎ (see Theorems 4 and 5) enable a reproduction of the arguments of the continuous inf-sup condition (see, e.g., [31,Section 2] or [2,Theorem 4.9]) to see that A ℎ satisfies a uniform inf-sup condition with respect to the norm on curl,ℎ × div,ℎ induced by · ,curl,1,ℎ and · div,1,ℎ .

Numerical tests
We present here the results of some numerical tests obtained with the DDR scheme (7.2) for the magnetostatics model (7.1), focusing on comparing outputs obtained using either the complements (2.5), hereafter denoted by (K), or the orthogonal complements of [31,29], denoted by (⊥). Both versions of the DDR complex, and related schemes, have been implemented in the HArDCore3D C++ framework (see https://github.com/jdroniou/HArDCore), using linear algebra facilities from the Eigen3 library (see http://eigen.tuxfamily.org) and the Intel MKL PARDISO library (see https://software.intel.com/en-us/mkl) for the resolution of the global sparse linear system. This solver proved to be the most efficient among those at our disposal. All tests were run on a 16inch 2019 MacBook Pro equipped with an 8-core Intel Core i9 processor (I9-9980HK) and 32Gb of RAM, and running macOS Big Sur version 11.5.1. We consider a constant permeability = 1, and the same exact smooth solution and mesh families as in [29,Section 4.4] for comparability. Figure 1 presents the errors, for various values of , computed in the relative discrete H(curl; Ω) × H(div; Ω) norm: In the case of the Koszul complements, Theorem 12 states that this error should decrease as O(ℎ +1 ) with the mesh size. No such estimate is known for the DDR scheme using orthogonal complements and, due to the lack of key properties of these complements (hierarchical inclusions, structure of traces), it is not clear whether the analysis carried out in the rest of this paper could be adapted to such complements. Nonetheless, the graphs in Figure 1 show that both schemes converge with an order + 1. The errors between (K) and (⊥) are essentially indistinguishable, except for ≥ 1 on tetrahedral meshes, where (⊥) leads to slightly larger errors than (K) -about twice as large on the finest mesh with = 3. The assembly of the (⊥)-DDR scheme requires, for any ∈ T ℎ ∪ F ℎ , to compute bases for the L 2 -orthogonal complements in P ℓ ( ) of G ℓ ( ) and R ℓ ( ), which is done by computing the kernels of local matrices through a full pivot LU algorithm [29,Section 5.1]. On the contrary, in the (K) version, explicit bases for G c,ℓ ( ) and R c,ℓ ( ) can be devised; even though these bases are then orthonormalised to ensure a better numerical stability of the scheme (especially on non-isotropic elements, see the discussion in [28, Section B.1.1] on this topic), the computational cost of creating the polynomial bases in (⊥) can be expected to be larger than in (K). Figure 2 compares the processor times for the two DDR schemes required for (a) the creation of the bases for local polynomial spaces and (b) the model construction (computation of the discrete operators, potentials, and L 2 -products, and global system assembly). We do not compare the linear system resolution times as they are very close for both schemes. In all the cases, the finest mesh of each sequence is considered; see Table 3. A profiling of the code shows that numerical integration is by far the most expensive operation. We therefore include in Figure 2 also a comparison between two integration strategies on general meshes: on one hand, the Homogeneous Numerical Integration of [20]; on the other hand, the use of standard quadratures on a simplicial subdivision of (nonsimplicial) elements. In the left column of Figure 2 we report the total CPU time, which constitutes the most reliable measure to assess performance. Since our code makes use of multi-threading, we also report, in the right column, wall-clock times, which are more representative of real-life performance on the selected architecture. Wall-clock times are subject to outside influences, such as the impact of other processes, and should therefore be regarded with caution.
As expected, when considering standard quadratures on element subdivisions, (K) polynomial bases are faster to create than (⊥) polynomial bases, but not by a large factor (this factor however becomes very large when +1 is required, which is not the case for the scheme (7.2), as the computation of (R +2 ( )) ⊥ (see (4.1)) necessitates to integrate polynomials of degree 2 + 4 over the elements). There is a more pronounced difference when comparing the time for model construction, which is mostly dedicated to the creation of the discrete vector calculus operators and potentials in curl,ℎ and div,ℎ (once these are created, assembling the global linear system itself takes only a small fraction of the total model construction time). Basis construction and model assembly times, on the other hand, basically even out between (K) and (⊥) when considering Homogeneous Numerical Integration, thereby showing the importance of efficient integral computation. Drawing more definitive conclusions is always difficult, as running times highly depend on specific implementation choices, and our implementation is designed for flexibility rather than for efficiency on one given model. The results presented in this section seem  to show, however, that the DDR complex using Koszul complements is not only theoretically better (as it allows for complete consistency analysis and error estimates), but also requires less computational resources, at least when efficient integration is not available in the codes at hand. The comparison of CPU times and wall clock times also confirms that the assembly step strongly benefits from parallel implementations.
To close this section, we briefly assess the evolution of construction and solution times with mesh refinement. On a linear problem such as the one considered here, it is expected that the solution time be larger than the construction time starting from a certain number of elements. To check whether this is the case, we consider the Voronoi mesh family "Voro-small-0" and the polynomial degree = 2 for the sequence. This test is representative of a worst-case scenario for the construction time, since mesh elements are genuinely polyhedral and it is not possible to optimise the construction using standard (reference element) techniques. The plots in Figure 3 show that the asymptotic behaviour of the construction times ("Bases" and "Model") scale linearly with the number of elements (with the former being essentially negligible with respect to the latter), whereas the solution time ("Solve") has a quadratic scaling. The solution time exceeds the construction time starting from the fourth mesh in the sequence, which has 729 elements (a small number for a three-dimensional computation). For meshes of real-life geometries, one can thus expect that the solution time will be the dominating cost (even if more efficient linear solvers become available at some point).

A Results on local spaces
This section collects miscellaneous results on the Koszul complements defined in (2.3) and (2.5), as well as on the trimmed spaces (2.12) obtained from the latter. The first result on traces of Raviart-Thomas and Nédélec functions is known on simplices, see e.g. [10, Proposition 2.3.3]; we however provide a proof on general polyhedra for sake of completeness.

Proof of (A.5). For all
where we have used the vector algebra identity ( × ) × = ( · ) − ( · ) ∀ , , ∈ R 3 (A.6) with = − , = , and = to pass to the second line; to pass to the third line, we have noticed that ( − ) · is constant on for the first term, and we have added ± inside the last parentheses and developed; the last line follows from an application of (A.6) with = − , = , and = . Since P ℓ−1 ( ) = R ℓ−1 ( ) ⊕ R c,ℓ−1 ( ) ⊂ R ℓ−1 ( ) ⊕ R c,ℓ ( ) = RT ℓ ( ), this concludes the proof. ⊓ ⊔ Lemma 9 (Norms of the inverses of local differential isomorphisms) The norms of the inverses of the isomorphisms defined in (2.8)-(2.10) satisfy, for all ∈ F ℎ or ∈ T ℎ , where, above, · denotes the norm of the corresponding isomorphism when its domain and co-domains are endowed with their L 2 -norms, and means that ≤ with depending only on the polynomial degree ℓ and on the mesh regularity parameter.

B Curl lifting
We prove here that the face curl, and element curl, liftings, detailed in Section 6.5, are well defined and satisfy the key properties (6.26) and (6.27).

B.1.1 Existence of
Owing to (6.22b), we look for = rot for some ∈ H 1 ( ). Using the property rot (rot ) = − div (grad ) = −Δ (which stems from (2.1)) and that rot (resp. ) is grad (resp. ) rotated by − /2 in the plane spanned by , we see that (6.22) reduces to the following Neumann problem on : Recalling that is the outer normal, in the plane spanned by , to on , we see that the compatibility condition of this Neumann problem simply amounts to the definition (3.19) of with = 1. There exists therefore a unique ∈ H 1 ( ) solution of this problem with ∫ = 0. Using as a test function in the weak formulation and applying Cauchy-Schwarz inequalities leads to
Since ≥ 0 is strictly positive on a ball, the mapping ( , ) ↦ → ∫ is an inner product on P ( ) and there exists therefore a unique ∈ P ( ) that satisfies this property. This establishes the existence of . Moreover, since = 1 on and · L 2 ( ) and · L 2 ( ) are uniformly equivalent on P ( ) (see the proof of [28, Lemma 1.25]), using = above leads to 2 L 2 ( ) ∫ 2 ≤ t, − L 2 ( ) (div ) −1 L 2 ( ) ||| ||| curl, + L 2 ( ) ℎ L 2 ( ) , where the conclusion follows from a triangle inequality along with the boundedness (4.23) of t, and the estimate (B.2) for the first factor, and Lemma 9 for the second factor. Simplifying, we obtain

B.1.3 Orthogonality property of curl,
We prove here (6.26). Notice first that, since vanishes on and rot grad = 0, by (6.22)  where the second equality follows from (B.4), and the conclusion has been obtained using an integration by parts. This proves that (6.26) holds for ∈ R ( ).

B.2.1 Existence of
Owing to (6.24b), we look for under the form of a potential gradient grad with ∈ H 1 ( ). Equations (6.24a) and (6.24c) then show that must solve the Neumann problem where we recall that is the outer normal to on . The compatibility condition of this problem is which holds true owing to (3.30) with = 1. There exists therefore a unique ∈ H 1 ( ) with ∫ = 0 solution to (B.5). Using as a test function in the weak formulation of (B.5) yields Using the Poincaré-Wirtinger and continuous trace inequalities as we did to obtain (B.2), and recalling that = grad , we infer (B.6) where the conclusion follows from (4.34).
(B.7) The equations (6.25a) and (6.25c) then lead to a curl-curl problem on , whose variational form is: Find ∈ (grad H 1 ( )) ⊥ such that ∫ curl · curl = ∥ ( ) is the set of functions on whose restriction to each face ∈ F belongs to H 1 /2 ( ), and whose tangential traces on the edges are weakly continuous (see [3, 15) and to the fact that (grad H 1 ( )) ⊥ is a closed subspace of H(curl; ), there exists a unique solution to (B.8). We now prove that satisfies (B.8) for all ∈ H(curl; ) = grad H 1 ( ) ⊕ (grad H 1 ( )) ⊥ , which amounts to showing that the right-hand side vanishes whenever = grad for some ∈ H 1 ( ). By density of smooth functions in H 1 ( ), we only need to prove this result for ∈ C ∞ ( ). Plugging = grad in the right-hand side of (B.8), the duality product can be written as standard integrals (since curl, ∈ L 2 ( ) for all ∈ F ) and, integrating by parts, we obtain where we have used (6.24a) to cancel the term in the first equality, and (6.24c) together with integrations by parts on each face in the second equality. Recalling (B.4) and that = 0 if 1 , 2 are the two faces of that share the edge , the right-hand side above vanishes, which shows that (B.8) indeed holds for = grad , and thus for all ∈ H(curl; ).
Since curl, = curl , applying this relation to a generic ∈ C ∞ ( ) and integrating by parts yields (6.25a); using then a generic ∈ C ∞ ( ) and again integrating by parts, we infer (6.25c).

B.2.3 Bound on curl,
We prove here the estimate (6.27). The estimate on curl curl, follows from (6.25a), (4.34) and (B.6). It remains to bound the L 2 -norm of curl, . To do so, we use provided by Lemma  Recalling the construction of each = , for ∈ F , we can extend into a polynomial ∈ P ( ) (for example, by making independent of the coordinate perpendicular to ). We then have, by (B.3), The smooth, compactly supported function can be extended in into such that 0 ≤ ≤ 1, has a compact support in a ball of radius ≃ ℎ that does not touch the faces in F \{ }, and | grad | ℎ −1 . Then, for each ∈ F , the chain rule yields where the second inequality follows from an inverse inequality and (B.12). We then set , = ∈F grad( ) ∈ C ∞ ( ).

C Notations
The notations used in the paper follow these rules: polynomial spaces, subspaces and projections are in curly letters; functions and operators with values in R 2 or R 3 are in boldface; the exponents indicate the maximum polynomial degree of the space or operator; full discrete gradient and curl, which need to be projected to define the operators in the DDR sequence, are in sans serif (and boldface since they are R 3 -valued); spaces, vectors and operators made of components attached to mesh entities of different dimensions are underlined. Table 4 lists the main notations used in the design and analysis of the DDR complex.