discrete de-Rham complex involving a discontinuous finite element space for velocities: the case of periodic straight triangular and Cartesian meshes

Vincent Perrier

► To cite this version:

Vincent Perrier. discrete de-Rham complex involving a discontinuous finite element space for velocities: the case of periodic straight triangular and Cartesian meshes. 2024. hal-04564069

HAL Id: hal-04564069

https://inria.hal.science/hal-04564069

Preprint submitted on 30 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# discrete de-Rham complex involving a discontinuous finite element space for velocities: the case of periodic straight triangular and Cartesian meshes

Vincent Perrier

► To cite this version:

Vincent Perrier. discrete de-Rham complex involving a discontinuous finite element space for velocities: the case of periodic straight triangular and Cartesian meshes. Annales Henri Lebesgue, 2025, 8, pp.417-452. 10.5802/ahl.239. hal-04564069v3

HAL Id: hal-04564069

https://inria.hal.science/hal-04564069v3

Submitted on 14 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
discrete de-Rham complex involving a discontinuous finite
element space for velocities: the case of periodic straight
triangular and Cartesian meshes

Vincent Perrier

Team Cagire, INRIA Bordeaux Sud-Ouest.

Laboratoire de Mathématiques et de leurs applications

Bâtiment IPRA, Université de Pau et des Pays de l'Adour,

Avenue de l'Université, 64013 Pau Cedex

April 30, 2024

Abstract

The aim of this article is to derive discontinuous finite elements vector spaces which can
be put in a discrete de-Rham complex for which an *harmonic gap property* may be proven.

First, discontinuous finite element spaces inspired by classical Nédélec or Raviart-Thomas conforming space are considered, and we prove that by relaxing the normal or tangential constraint, discontinuous spaces ensuring the harmonic gap property can be built.

Then the triangular case is addressed, for which we prove that such a property holds for
the classical discontinuous finite element space for vectors.

On Cartesian meshes, this result does not hold for the classical discontinuous finite element space for vectors. We then show how to use the de-Rham complex found for triangular meshes for enriching the finite element space on Cartesian meshes in order to recover a de-Rham complex, on which the same harmonic gap property is proven.

Contents

<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Introduction</td>
      <td>2</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Notations</td>
      <td>4</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Mesh notations</td>
      <td>4</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Finite element space notations</td>
      <td>5</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Finite element spaces inspired by the conformal case</td>
      <td>6</td>
    </tr>
    <tr>
      <td>4</td>
      <td>The triangular mesh case</td>
      <td>8</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Dimension of the finite elements spaces</td>
      <td>8</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Study of ∇⊥</td>
      <td>13</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Discrete divergence free polynomials on the reference cell</td>
      <td>13</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>ker (∇⊗') and Range (∇⊥)</td>
      <td>16</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Study of ∇⊗'</td>
      <td>17</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Summary on the de-Rham complex</td>
      <td>18</td>
    </tr>
    <tr>
      <td>5</td>
      <td>The case of Cartesian meshes</td>
      <td>20</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Why the Cartesian case is more complicated</td>
      <td>20</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Finite element space definition</td>
      <td>20</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Dimension of the finite element spaces</td>
      <td>20</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Study of ∇⊥⊗'</td>
      <td>21</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Discrete divergence free polynomials on the reference cell</td>
      <td>21</td>
    </tr>
  </tbody>
</table>
discrete de-Rham complex involving a discontinuous finite
element space for velocities: the case of periodic straight
triangular and Cartesian meshes

Vincent Perrier

Team Cagire, INRIA Bordeaux Sud-Ouest.

Laboratoire de Mathématiques et de leurs applications

Bâtiment IPRA, Université de Pau et des Pays de l'Adour,

Avenue de l'Université, 64013 Pau Cedex

January 14, 2025

Abstract

The aim of this article is to derive discontinuous finite elements vector spaces which can
be put in a discrete de-Rham complex for which the matching between the continuous and
discrete cohomology spaces can be proven for periodic meshes.

Then the triangular case is addressed, for which we prove that the same property holds
for the classical discontinuous finite element space for vectors.

On Cartesian meshes, this result does not hold for the classical discontinuous finite element space for vectors. We then show how to use the de-Rham complex found for triangular meshes for enriching the finite element space on Cartesian meshes in order to recover a de-Rham complex, on which the same property is proven.

Contents

<table>
  <tbody>
    <tr>
      <td><b>1 Introduction</b></td>
      <td><b>2</b></td>
    </tr>
    <tr>
      <td><b>2 Notations</b></td>
      <td><b>4</b></td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;2.1 Mesh notations</td>
      <td>4</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;2.2 Finite element space notations</td>
      <td>5</td>
    </tr>
    <tr>
      <td><b>3 Finite element spaces inspired by the conformal case</b></td>
      <td><b>6</b></td>
    </tr>
    <tr>
      <td><b>4 The triangular mesh case</b></td>
      <td><b>8</b></td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;4.1 Dimension of the finite elements spaces</td>
      <td>8</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;4.2 Study of ∇⊥</td>
      <td>13</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;4.3 Discrete divergence free polynomials on the reference cell</td>
      <td>14</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;4.4 ker (∇⊗') and Range (∇⊥)</td>
      <td>16</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;4.5 Study of ∇⊗'</td>
      <td>17</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;4.6 Summary on the de-Rham complex</td>
      <td>18</td>
    </tr>
    <tr>
      <td><b>5 The case of Cartesian meshes</b></td>
      <td><b>19</b></td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;5.1 Why the Cartesian case is more complicated</td>
      <td>19</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;5.2 Finite element space definition</td>
      <td>19</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;5.3 Dimension of the finite element spaces</td>
      <td>21</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;5.4 Study of ∇⊥</td>
      <td>21</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;5.5 Discrete divergence free polynomials on the reference cell</td>
      <td>22</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;5.6 ker (∇⊗') and Range (∇⊥)</td>
      <td>24</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;5.7 Study of ∇⊗'</td>
      <td>25</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;5.8 Summary on the de-Rham complex</td>
      <td>25</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>5.6</td>
      <td>ker(&nabla;<sub>&#x220B;</sub>,') and Range(&nabla;<sup>&perp;</sup>)</td>
      <td>23</td>
    </tr>
    <tr>
      <td>5.7</td>
      <td>Study of &nabla;<sub>&#x220B;</sub>,</td>
      <td>24</td>
    </tr>
    <tr>
      <td>5.8</td>
      <td>Summary on the de-Rham complex</td>
      <td>24</td>
    </tr>
    <tr>
      <td><b>6</b></td>
      <td><b>Conclusion</b></td>
      <td><b>25</b></td>
    </tr>
    <tr>
      <td><b>A</b></td>
      <td><b>Proof for low order quads</b></td>
      <td><b>27</b></td>
    </tr>
  </tbody>
</table>

# 1 Introduction

In this article, we are interested in the de-Rham complex on a two dimensional space $\Omega$:

$$ \Lambda^0(\Omega) \xrightarrow{d} \Lambda^1(\Omega) \xrightarrow{d} \Lambda^2(\Omega), \qquad (1) $$

where $\Lambda^k(\Omega)$ is the set of $k$-differential forms, and $d$ is the exterior derivative. For the sake of simplicity, we suppose that $\Omega$ is the two dimensional torus $T^2$.

In the context of partial differential equations, it is usually convenient to translate the multilinear forms of (1) in terms of proxies. This is achieved by choosing a basis of $\Lambda^1(\Omega)$. If $(e_1, e_2)$ is an orthogonal basis, its dual basis is denoted by $(dx^1, dx^2)$, and this leads usually to the two following situations

* When the basis $(dx^1, dx^2)$ is used for $\Lambda^1$, and $dx^1 \wedge dx^2$ is used for $\Lambda^2$, the exterior derivative between $\Lambda^0$ and $\Lambda^1$ gives a gradient operator on the proxies. In this case, the 0-forms of (1) maps to the set $A$ of the scalar potentials, the 1-forms to a set of vectors $B$, and the 2-forms to a set of scalars. The exterior derivative between $A$ and $B$ is the gradient, $\nabla$:

$$ \begin{array}{rcl} \nabla : & A & \mapsto \mathbf{B} \\ & a & \mapsto \mathbf{b} = (\partial_x a, \partial_y a)^T \end{array} $$

whereas the exterior derivative between $B$ and $C$ is the rotated divergence (or scalar curl, which is obtained by taking the $z$ component of the classical three dimensional curl)

$$ \begin{array}{rcl} \nabla^{\perp} : & \mathbf{B} & \mapsto C \\ & \mathbf{b} & \mapsto c = -\partial_y \mathbf{b}_x + \partial_x \mathbf{b}_y, \end{array} $$

and the diagram (1) can be rewritten in term of proxies as

$$ A \xrightarrow{\nabla} \mathbf{B} \xrightarrow{\nabla^{\perp}} C. \qquad (2) $$

* When the basis $(-dx^2, dx^1)$ is used for $\Lambda^1$, and $dx^1 \wedge dx^2$ is used for $\Lambda^2$, the exterior derivative between $\Lambda^0$ and $\Lambda^1$ gives the rotated gradient (which can also be seen as a curl) on the proxies. In this case, the 0-forms of (1) maps to the set $A$ of scalars (which may be seen as potential vectors by taking the curl of a vector which would have only a $z$ component), the 1-forms to a set of vectors $B$, and the 2-forms to a set of scalars $C$. The exterior derivative between $A$ and $B$ is the rotated gradient, $\nabla^{\perp}$:

$$ \begin{array}{rcl} \nabla^{\perp} : & A & \mapsto \mathbf{B} \\ & a & \mapsto \mathbf{b} = (-\partial_y a, \partial_x a)^T \end{array} $$

whereas the exterior derivative between $B$ and $C$ is the opposite of the two-dimensional divergence

$$ \begin{array}{rcl} \nabla : & \mathbf{B} & \mapsto C \\ & \mathbf{b} & \mapsto c = \partial_x \mathbf{b}_x + \partial_y \mathbf{b}_y, \end{array} $$

and the diagram (1) can be rewritten as

$$ A \xrightarrow{\nabla^{\perp}} \mathbf{B} \xrightarrow{-\nabla} C, \qquad (3) $$

For the sake of simplicity, we will address rather the diagram

$$ A \xrightarrow{\nabla^{\perp}} \mathbf{B} \xrightarrow{\nabla} C, \qquad (4) $$
# 1 Introduction

In this article, we are interested in the de-Rham complex on a two dimensional space $\Omega$:

$$ \Lambda^0(\Omega) \xrightarrow{d} \Lambda^1(\Omega) \xrightarrow{d} \Lambda^2(\Omega), \quad (1) $$

where $\Lambda^k(\Omega)$ is the set of $k$-differential forms, and $d$ is the exterior derivative. For the sake of simplicity, we suppose that $\Omega$ is the two dimensional torus $T^2$.

In the context of partial differential equations, it is usually convenient to translate the multilinear forms of (1) in terms of proxies. This is achieved by choosing a basis of $\Lambda^1(\Omega)$. If $(e_1, e_2)$ is an orthogonal basis, its dual basis is denoted by $(dx^1, dx^2)$, and this leads usually to the two following situations

* When the basis $(dx^1, dx^2)$ is used for $\Lambda^1$, and $dx^1 \wedge dx^2$ is used for $\Lambda^2$, the exterior derivative between $\Lambda^0$ and $\Lambda^1$ gives a gradient operator on the proxies. In this case, the 0-forms of (1) maps to the set $A$ of the scalar potentials, the 1-forms to a set of vectors $B$, and the 2-forms to a set of scalars $C$. The exterior derivative between $A$ and $B$ is the gradient, $\nabla$:

$$ \begin{aligned} \nabla : A &\mapsto B \\ a &\mapsto b = (\partial_x a, \partial_y a)^T \end{aligned} $$

whereas the exterior derivative between $B$ and $C$ is the rotated divergence (or scalar curl, which is obtained by taking the $z$ component of the classical three dimensional curl)

$$ \begin{aligned} \nabla^{\perp} : B &\mapsto C \\ b &\mapsto c = -\partial_y b_x + \partial_x b_y, \end{aligned} $$

and the diagram (1) can be rewritten in term of proxies as

$$ A \xrightarrow{\nabla} B \xrightarrow{\nabla^{\perp}} C. \quad (2) $$

* When the basis $(-dx^2, dx^1)$ is used for $\Lambda^1$, and $dx^1 \wedge dx^2$ is used for $\Lambda^2$, the exterior derivative between $\Lambda^0$ and $\Lambda^1$ gives the rotated gradient (which can also be seen as a curl) on the proxies. In this case, the 0-forms of (1) maps to the set $A$ of scalars (which may be seen as potential vectors by taking the curl of a vector which would have only a $z$ component), the 1-forms to a set of vectors $B$, and the 2-forms to a set of scalars $C$. The exterior derivative between $A$ and $B$ is the rotated gradient, $\nabla^{\perp}$:

$$ \begin{aligned} \nabla^{\perp} : A &\mapsto B \\ a &\mapsto b = (-\partial_y a, \partial_x a)^T \end{aligned} $$

whereas the exterior derivative between $B$ and $C$ is the opposite of the two-dimensional divergence

$$ \begin{aligned} \nabla : B &\mapsto C \\ b &\mapsto c = \partial_x b_x + \partial_y b_y, \end{aligned} $$

and the diagram (1) can be rewritten as

$$ A \xrightarrow{\nabla^{\perp}} B \xrightarrow{-\nabla^{\perp}} C, \quad (3) $$

For the sake of simplicity, we will address rather the diagram

$$ A \xrightarrow{\nabla^{\perp}} B \xrightarrow{\nabla^{\perp}} C, \quad (4) $$

A natural question that arises when considering diagrams such as (2) and (4) is whether the sequence is exact, which can be summarised for (2) as:
A natural question that arises when considering diagrams such as (2) and (4) is whether the sequence is exact, which can be summarised for (2) as :

* Is the kernel of $\nabla$ reduced to 0?

* Is $(\nabla^{\perp}.)$ full rank?

* Do we have $\text{Range}(\nabla) = \text{ker}(\nabla^{\perp}.)$?

In general, the answer of the previous questions is "no", however, it is nearly "yes" in the sense that the dimension of the vectorial spaces $\text{ker} \nabla$, $C/\text{Range}(\nabla^{\perp}.)$ and $\text{ker} \nabla^{\perp}/\text{Range}(\nabla)$ (called the *cohomology spaces*) are finite. Following the Hodge theory, these dimensions equal to the zeroth, first and second Betti numbers (denoted by $b_0$, $b_1$ and $b_2$), which are characteristics of the topology of the domain $\Omega$. We are interested in this article in two-dimensional periodic domain, namely a torus, for which we have $b_0 = b_2 = 1$ and $b_1 = 2$. Also, $\text{ker} \nabla$ and $C/\text{Range}(\nabla^{\perp}.)$ match with the uniform functions (which is a one dimensional vector space, and so is consistent with $b_0 = b_2 = 1$), whereas $\text{ker}(\nabla^{\perp}.)/\text{Range}(\nabla)$ matches with the uniform vectors (which is a two dimensional vector space, so is consistent with $b_1 = 2$).

We are now interested in the discrete counterpart of these properties: once $\Omega$ is discretised by a mesh, do we have a discrete counterpart of (2) and (4)? This indeed exists in the conforming finite element context [6]. For example, for triangular meshes, the following discrete version of (4) involving the space of continuous finite elements $P_k$, the space of Raviart-Thomas finite elements $RT_k$ and the space of discontinuous finite elements $dP_{k-1}$

$$ P_{k+1} \xrightarrow{\nabla^{\perp}} RT_{k+1} \xrightarrow{\nabla} dP_k, \quad (5) $$

whereas a discrete version of (2) can also be derived

$$ P_{k+1} \xrightarrow{\nabla} N_{k+1} \xrightarrow{\nabla^{\perp}} dP_k, \quad (6) $$

where $N_k$ is the space of two dimensional triangular Nédélec first species finite elements. Several properties of the discrete diagrams (5) and (6) are important (see [1, Chap. 5.2.2]):

* the approximation property: this property is ensured if the discrete spaces are correctly approximating the continuous spaces.

* the subcomplex property: this property is a compatibility property between the discrete complex (5) and (6), which should be a subcomplex respectively of the continuous complexes (4) and (2). This is a property that can be ensured only in the conformal case, but not in the discontinuous finite element case in which we are interested in this article.

* the bounded cochain projection property: this property is a commutation property of the diagram involving both the continuous complex, the discrete complex, and the natural projection operators between the discrete and continuous spaces.

These three properties induce another one, the *harmonic gap property* [1, Chap. 5.2.3], which may be summarised as

**Definition 1 (Harmonic gap property).** A discrete diagram ensures the harmonic gap property if the discrete and continuous cohomology spaces are isomorphic.

The finite element exterior calculus has been thoroughly addressed over the last thirty years, first in the electromagnetism context [10, 11, 20, 21], and then extended to the slightly more abstract Hodge Laplacian problem [4, 5], and led to a quite complete theory for conforming finite elements on classical cells (quads, triangles, hexa and tetrahedra) [1]. This type of approximation was extended to polytopal meshes see e.g. [7, 8, 9, 26] for the "Compatible Discrete Operators" framework or [17] for the "Hybrid High Order" method for citing few of these methods. For the classical discontinuous Galerkin methods, as far as we know, few work was considered see however e.g. [22] for recent advances on this topic for the Hodge Laplacian.

In this article, we are interested in finding a discrete counterpart of (2) and (4) when $\Omega$ is meshed with a triangular mesh or with a Cartesian mesh, and with a particular constraint:
* Is the kernel of $\nabla$ reduced to 0?

* Is $(\nabla^{\perp}.)$ full rank?

* Do we have $\text{Range}(\nabla) = \text{ker}(\nabla^{\perp}.)$?

In general, the answer of the previous questions is "no", however, it is nearly "yes" in the sense that the dimension of the vectorial spaces $\text{ker} \nabla$, $C/\text{Range}(\nabla^{\perp}.)$ and $\text{ker} \nabla^{\perp}. / \text{Range}(\nabla)$ (called the cohomology spaces) are finite. Following the Hodge theory, these dimensions equal to the zeroth, first and second Betti numbers (denoted by $b_0$, $b_1$ and $b_2$), which are characteristics of the topology of the domain $\Omega$. We are interested in this article in two-dimensional periodic domain, namely a torus, for which we have $b_0 = b_2 = 1$ and $b_1 = 2$. Also, $\text{ker} \nabla$ and $C/\text{Range}(\nabla^{\perp}.)$ match with the uniform functions (which is a one dimensional vector space, and so is consistent with $b_0 = b_2 = 1$), whereas $\text{ker}(\nabla^{\perp}.) / \text{Range}(\nabla)$ matches with the uniform vectors (which is a two dimensional vector space, so is consistent with $b_1 = 2$).

We are now interested in the discrete counterpart of these properties: once $\Omega$ is discretised by a mesh, do we have a discrete counterpart of (2) and (4)? This indeed exists in the conforming finite element context [6]. For example, for triangular meshes, the following discrete version of (4) involving the space of continuous finite elements $\mathbb{P}_k$, the space of Raviart-Thomas finite elements $\mathbb{RT}_k$ and the space of discontinuous finite elements $\mathbb{dP}_{k-1}$

$$ \mathbb{P}_{k+1} \xrightarrow{\nabla^{\perp}} \mathbb{RT}_{k+1} \xrightarrow{\nabla} \mathbb{dP}_k, \quad (5) $$

whereas a discrete version of (2) can also be derived

$$ \mathbb{P}_{k+1} \xrightarrow{\nabla} \mathbb{N}_{k+1} \xrightarrow{\nabla^{\perp}.} \mathbb{dP}_k, \quad (6) $$

where $\mathbb{N}_k$ is the space of two dimensional triangular Nédélec first species finite elements. Several properties of the discrete diagrams (5) and (6) are important (see [1, Chap. 5.2.2]):

* the approximation property: this property is ensured if the discrete spaces are correctly approximating the continuous spaces.

* the subcomplex property: this property is a compatibility property between the discrete complex (5) and (6), which should be a subcomplex respectively of the continuous complexes (4) and (2), this means for example for the diagram (5), $\nabla^{\perp}(\mathbb{P}_{k+1}) \subset \mathbb{RT}_{k+1}$ and $\nabla \cdot (\mathbb{RT}_{k+1}) \subset \mathbb{dP}_k$ are ensured.

* the bounded cochain projection property: this property means the existence of projection operators, e.g. between the continuous and discrete spaces of (5), that commutes with the exterior derivative, and that is bounded. This property is usually not easy to address, as the canonical interpolant are not bounded, see [4, Section 5.4] for an example of construction of a bounded cochain projection in the conformal case.

The second and third properties, combined with an additional approximation property, induce another property, the isomorphism of cohomology [1, Theorem 5.1]. Another property, the gap between harmonic forms [1, Theorem 5.2] controls the gap between the continuous and discrete harmonic forms. All these properties depend on the definition of the bounded cochain projection, which is not yet defined for the spaces we wish to address, and this is why we focus on a simplified case, the periodic case, for which the cohomology spaces are explicit, and match respectively for their proxies with constant scalar, constant vectors, and constant scalars. In this article, we would like to address the following proposition

**Proposition 1.** The discrete cohomology space matches exactly with the continuous one.

which is expected to hold for periodic meshes.

The finite element exterior calculus has been thoroughly addressed over the last thirty years, first in the electromagnetism context [13, 14, 23, 24], and then extended to the slightly more abstract Hodge Laplacian problem [4, 5], and led to a quite complete theory for conforming finite elements on classical cells (quads, triangles, hexa and tetrahedra) [1]. This type of approximation was extended to polytopal meshes see e.g. [10, 11, 12, 30] for the "Compatible Discrete Operators" framework or [20] for the "Hybrid High Order" method for citing few of these methods. For the classical discontinuous Galerkin methods, as far as we know, few
we want the space of vectors **B** to be discretised with discontinuous finite elements. More precisely, we aim at addressing the problem of the *harmonic gap property* for these discontinuous approximation spaces, but without addressing the approximation, subcomplex or bounded cochain projection properties, which are complicated to address in the nonconformal case, and out of the scope of this paper. As the approximation is nonconforming, the classical differential operators cannot be considered, and discrete differential operators shall be defined. For the sake of simplicity, we will consider differential operators matching with the derivation in the sense of distributions, leading to approximation spaces that are Cartesian products of approximation spaces on different entities of the mesh, similar to the Hybrid High Order framework [17]. The complexes considered are similar to the ones of [25], but include a lower number of degrees of freedom.

The article is organised as follows. In section 2, the notations for the mesh and the finite element space and the discrete differential operators are given. Some enumeration properties of the mesh are also proven in this section. Then in section 3, we recall the results of [25] for a choice of vector finite elements inspired by the conformal case (6),(5). Then, in section 4, we prove that if $\Omega$ is meshed with triangles and if **B** is approximated by the usual discontinuous finite element space, then it can be put in a discrete diagram similar to (2) and (4) where the space $A$ is approximated by the continuous finite element space. The harmonic gap property is proven for this discrete diagram. Then in section 5, the same problem is addressed for Cartesian meshes. We first prove that the harmonic gap property fails for $k = 0$ with the classical piecewise constant finite element vector space. Inspired by the diagram that holds on triangles, we prove that by enriching the classical discontinuous finite element space, the harmonic gap property can be recovered. section 6 is the conclusion.

# 2 Notations

## 2.1 Mesh notations

We denote by $\mathcal{P}$ the set of points of the mesh, by $\mathcal{C}$ the set of cells of the mesh, and by $\mathcal{F}$ the set of the faces of the mesh. For a given entity, for example a cell $c$, we denote by $\mathcal{F}(c)$ the set of faces neighbouring the cell $c$.

Each face is supposed to be oriented, and we denote by $\mathbf{n}_f$ the unit normal to the face $f \in \mathcal{F}$ that is positive. If $\mathbf{u}$ is a vector that is discontinuous through the face $f$, then its jump $[\cdot]$ is defined as

$$[\mathbf{u} \cdot \mathbf{n}_f] = \mathbf{u}_R \cdot \mathbf{n}_f - \mathbf{u}_L \cdot \mathbf{n}_f,$$

where $\mathbf{u}_L$ is the value on the left and $\mathbf{u}_R$ is the value on the right.

**Proposition 1** (Triangular mesh of a torus). For a triangular mesh, if $N$ denotes the number of cells, then

$$\left\{ \begin{array}{l} \#C = N \\ \#\mathcal{F} = \frac{3N}{2} \\ \#\mathcal{P} = \frac{N}{2} \end{array} \right.$$

*Proof.* We remark that the following sum

$$\sum_{f \in \mathcal{F}} \sum_{c \in \mathcal{C}(f)} 1,$$

can be computed in two manners: on one hand, we have two cells per face, so this sum is equal to $2\#\mathcal{F}$. On the other hand, when doing this sum, each cell is visited 3 times (because each cell has three faces), and so the sum is equal to $3N$. This gives $\#\mathcal{F} = \frac{3N}{2}$.

The Euler formula states that

$$\#\mathcal{P} - \#\mathcal{F} + \#\mathcal{C} = 2(1 - g),$$

where $g$ is the genus of the surface. As we are dealing with a two-dimensional domain, with periodic boundary conditions, this is a torus in three dimensions, so that $g = 1$. This leads to

$$\#\mathcal{P} = \#\mathcal{F} - \#\mathcal{C} = \frac{3N}{2} - N = \frac{N}{2}.$$
work was considered see however e.g. [25] for recent advances on this topic for the Hodge Laplacian.

In this article, we are interested in finding a discrete counterpart of (2) and (4) when $\Omega$ is meshed with a triangular mesh or with a Cartesian mesh, and with a particular constraint: we want the space of vectors $\mathbf{B}$ to be discretised with discontinuous finite elements. This constraint is in fact motivated by recent results [22, 27], which suggests that not only the classical discontinuous approximation space on triangles has a structure that allows the preservation of curl constraints for hyperbolic systems, but also that a similar structure exists on quadrangular meshes, see e.g. [28] for the low order quadrangular case. Note that currently, the problem of finding curl or divergence preserving schemes is usually addressed with staggered schemes [34, 7, 8, 9] (based on the discrete de-Rham complex of [1]), which makes the task of limiting for shocks while remaining conservative difficult, because the degrees of freedom are spread on all the entities of the mesh, and also because the notion of local conservation is hard to define.

This article is focused on addressing Proposition 1 for these discontinuous approximation spaces, but without addressing bounded cochain projection property, which are complicated to address in the nonconformal case, and out of the scope of this paper. As the approximation is nonconforming, the classical differential operators cannot be considered, and discrete differential operators shall be defined. For the sake of simplicity, we will consider differential operators matching with the derivation in the sense of distributions, leading to approximation spaces that are Cartesian products of approximation spaces on different entities of the mesh, similar to the Hybrid High Order framework [20]. The complexes considered are similar to the ones of [29], but include a lower number of degrees of freedom.

The article is organised as follows. In section 2, the notations for the mesh and the finite element space and the discrete differential operators are given. Some enumeration properties of the mesh are also proven in this section. Then in section 3, we recall the results of [29] for a choice of vector finite elements inspired by the conformal case (6),(5). Then, in section 4, we prove that if $\Omega$ is meshed with triangles and if $\mathbf{B}$ is approximated by the usual discontinuous finite element space, then it can be put in a discrete diagram similar to (2) and (4) where the space $A$ is approximated by the continuous finite element space. The Proposition 1 is proven for this discrete diagram. Then in section 5, the same problem is addressed for Cartesian meshes. We first prove that Proposition 1 fails for $k = 0$ with the classical piecewise constant finite element vector space. Inspired by the diagram that holds on triangles, we prove that by enriching the classical discontinuous finite element space, Proposition 1 can be recovered. section 6 is the conclusion.

# 2 Notations

## 2.1 Mesh notations

We denote by $\mathcal{P}$ the set of points of the mesh, by $\mathcal{C}$ the set of cells of the mesh, and by $\mathcal{F}$ the set of the faces of the mesh. For a given entity, for example a cell $c$, we denote by $\mathcal{F}(c)$ the set of faces neighbouring the cell $c$, and by $\mathcal{C}(f)$ the set of cells neighbouring the face $f$.

Each face joining points $P$ and $Q$ is supposed to be oriented, and we denote by $\mathbf{n}_f$ the unit normal to the face $f \in \mathcal{F}$ that is positive, namely such that the angle between the vectors $\mathbf{n}_f$ and $\vec{PQ}$ is positive. Then the neighbouring cell of this face such that the normal $\mathbf{n}_f$ is inward is the right cell, and the other is the left cell. If $\mathbf{u}$ is a vector that is discontinuous through the face $f$, then its jump $[\cdot]$ is defined as

$$[\mathbf{u} \cdot \mathbf{n}_f] = \mathbf{u}_R \cdot \mathbf{n}_f - \mathbf{u}_L \cdot \mathbf{n}_f,$$

where $\mathbf{u}_L$ is the value on the left and $\mathbf{u}_R$ is the value on the right.

**Proposition 2** (Triangular mesh of a torus). For a triangular mesh, if $N$ denotes the number of cells, then

$$ \left\{ \begin{array}{l} \#C = N \\ \#\mathcal{F} = \frac{3N}{2} \\ \#\mathcal{P} = \frac{N}{2} \end{array} \right. . $$
□

**Proposition 2** (Cartesian mesh of a torus). For a Cartesian mesh with periodicity, if $N$ denotes the number of cells, then

$$ \begin{cases} \#C = N \\ \#\mathcal{F} = 2N \\ \#\mathcal{P} = N. \end{cases} $$

*Proof.* We remark that the following sum

$$ \sum_{f \in \mathcal{F}} \sum_{c \in \mathcal{C}(f)} 1, $$

can be computed in two manners: on one hand, we have two cells per face, so this sum is equal to $2\#\mathcal{F}$. On the other hand, when doing this sum, each cell is visited 4 times (because each cell has four faces), and so the sum is equal to $4N$. This means that

$$ \#\mathcal{F} = 2N. $$

We are now interested in the following sum

$$ \sum_{p \in \mathcal{P}} \sum_{c \in \mathcal{C}(p)} 1, $$

which is both equal to four times the number of points, but also four times the number of cells. This means that $\#\mathcal{P} = N$. □

## 2.2 Finite element space notations

In this article, we will consider continuous and discontinuous finite element spaces on faces and cells. We adopt notations close of the ones proposed in [6]. We will denote by $\mathbb{P}_k$ the continuous finite element space on triangles. If the finite element space is discontinuous, we will denote it by $d\mathbb{P}_k$. We will also consider vectorial finite element space, $d\mathbb{P}_k$. Last, when needed, we will have to consider finite element spaces on entities of the mesh that are not the cells. In this case, we will then denote by a parenthesis indicating on which entity of the mesh the finite element space is defined. For example, $d\mathbb{P}_k(C)$ is the discontinuous finite element space of degree $k$ defined on the cells, whereas $d\mathbb{P}_k(\mathcal{F})$ is the discontinuous finite element space of degree $k$ defined on the faces.

The continuous and discontinuous finite element spaces are equipped with the classical $L^2$ scalar product and its induced norm. In this article, we will also need to deal with Cartesian products of finite element spaces of type $d\mathbb{P}_i(C) \times d\mathbb{P}_j(\mathcal{F})$, on which we will use the following scalar product

$$ \langle p|q \rangle_{[d\mathbb{P}_i(C) \times d\mathbb{P}_j(\mathcal{F})]} = \sum_{c \in \mathcal{C}} \int_c p_c q_c + \sum_{f \in \mathcal{F}} \int_f p_f q_f. \quad (7) $$

We define the same type of notations by replacing $\mathbb{P}$ by $\mathcal{Q}$ for the case of Cartesian meshes. On Cartesian meshes, we will use enriched versions of $d\mathbb{Q}_k$. For this, we define

$$ \mathcal{Q}_{i,j} = \{p \in \mathbb{R}[x,y] \quad d_x^\circ p \le i \text{ and } d_y^\circ p \le j\}, $$

and can define the following cellwise continuous vectorial finite element space on quads

$$ \widehat{d\mathbb{Q}}_{k}^{\text{div}}(\mathcal{C}) = [(d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k+1,k-1}) \times (d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k-1,k+1})] \oplus \text{Vec} \left( \begin{array}{c} -x^{k+1}y^k \\ x^k y^{k+1} \end{array} \right) \quad (8) $$

that will be suited for the curl/div diagram (4), and the following vectorial finite element space

$$ \widehat{d\mathbb{Q}}_{k}^{\text{curl}}(\mathcal{C}) = [(d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k-1,k+1}) \times (d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k+1,k-1})] \oplus \text{Vec} \left( \begin{array}{c} x^k y^{k+1} \\ x^{k+1} y^k \end{array} \right) \quad (9) $$

that will be suited for the grad/curl diagram (2), and which is nothing but a $\pi/2$ rotation of $\widehat{d\mathbb{Q}}_{k}^{\text{div}}(\mathcal{C})$. Note that the space (8) is the discontinuous version of the space $\mathcal{S}_r$ defined in
Proof. We remark that the following sum

$$ \sum_{f \in \mathcal{F}_c} \sum_{c \in \mathcal{C}(f)} 1, $$

can be computed in two manners: on one hand, we have two cells per face, so this sum is equal to $2\#\mathcal{F}$. On the other hand, when doing this sum, each cell is visited 3 times (because each cell has three faces), and so the sum is equal to $3N$. This gives $\#\mathcal{F} = \frac{3N}{2}$.

The Euler formula states that

$$ \#\mathcal{P} - \#\mathcal{F} + \#\mathcal{C} = 2(1 - g), $$

where $g$ is the genus of the surface. As we are dealing with a two-dimensional domain, with periodic boundary conditions, this is a torus in three dimensions, so that $g = 1$. This leads to

$$ \#\mathcal{P} = \#\mathcal{F} - \#\mathcal{C} = \frac{3N}{2} - N = \frac{N}{2}. $$

**Proposition 3** (Cartesian mesh of a torus). For a Cartesian mesh with periodicity, if $N$ denotes the number of cells, then

$$ \left\{ \begin{array}{l} \#\mathcal{C} = N \\ \#\mathcal{F} = 2N \\ \#\mathcal{P} = N. \end{array} \right. $$

Proof. We remark that the following sum

$$ \sum_{f \in \mathcal{F}_c} \sum_{c \in \mathcal{C}(f)} 1, $$

can be computed in two manners: on one hand, we have two cells per face, so this sum is equal to $2\#\mathcal{F}$. On the other hand, when doing this sum, each cell is visited 4 times (because each cell has four faces), and so the sum is equal to $4N$. This means that

$$ \#\mathcal{F} = 2N. $$

We are now interested in the following sum

$$ \sum_{p \in \mathcal{P}_c} \sum_{c \in \mathcal{C}(p)} 1, $$

which is both equal to four times the number of points, but also four times the number of cells. This means that $\#\mathcal{P} = N$. □

## 2.2 Finite element space notations

In this article, we will consider continuous and discontinuous finite element spaces on faces and cells. We adopt notations close of the ones proposed in [6]. We will denote by $\mathbb{P}_k$ the continuous finite element space on triangles. If the finite element space is discontinuous, we will denote it by $d\mathbb{P}_k$. We will also consider vectorial finite element space, $d\mathbb{P}_k$. Last, when needed, we will have to consider finite element spaces on entities of the mesh that are not the cells. In this case, we will then denote by a parenthesis indicating on which entity of the mesh the finite element space is defined. For example, $d\mathbb{P}_k(C)$ is the discontinuous finite element space of degree $k$ defined on the cells, whereas $d\mathbb{P}_k(\mathcal{F})$ is the discontinuous finite element space of degree $k$ defined on the faces.

The continuous and discontinuous finite element spaces are equipped with the classical $L^2$ scalar product and its induced norm. In this article, we will also need to deal with Cartesian products of finite element spaces of type $d\mathbb{P}_i(C) \times d\mathbb{P}_j(\mathcal{F})$, on which we will use the following scalar product

$$ \langle p|q \rangle_{[d\mathbb{P}_i(C) \times d\mathbb{P}_j(\mathcal{F})]} = \sum_{c \in \mathcal{C}} \int_C p_c q_c + \sum_{f \in \mathcal{F}} \int_f p_f q_f. \quad (7) $$

We define the same type of notations by replacing $\mathbb{P}$ by $\mathbb{Q}$ for the case of Cartesian meshes.
[2, p. 2432] for ensuring optimal approximation of vectors on general quadrangular meshes. It is clear that the cellwise divergence of $\widehat{dQ}_k^{\text{div}}(C)$ or the cellwise curl of $\widehat{dQ}_k^{\text{curl}}(C)$ map to the following finite element space

$$ \widehat{dQ}_{k-1}(C) := dQ_{k-1}(C) + dQ_{k,k-1}(C) + dQ_{k-1,k}(C). $$

Last, we will denote by $\mathbb{K}$ the space of constant elements of the discretisation of the space $A$, $\mathbb{K}$ the space of constant vectors of the discretisation of the space $B$ and $k$ the space of constant elements of $d\mathbb{P}_k(\mathcal{F})$.

### 3 Finite element spaces inspired by the conformal case

In the conformal case, it is known that the harmonic gap property is ensured for the following complexes

$$ \left\{ \begin{array}{ccccc} \mathbb{P}_{k+1} & \xrightarrow{\nabla^{\perp}} & \mathbf{RT}_{k+1}^{\Delta}(C) & \xrightarrow{\nabla} & d\mathbb{P}_k(C) \\ \mathbb{Q}_{k+1} & \xrightarrow{\nabla^{\perp}} & \mathbf{RT}_{k+1}^{\Box}(C) & \xrightarrow{\nabla} & d\mathbb{Q}_k(C) \end{array} \right. , \quad (10) $$

on both the triangular and quadrangular case. The trace of the Raviart-Thomas finite element spaces $\mathbf{RT}_{k+1}^{\Delta}$ and $\mathbf{RT}_{k+1}^{\Box}$ are known to be of degree $k$ in both the triangular and quadrangular case. Therefore, by relaxing the normal continuity constraint, the following complexes may be considered

$$ \left\{ \begin{array}{ccccc} \mathbb{P}_{k+1} & \xrightarrow{\nabla^{\perp}} & \mathbf{dRT}_{k+1}^{\Delta}(C) & \xrightarrow{\nabla_{\mathcal{D}'}} & d\mathbb{P}_k(C) \times d\mathbb{P}_k(\mathcal{F}) \\ \mathbb{Q}_{k+1} & \xrightarrow{\nabla^{\perp}} & \mathbf{dRT}_{k+1}^{\Box}(C) & \xrightarrow{\nabla_{\mathcal{D}'}} & d\mathbb{Q}_k(C) \times d\mathbb{P}_k(\mathcal{F}). \end{array} \right. \quad (11) $$

The discrete maps are defined as follows:

* $\nabla^{\perp}$ is the classical $\nabla^{\perp}$ operator:

$$ \forall c \in \mathcal{C} \quad \forall p \in \mathbb{P}_{k+1} \quad \nabla^{\perp}(p)|_c = \nabla^{\perp}(p)|_c. $$

* $\nabla_{\mathcal{D}'}$ is the divergence where the derivation is taken in the sense of distributions:

$$ \forall \mathbf{u} \in d\mathbb{P}_k \quad \left\{ \begin{array}{ll} \forall c \in \mathcal{C} & \nabla_{\mathcal{D}'}(\mathbf{u})|_c = \nabla \cdot (\mathbf{u}|_c) \\ \forall f \in \mathcal{F} & \nabla_{\mathcal{D}'}(\mathbf{u})|_f = [\mathbf{u} \cdot \mathbf{n}_f]. \end{array} \right. $$

We first compute the dimension of each of the finite element spaces of (11)

**Proposition 3 (Dimension of the finite element spaces).** If the mesh is triangular and periodic, then

$$ \left\{ \begin{array}{l} \dim \mathbb{P}_{k+1} = \frac{N(k+1)^2}{2} \\ \dim \mathbf{dRT}_{k+1}^{\Delta}(C) = N(k+1)(k+3) \\ \dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(C)) = \frac{N(k+1)(k+5)}{2}. \end{array} \right. $$

whereas for a Cartesian periodic mesh,

$$ \left\{ \begin{array}{l} \dim \mathbb{Q}_{k+1} = N(k+1)^2 \\ \dim \mathbf{dRT}_{k+1}^{\Box} = 2N(k+2)(k+1) \\ \dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(C)) = N(k+1)(k+3). \end{array} \right. $$

*Proof.* We first address the triangular case. A $\mathbb{P}_{k+1}$ continuous finite element space has

* 1 degree of freedom on each point.

* $k$ degrees of freedom on each face.

* $\frac{k(k-1)}{2}$ degrees of freedom inside each cell.
On Cartesian meshes, we will use enriched versions of $d\mathbb{Q}_k$. For this, we define

$$ \mathbb{Q}_{i,j} = \{p \in \mathbb{R}[x, y] \quad d_x^\circ p \le i \text{ and } d_y^\circ p \le j\}, $$

where $d_x^\circ$ (resp. $d_y^\circ$) is the degree in $x$ (resp. $y$), and can define the following cellwise continuous vectorial finite element space on quads

$$ \widehat{d\mathbb{Q}}_k^{\text{div}}(C) = \left( \begin{array}{c} d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k+1,k-1} \\ d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k-1,k+1} \end{array} \right) \oplus \text{Vec} \left( \begin{array}{c} -x^{k+1}y^k \\ x^k y^{k+1} \end{array} \right) \quad (8) $$

that will be suited for the curl/div diagram (4), and the following vectorial finite element space

$$ \widehat{d\mathbb{Q}}_k^{\text{curl}}(C) = \left( \begin{array}{c} d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k-1,k+1} \\ d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k+1,k-1} \end{array} \right) \oplus \text{Vec} \left( \begin{array}{c} x^k y^{k+1} \\ x^{k+1} y^k \end{array} \right) \quad (9) $$

that will be suited for the grad/curl diagram (2), and which is nothing but a $\pi/2$ rotation of $\widehat{d\mathbb{Q}}_k^{\text{div}}(C)$. Note that the space (8) is the discontinuous version of the space $S_r$ defined in [2, p. 2432] for ensuring optimal approximation of vectors on general quadrangular meshes. It is clear that the cellwise divergence of $\widehat{d\mathbb{Q}}_k^{\text{div}}(C)$ or the cellwise curl of $\widehat{d\mathbb{Q}}_k^{\text{curl}}(C)$ map to the following finite element space

$$ \widehat{d\mathbb{Q}}_{k-1}(C) := d\mathbb{Q}_{k-1}(C) + d\mathbb{Q}_{k,k-1}(C) + d\mathbb{Q}_{k-1,k}(C). $$

Last, we will denote by $\mathbb{K}$ the space of constant elements of the discretisation of the space $A$, $\mathbb{K}$ the space of constant vectors of the discretisation of the space $B$ and $k$ the space of constant elements of $d\mathbb{P}_k(\mathcal{F})$.

## 3 Finite element spaces inspired by the conformal case

In the conformal case, it is known that the Proposition 1 is ensured for the following complexes

$$ \left\{ \begin{array}{ccccc} \mathbb{P}_{k+1} & \xrightarrow{\nabla^{\perp}} & \mathbf{RT}_{k+1}^{\Delta}(C) & \xrightarrow{\nabla} & d\mathbb{P}_k(C) \\ \mathbb{Q}_{k+1} & \xrightarrow{\nabla^{\perp}} & \mathbf{RT}_{k+1}^{\Box}(C) & \xrightarrow{\nabla} & d\mathbb{Q}_k(C) \end{array} \right. , \quad (10) $$

on both the triangular and quadrangular case. The trace of the Raviart-Thomas finite element spaces $\mathbf{RT}_{k+1}^{\Delta}$ and $\mathbf{RT}_{k+1}^{\Box}$ are known to be of degree $k$ in both the triangular and quadrangular case. Therefore, by relaxing the normal continuity constraint, the following complexes may be considered

$$ \left\{ \begin{array}{ccccc} \mathbb{P}_{k+1} & \xrightarrow{\nabla^{\perp}} & \mathbf{dRT}_{k+1}^{\Delta}(C) & \xrightarrow{\nabla_{\mathcal{D}'}} & d\mathbb{P}_k(C) \times d\mathbb{P}_k(\mathcal{F}) \\ \mathbb{Q}_{k+1} & \xrightarrow{\nabla^{\perp}} & \mathbf{dRT}_{k+1}^{\Box}(C) & \xrightarrow{\nabla_{\mathcal{D}'}} & d\mathbb{Q}_k(C) \times d\mathbb{P}_k(\mathcal{F}). \end{array} \right. \quad (11) $$

The discrete maps are defined as follows:

* $\nabla^{\perp}$ is the classical $\nabla^{\perp}$ operator:

$$ \forall c \in \mathcal{C} \quad \forall p \in \mathbb{P}_{k+1} \quad \nabla^{\perp}(p)|_c = \nabla^{\perp}(p|_c). $$

* $\nabla_{\mathcal{D}'}$ is the divergence where the derivation is taken in the sense of distributions:

$$ \forall \mathbf{u} \in d\mathbb{P}_k \quad \left\{ \begin{array}{ll} \forall c \in \mathcal{C} & \nabla_{\mathcal{D}'}(\mathbf{u})|_c = \nabla \cdot (\mathbf{u}|_c) \\ \forall f \in \mathcal{F} & \nabla_{\mathcal{D}'}(\mathbf{u})|_f = [\mathbf{u} \cdot \mathbf{n}_f]. \end{array} \right. $$

We first compute the dimension of each of the finite element spaces of (11)
Adding all these degrees of freedom leads to

$$
\begin{align*}
\dim \mathbb{P}_{k+1} &= 1 \times \#\mathcal{P} + k\#\mathcal{F} + \frac{k(k-1)}{2} \#\mathcal{C} \\
&= \frac{N}{2} + k \frac{3N}{2} + \frac{k(k-1)^2}{2} N \\
&= \frac{N}{2} (1 + 3k + k(k-1)) \\
&= \frac{N}{2} (k^2 + 2k + 1) \\
\dim \mathbb{P}_{k+1} &= \frac{N(k+1)^2}{2}.
\end{align*}
$$

Then the $(k+1)$th order Raviart-Thomas simplicial finite element is known for having $(k+1)(k+3)$ degrees of freedom (see e.g. [18, Lemma 14.6 p.137]<fnref n="1" />). As the space is discontinuous, this gives

$$
\dim \mathbf{dRT}_{k+1}^{\Delta} = N(k+1)(k+3).
$$

It remains to compute the dimension of $d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(\mathcal{C})$

$$
\begin{align*}
\dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(\mathcal{C})) &= (k+1)\#\mathcal{F} + \frac{(k+1)(k+2)}{2} \#\mathcal{C} \\
&= (k+1) \frac{3N}{2} + \frac{(k+1)(k+2)}{2} N \\
&= \frac{N}{2} (k+1)(k+5).
\end{align*}
$$

We are now interested in the dimension of the finite element spaces for the quadrangular
mesh. We begin by computing the dimension of $\mathbb{Q}_{k+1}$. An element of $\mathbb{Q}_{k+1}$ has

• 1 degree of freedom at each point,

- $k$ degrees of freedom at each face,

- $k^2$ degrees of freedom inside each cell.

Summing all these degrees of freedom and using Proposition 2 gives

$$
\dim \mathbb{Q}_{k+1} = \#\mathcal{P} + k\#\mathcal{F} + k^2\#\mathcal{C} = N + k(2N) + k^2N = N(k+1)^2.
$$

Then the $(k+1)$th order Raviart-Thomas quadrangular finite element is known for having $2(k+1)(k+2)$ degrees of freedom (see e.g. [18, Section 14.5.2 p.142] As the space is discontinuous, this gives

$$
\dim \mathbf{dRT}_{k+1}^{\square} = N(k+1)(k+3).
$$

It remains to compute the dimension of $d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(\mathcal{C})$

$$
\begin{align*}
\dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(\mathcal{C})) &= (k+1)\#\mathcal{F} + (k+1)^2 \#\mathcal{C} \\
&= (k+1)2N + (k+1)^2 N \\
&= N(k+1)(k+3).
\end{align*}
$$

<figure><img src="image_2.png" /></figure>

Properties of the complex (11) was addressed in a more general framework in [25], and
lead in dimension 2 to the following proposition

**Proposition 4.** The discrete diagram (11) ensures the harmonic gap property given in **Def**-
**ition 1.** Moreover, for triangles:

$$
\left\{
\begin{array}{l}
\mathbb{P}_{k+1}/\mathbb{K} = \ker (\nabla^{\perp}) \\
(d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(\mathcal{C}))/k = \text{Range}(\nabla_{\mathcal{D}'\cdot}),
\end{array}
\right.
$$

and for quadrangles

$$
\left\{
\begin{array}{l}
\mathbb{Q}_{k+1}/\mathbb{K} = \ker (\nabla^{\perp}) \\
(d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(\mathcal{C}))/k = \text{Range}(\nabla_{\mathcal{D}'\cdot}).
\end{array}
\right.
$$

<footnote n="1">¹Note that regarding the notations, what is denoted by $\mathbf{RT}_k^{\Delta}$ in [18] is denoted here $\mathbf{RT}_{k+1}^{\Delta}$.</footnote>
**Proposition 4** (Dimension of the finite element spaces). If the mesh is triangular and periodic, then

$$ egin{cases} 
\dim \mathbb{P}_{k+1} = \frac{N(k+1)^2}{2} \\ 
\dim \mathbf{dRT}_{k+1}^{\Delta}(\mathcal{C}) = N(k+1)(k+3) \\ 
\dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(\mathcal{C})) = \frac{N(k+1)(k+5)}{2}. 
\end{cases} $$

whereas for a Cartesian periodic mesh,

$$ egin{cases} 
\dim \mathbb{Q}_{k+1} = N(k+1)^2 \\ 
\dim \mathbf{dRT}_{k+1}^{\square} = 2N(k+2)(k+1) \\ 
\dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(\mathcal{C})) = N(k+1)(k+3). 
\end{cases} $$

*Proof.* We first address the triangular case. A $\mathbb{P}_{k+1}$ continuous finite element space has

* 1 degree of freedom on each point.

* k degrees of freedom on each face.

* $\frac{k(k-1)}{2}$ degrees of freedom inside each cell.

Adding all these degrees of freedom leads to

$$ egin{aligned} 
\dim \mathbb{P}_{k+1} &= 1 \times \#\mathcal{P} + k\#\mathcal{F} + \frac{k(k-1)}{2} \#\mathcal{C} \\ 
&= \frac{N}{2} + k \frac{3N}{2} + \frac{k(k-1)}{2} N \\ 
&= \frac{N}{2} (1 + 3k + k(k-1)) \\ 
&= \frac{N}{2} (k^2 + 2k + 1) \\ 
\dim \mathbb{P}_{k+1} &= \frac{N(k+1)^2}{2}. 
\end{aligned} $$

Then the $(k+1)$th order Raviart-Thomas simplicial finite element is known for having $(k+1)(k+3)$ degrees of freedom (see e.g. [21, Lemma 14.6 p.137]<fnref n="1" />). As the space is discontinuous, this gives

$$ \dim \mathbf{dRT}_{k+1}^{\Delta} = N(k+1)(k+3). $$

It remains to compute the dimension of $d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(\mathcal{C})$

$$ egin{aligned} 
\dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(\mathcal{C})) &= (k+1)\#\mathcal{F} + \frac{(k+1)(k+2)}{2} \#\mathcal{C} \\ 
&= (k+1) \frac{3N}{2} + \frac{(k+1)(k+2)}{2} N \\ 
&= \frac{N}{2} (k+1)(k+5). 
\end{aligned} $$

We are now interested in the dimension of the finite element spaces for the quadrangular mesh. We begin by computing the dimension of $\mathbb{Q}_{k+1}$. An element of $\mathbb{Q}_{k+1}$ has

* 1 degree of freedom at each point,

* k degrees of freedom at each face,

* $k^2$ degrees of freedom inside each cell.

Summing all these degrees of freedom and using Proposition 3 gives

$$ \dim \mathbb{Q}_{k+1} = \#\mathcal{P} + k\#\mathcal{F} + k^2\#\mathcal{C} = N + k(2N) + k^2N = N(k+1)^2. $$

Then the $(k+1)$th order Raviart-Thomas quadrangular finite element is known for having $2(k+1)(k+2)$ degrees of freedom (see e.g. [21, Section 14.5.2 p.142] As the space is discontinuous, this gives

$$ \dim \mathbf{dRT}_{k+1}^{\square} = N(k+1)(k+3). $$

<footnote n="1">¹Note that regarding the notations, what is denoted by $\mathbf{RT}_k^\Delta$ in [21] is denoted here $\mathbf{RT}_{k+1}^\Delta$.</footnote>
The location of the degrees of freedom for this discrete de-Rham complex for Cartesian meshes is summarised in Figure 1, and in Figure 3 for triangles.

By changing the representation of the linear forms, which is equivalent to rotating of $\pi/2$ the vector spaces, the following proposition is also obtained:

**Proposition 5.** The discrete diagram

$$ \left\{ \begin{array}{ccc} \mathbb{P}_{k+1} & \xrightarrow{\nabla} & \mathbf{d}\mathbf{N}_{k+1}^{\Delta}(\mathcal{C}) \\ \mathbb{Q}_{k+1} & \xrightarrow{\nabla} & \mathbf{d}\mathbf{N}_{k+1}^{\Box}(\mathcal{C}) \end{array} \middle| \begin{array}{ccc} \nabla^{\perp}_{\mathcal{D}'} & \xrightarrow{\nabla^{\perp}_{\mathcal{D}'}} & \mathbf{d}\mathbb{P}_k(\mathcal{C}) \times \mathbf{d}\mathbb{P}_k(\mathcal{F}) \\ \nabla^{\perp}_{\mathcal{D}'} & \xrightarrow{\nabla^{\perp}_{\mathcal{D}'}} & \mathbf{d}\mathbb{Q}_k(\mathcal{C}) \times \mathbf{d}\mathbb{P}_k(\mathcal{F}) \end{array} \right. $$

where $\nabla^{\perp}_{\mathcal{D}'}$ is $\nabla^{\perp}$. in the sense of distributions, ensures the harmonic gap property given in Definition 1. Moreover, for triangles:

$$ \left\{ \begin{array}{l} \mathbb{P}_{k+1}/\mathbb{K} = \ker(\nabla) \\ (\mathbf{d}\mathbb{P}_k(\mathcal{F}) \times \mathbf{d}\mathbb{P}_k(\mathcal{C}))/k = \text{Range}(\nabla^{\perp}_{\mathcal{D}'}) \end{array} \right. , $$

and for quadrangles

$$ \left\{ \begin{array}{l} \mathbb{Q}_{k+1}/\mathbb{K} = \ker(\nabla) \\ (\mathbf{d}\mathbb{P}_k(\mathcal{F}) \times \mathbf{d}\mathbb{Q}_k(\mathcal{C}))/k = \text{Range}(\nabla^{\perp}_{\mathcal{D}'}) \end{array} \right. . $$

The location of the degrees of freedom for this discrete de-Rham complex for Cartesian meshes is summarised in Figure 2 and in Figure 4 for triangles.

In this section, results of [25] for discontinuous finite element spaces for vectors have been recalled to ensure the harmonic gap property. Still, as the spaces $dRT$ and $dN$ are obtained by relaxing the normal continuity constraint of the classical conformal finite element spaces $RT$ and $N$, their number of degrees of freedom are not optimal. In the following sections, we will try to develop vector finite element spaces with a lower number of degrees of freedom for which the harmonic gap property holds also, by beginning by the triangular meshes case.

# 4 The triangular mesh case

In this section, we are interested in the following diagram

$$ \mathbb{P}_{k+1} \xrightarrow{\nabla^{\perp}} \mathbf{d}\mathbb{P}_k(\mathcal{C}) \xrightarrow{\nabla^{\perp}_{\mathcal{D}'}} \mathbf{d}\mathbb{P}_k(\mathcal{F}) \times \mathbf{d}\mathbb{P}_{k-1}(\mathcal{C}). \quad (12) $$

## 4.1 Dimension of the finite element spaces

We first compute the dimension of each of the finite element spaces involved in (12), induced by the number of faces, points and cells that was computed in Proposition 1.

**Proposition 6 (Dimension of the finite element spaces).** If the mesh is triangular and periodic, then

$$ \left\{ \begin{array}{l} \dim \mathbb{P}_{k+1} = \frac{N(k+1)^2}{2} \\ \dim \mathbf{d}\mathbb{P}_k(\mathcal{C}) = N(k+1)(k+2) \\ \dim (\mathbf{d}\mathbb{P}_k(\mathcal{F}) \times \mathbf{d}\mathbb{P}_{k-1}(\mathcal{C})) = \frac{(k+1)(k+3)N}{2} \end{array} \right. . $$

*Proof.* The dimension of $\mathbb{P}_{k+1}$ was already proven in Proposition 3. We are now interested in the dimension of $\mathbf{d}\mathbb{P}_k(\mathcal{C})$. This space is a vector space, and so is composed of two components, each of these components having $\frac{(k+1)(k+2)}{2}$ degrees of freedom on each cell. This gives

$$ \dim \mathbf{d}\mathbb{P}_k(\mathcal{C}) = 2N \times \frac{(k+1)(k+2)}{2} = N(k+1)(k+2). $$

We are finally interested in the dimension of $\mathbf{d}\mathbb{P}_k(\mathcal{F}) \times \mathbf{d}\mathbb{P}_{k-1}(\mathcal{C})$. This finite element space includes $k+1$ degrees of freedom on each face, and $\frac{k(k+1)}{2}$ degrees of freedom on each cell.
It remains to compute the dimension of $d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(\mathcal{C})$

$$ \begin{aligned} \dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(\mathcal{C})) &= (k+1)\#\mathcal{F} + (k+1)^2 \#\mathcal{C} \\ &= (k+1)2N + (k+1)^2 N \\ &= N(k+1)(k+3). \end{aligned} $$

<figure><img src="image_2.png" /></figure>

Properties of the complex (11) was addressed in a more general framework in [29], and lead in dimension 2 to the following proposition

**Proposition 5.** The discrete diagram (11) ensures the Proposition 1. Moreover, for triangles:

$$ \left\{ \begin{array}{l} \mathbb{P}_{k+1}/\mathbb{K} = \ker (\nabla^{\perp}) \\ (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(\mathcal{C})) / k = \text{Range}(\nabla_{\mathcal{D}'}^{\perp}), \end{array} \right. $$

and for quadrangles

$$ \left\{ \begin{array}{l} \mathbb{Q}_{k+1}/\mathbb{K} = \ker (\nabla^{\perp}) \\ (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(\mathcal{C})) / k = \text{Range}(\nabla_{\mathcal{D}'}^{\perp}). \end{array} \right. $$

The location of the degrees of freedom for this discrete de-Rham complex for Cartesian meshes is summarised in Figure 1, and in Figure 3 for triangles.

By changing the representation of the linear forms, which is equivalent to rotating of $\pi/2$ the vector spaces, the following proposition is also obtained:

**Proposition 6.** The discrete diagram

$$ \left\{ \begin{array}{cccc} \mathbb{P}_{k+1} & \xrightarrow{\nabla} & \mathbf{d}\mathbb{N}_{k+1}^{\Delta}(\mathcal{C}) & \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} & d\mathbb{P}_k(\mathcal{C}) \times d\mathbb{P}_k(\mathcal{F}) \\ \mathbb{Q}_{k+1} & \xrightarrow{\nabla} & \mathbf{d}\mathbb{N}_{k+1}^{\Box}(\mathcal{C}) & \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} & d\mathbb{Q}_k(\mathcal{C}) \times d\mathbb{P}_k(\mathcal{F}). \end{array} \right. $$

where $\nabla_{\mathcal{D}'}^{\perp}$ is $\nabla^{\perp}$ in the sense of distributions, ensures the Proposition 1. Moreover, for triangles:

$$ \left\{ \begin{array}{l} \mathbb{P}_{k+1}/\mathbb{K} = \ker (\nabla) \\ (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_k(\mathcal{C})) / k = \text{Range}(\nabla_{\mathcal{D}'}^{\perp}), \end{array} \right. $$

and for quadrangles

$$ \left\{ \begin{array}{l} \mathbb{Q}_{k+1}/\mathbb{K} = \ker (\nabla) \\ (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{Q}_k(\mathcal{C})) / k = \text{Range}(\nabla_{\mathcal{D}'}^{\perp}). \end{array} \right. $$

The location of the degrees of freedom for this discrete de-Rham complex for Cartesian meshes is summarised in Figure 2 and in Figure 4 for triangles.

In this section, results of [29] for discontinuous finite element spaces for vectors have been recalled to ensure the Proposition 1. Still, as the spaces $d\mathbb{R}^T$ and $d\mathbb{N}$ are obtained by relaxing the normal continuity constraint of the classical conformal finite element spaces $\mathbb{R}^T$ and $\mathbb{N}$, their number of degrees of freedom are not optimal. In the following sections, we will try to develop vector finite element spaces with a lower number of degrees of freedom for which the Proposition 1 holds also, by beginning by the triangular meshes case.

# 4 The triangular mesh case

In this section, we are interested in the following diagram

$$ \mathbb{P}_{k+1} \xrightarrow{\nabla^{\perp}} \mathbf{d}\mathbb{P}_k(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_{k-1}(\mathcal{C}). \quad (12) $$

## 4.1 Dimension of the finite elements spaces

We first compute the dimension of each of the finite element spaces involved in (12), induced by the number of faces, points and cells that was computed in Proposition 2.
<figure><img src="image_3.png" /><figcaption>Figure 1: Representation of the degrees of freedom of the finite element spaces involved in the curl/div de-Rham complex for Cartesian meshes for $k = 0, 1, 2$ and $3$. Points denote scalar degrees of freedom, whereas arrows denote vectorial degrees of freedom. Both scalar and vectorial volume degrees of freedom are represented in blue, whereas the degrees of freedom in the face finite element space are represented in red.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 1: Representation of the degrees of freedom of the finite element spaces involved in the curl/div de-Rham complex for Cartesian meshes for $k = 0, 1, 2$ and $3$. Points denote scalar degrees of freedom, whereas arrows denote vectorial degrees of freedom. Both scalar and vectorial volume degrees of freedom are represented in blue, whereas the degrees of freedom in the face finite element space are represented in red.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Representation of the degrees of freedom of the finite element spaces involved in the grad/curl de-Rham complex for Cartesian meshes for $k = 0, 1, 2$ and $3$. Same code for colors as in Figure 1 is used.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Representation of the degrees of freedom of the finite element spaces involved in the grad/curl de-Rham complex for Cartesian meshes for $k = 0, 1, 2$ and $3$. Same code for colors as in Figure 1 is used.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Representation of the finite element spaces involved in the curl/div de-Rham complex for triangular meshes for $k = 0, 1, 2$ and $3$. Same code for colors as in Figure 1 is used.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Representation of the finite element spaces involved in the curl/div de-Rham complex for triangular meshes for $k = 0, 1, 2$ and $3$. Same code for colors as in Figure 1 is used.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Representation of the finite element spaces involved in the grad/curl de-Rham complex for triangular meshes for $k = 0, 1, 2$ and $3$. Same code for as in Figure 1 is used.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Representation of the finite element spaces involved in the grad/curl de-Rham complex for triangular meshes for $k = 0, 1, 2$ and $3$. Same code for as in Figure 1 is used.</figcaption></figure>
Adding all these degrees of freedom leads to

$$
\begin{align*}
\dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_{k-1}(\mathcal{C})) &= (k+1)\#\mathcal{F} + \frac{k(k+1)}{2}\#\mathcal{C} \\
&= (k+1)\frac{3N}{2} + \frac{k(k+1)}{2}N \\
\dim (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_{k-1}(\mathcal{C})) &= \frac{(k+1)(k+3)N}{2},
\end{align*}
$$

which ends the proof of this proposition.

<figure><img src="image_7.png" /></figure>

4.2 Study of ∇ ⊥

We are now interested in the study of the ∇ ⊥ operator.

**Proposition 7** ($\nabla^{\perp}$ in the triangular case). We have

$$
\left\{
\begin{array}{l}
\dim (\ker \nabla^{\perp}) = 1 \\
\rank (\nabla^{\perp}) = \frac{N(k+1)^2}{2} - 1.
\end{array}
\right.
$$

Proof. Suppose that a $\psi \in \mathbb{P}_{k+1}$ is such that $\nabla^{\perp}\psi = 0$. Then $\partial_x\psi = \partial_y\psi = 0$, so that $\psi$ is piecewise constant. But as $\psi$ is continuous, it is actually uniform:

$$
\ker \nabla^{\perp} = \mathbb{K}.
$$

This gives dim ker ∇ ⊥ = 1, as we are working on a domain with a single connected component.
Applying the rank-nullity theorem in the triangular case

$$
\dim \mathbb{P}_k = \dim \ker \nabla^{\perp} + \mathrm{rank} \left( \nabla^{\perp} \right),
$$

leads to

$$
\mathrm{rank} \left( \nabla^{\perp} \right) = \dim \mathbb{P}_k - 1 = \frac{N(k+1)^2}{2} - 1.
$$

In the same manner, we get, in the quadrangular case:

$$
\mathrm{rank} \left( \nabla^{\perp} \right) = N(k+1)^2 - 1.
$$

<figure><img src="image_8.png" /></figure>

4.3 Discrete divergence free polynomials on the reference cell

We consider the following application

$$
\mathbf{C}_k^{\partial} : \psi \in d\mathbb{P}_{k+1}(\hat{K}) \mapsto \mathrm{Tr}(\nabla^{\perp}\psi) \in \mathcal{R}_k(\partial\hat{K}) \quad (13)
$$

**Proposition 8.** We denote by $k$ the constant elements of $R_k(\partial \hat{K})$. Then

$$
R_k(\partial \hat{K}) = \mathrm{Range} \, \mathbf{C}_k^{\partial} \oplus \mathbf{k},
$$

where the sum is orthogonal.

Proof. We denote by $c$ an element of $k$. We denote also by $c$ the function equal to $c$ on $\hat{K}$.
We also denote by $u$ an element of $\mathrm{Range} \, \mathbf{C}^{\partial}$. Then a $\psi$ exists such that $u = \mathbf{C}^{\partial}_k(\psi)$. Then

$$
\begin{align*}
\int_{\partial \hat{K}} u c &= \int_{\partial \hat{K}} \mathbf{Tr} (\nabla^{\perp} \psi) c \\
&= \int_{\hat{K}} \nabla \cdot (c \nabla^{\perp} \psi) \\
&= \int_{\hat{K}} c \nabla \cdot (\nabla^{\perp} \psi) + \int_{\hat{K}} \nabla c \cdot \nabla^{\perp} \psi \\
\int_{\partial \hat{K}} u c &= 0,
\end{align*}
$$
**Proposition 7** (Dimension of the finite element spaces). If the mesh is triangular and periodic, then

$$ egin{cases} \dim \mathbb{P}_{k+1} = \frac{N(k+1)^2}{2} \\ \dim \mathbf{d}\mathbb{P}_k(C) = N(k+1)(k+2) \\ \dim (\mathrm{d}\mathbb{P}_k(\mathcal{F}) \times \mathrm{d}\mathbb{P}_{k-1}(C)) = \frac{(k+1)(k+3)N}{2}. \end{cases} $$

*Proof.* The dimension of $\mathbb{P}_{k+1}$ was already proven in Proposition 4. We are now interested in the dimension of $\mathbf{d}\mathbb{P}_k(C)$. This space is a vector space, and so is composed of two components, each of these components having $\frac{(k+1)(k+2)}{2}$ degrees of freedom on each cell. This gives

$$ \dim \mathbf{d}\mathbb{P}_k(C) = 2N \times \frac{(k+1)(k+2)}{2} = N(k+1)(k+2). $$

We are finally interested in the dimension of $\mathrm{d}\mathbb{P}_k(\mathcal{F}) \times \mathrm{d}\mathbb{P}_{k-1}(C)$. This finite element space includes $k+1$ degrees of freedom on each face, and $\frac{k(k+1)}{2}$ degrees of freedom on each cell. Adding all these degrees of freedom leads to

$$ \begin{aligned} \dim (\mathrm{d}\mathbb{P}_k(\mathcal{F}) \times \mathrm{d}\mathbb{P}_{k-1}(C)) &= (k+1)\#\mathcal{F} + \frac{k(k+1)}{2}\#C \\ &= (k+1)\frac{3N}{2} + \frac{k(k+1)}{2}N \\ \dim (\mathrm{d}\mathbb{P}_k(\mathcal{F}) \times \mathrm{d}\mathbb{P}_{k-1}(C)) &= \frac{(k+1)(k+3)N}{2}, \end{aligned} $$

which ends the proof of this proposition.

<figure><img src="image_7.png" /></figure>

## 4.2 Study of $\nabla^\perp$

We are now interested in the study of the $\nabla^\perp$ operator.

**Proposition 8** ($\nabla^\perp$ in the triangular case). We have

$$ \begin{cases} \dim (\ker \nabla^\perp) = 1 \\ \rank (\nabla^\perp) = \frac{N(k+1)^2}{2} - 1. \end{cases} $$

*Proof.* Suppose that a $\psi \in \mathbb{P}_{k+1}$ is such that $\nabla^\perp \psi = 0$. Then $\partial_x \psi = \partial_y \psi = 0$, so that $\psi$ is piecewise constant. But as $\psi$ is continuous, it is actually uniform:

$$ \ker \nabla^\perp = \mathbb{K}. $$

This gives $\dim \ker \nabla^\perp = 1$, as we are working on a domain with a single connected component. Applying the rank-nullity theorem in the triangular case

$$ \dim \mathbb{P}_k = \dim \ker \nabla^\perp + \rank (\nabla^\perp), $$

leads to

$$ \rank (\nabla^\perp) = \dim \mathbb{P}_k - 1 = \frac{N(k+1)^2}{2} - 1. $$

In the same manner, we get, in the quadrangular case:

$$ \rank (\nabla^\perp) = N(k+1)^2 - 1. $$

<figure><img src="image_8.png" /></figure>
because $c$ is constant and $\nabla \cdot (\nabla^{\perp} \psi) = 0$. We thus have proven that the sum

$$ \text{Range} \mathbf{C}_k^{\partial} + k, $$

is direct and orthogonal.

We are now interested in the study of the kernel of $\mathbf{C}_k^{\partial}$. Suppose that an element $\psi$ is such that $\mathbf{C}_k^{\partial}(\psi) = 0$. We consider the classical Lagrange basis of $d\mathbb{P}_{k+1}(\hat{K})$. Then $\psi$ is such that its value on the boundary of $\hat{K}$ is constant, and may take any value on the degrees of freedom matching with the interior nodes. This means that

$$ \dim (\ker \mathbf{C}_k^{\partial}) = 1 + \frac{k(k-1)}{2}. $$

Using the rank-nullity theorem gives

$$ \begin{aligned} \operatorname{rank} (\mathbf{C}_k^{\partial}) &= \dim d\mathbb{P}_{k+1} - \left(1 + \frac{k(k-1)}{2}\right) \\ &= \frac{(k+2)(k+3)}{2} - 1 - \frac{k(k-1)}{2} \\ &= \frac{k^2 + 5k + 6 - k^2 + k - 2}{2} \\ &= \frac{6k + 4}{2} \end{aligned} $$

$$ \operatorname{rank} (\mathbf{C}_k^{\partial}) = 3k + 2. $$

We also know that $\dim R_k(\partial \hat{K}) = 3(k+1)$. We have then $\text{Range} \mathbf{C}_k^{\partial} \oplus k \subset R_k(\partial \hat{K})$, and $\dim (\text{Range} \mathbf{C}_k^{\partial} \oplus k) = \dim R_k(\partial \hat{K})$, so that $\text{Range} \mathbf{C}_k^{\partial} \oplus k = R_k(\partial \hat{K})$, which ends the proof. $\square$

**Proposition 9** (Decomposition of divergence free elements). We denote by $\mathcal{L}^{f,i}$ the Legendre polynomial of degree $i$ on the face $f$ of $\hat{K}$, normalised such that

$$ \int_{\partial \hat{K}} (\mathcal{L}^{f,i})^2 = 1. $$

Suppose that $\mathbf{u} \in d\mathbb{P}_k(\hat{K})$ is divergence free. Then $\mathbf{u}$ can be uniquely decomposed as

$$ \mathbf{u} = \bar{\mathbf{v}}_{\mathbf{u}}^{0} + \bar{\mathbf{v}}_{\mathbf{u}}^{1} + \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \bar{\mathbf{v}}_{\mathbf{u}}^{f,i}, \quad (14) $$

where

* $\mathbf{v}_{\mathbf{u}}^0$ is in the set

$$ \Phi_k = \{ \mathbf{u} \in d\mathbb{P}_k \quad \nabla \cdot \mathbf{u} = 0 \quad \mathrm{Tr}(\mathbf{u}) = 0 \}. $$

* $\mathbf{v}_{\mathbf{u}}^1$ is constant.

* $\bar{\mathbf{v}}_{\mathbf{u}}^{f,i}$ is orthogonal to $\Phi_k$ and such that

$$ - \nabla \cdot \bar{\mathbf{v}}_{\mathbf{u}}^{f,i} = 0. $$

$$ - \mathrm{Tr}(\bar{\mathbf{v}}_{\mathbf{u}}^{f,i}) \text{ is orthogonal to all } \mathcal{L}^{g,j} \text{ for } \{g, j\} \neq \{f, i\}. $$

Note that the sum over the integers $i$ in the sum of $\bar{\mathbf{v}}_{\mathbf{u}}^{f,i}$ in (14) begins at 1 and not 0 for excluding the ones that would have a constant trace. Note also that the set $\Phi_k$ is the same as in the decomposition used in [13, Proposition 3.1] for the derivation of classical conformal divergence free finite elements of Raviart-Thomas [27, 28] or Brezzi-Douglas-Marini [12] types (these are respectively referred as RT and BDM in [6]).
## 4.3 Discrete divergence free polynomials on the reference cell

We consider the following application

$$ \mathbf{C}_k^{\partial}: \psi \in \mathrm{d}\mathbb{P}_{k+1}(\hat{K}) \mapsto \mathrm{Tr}(\nabla^{\perp} \psi) \in \mathrm{d}\mathbb{P}_k(\partial\hat{K}), \quad (13) $$

where $\hat{K}$ is the reference triangle.

**Proposition 9.** We denote by $k$ the constant elements of $\mathrm{d}\mathbb{P}_k(\partial\hat{K})$. Then

$$ \mathrm{d}\mathbb{P}_k(\partial\hat{K}) = \mathrm{Range} \, \mathbf{C}_k^{\partial} \oplus k, $$

where the sum is orthogonal.

*Proof.* We denote by $c$ an element of $k$. We denote also by $c$ the function equal to $c$ on $\hat{K}$. We also denote by $u$ an element of Range $\mathbf{C}^{\partial}$. Then a $\psi$ exists such that $u = \mathbf{C}_k^{\partial}(\psi)$. Then

$$ \begin{aligned} \int_{\partial \hat{K}} uc &= \int_{\partial \hat{K}} \mathrm{Tr}(\nabla^{\perp} \psi) c \\ &= \int_{\hat{K}} \nabla \cdot (c \nabla^{\perp} \psi) \\ &= \int_{\hat{K}} c \nabla \cdot (\nabla^{\perp} \psi) + \int_{\hat{K}} \nabla c \cdot \nabla^{\perp} \psi \\ \int_{\partial \hat{K}} uc &= 0, \end{aligned} $$

because $c$ is constant and $\nabla \cdot (\nabla^{\perp} \psi) = 0$. We thus have proven that the sum

$$ \mathrm{Range} \, \mathbf{C}_k^{\partial} + k, $$

is direct and orthogonal.

We are now interested in the study of the kernel of $\mathbf{C}_k^{\partial}$. Suppose that an element $\psi$ is such that $\mathbf{C}_k^{\partial}(\psi) = 0$. We consider the classical Lagrange basis of $\mathrm{d}\mathbb{P}_{k+1}(\hat{K})$. Then $\psi$ is such that its value on the boundary of $\hat{K}$ is constant, and may take any value on the degrees of freedom matching with the interior nodes. This means that

$$ \dim(\ker \mathbf{C}_k^{\partial}) = 1 + \frac{k(k-1)}{2}. $$

Using the rank-nullity theorem gives

$$ \begin{aligned} \mathrm{rank} (\mathbf{C}_k^{\partial}) &= \dim \mathrm{d}\mathbb{P}_{k+1} - \left(1 + \frac{k(k-1)}{2}\right) \\ &= \frac{(k+2)(k+3)}{2} - 1 - \frac{k(k-1)}{2} \\ &= \frac{k^2 + 5k + 6 - k^2 + k - 2}{2} \\ &= \frac{6k + 4}{2} \end{aligned} $$

$$ \mathrm{rank} (\mathbf{C}_k^{\partial}) = 3k + 2. $$

We also know that $\dim \mathrm{d}\mathbb{P}_k(\partial\hat{K}) = 3(k+1)$. We have then $\mathrm{Range} \, \mathbf{C}_k^{\partial} \oplus k \subset \mathrm{d}\mathbb{P}_k(\partial\hat{K})$, and $\dim(\mathrm{Range} \, \mathbf{C}_k^{\partial} \oplus k) = \dim \mathrm{d}\mathbb{P}_k(\partial\hat{K})$, so that $\mathrm{Range} \, \mathbf{C}_k^{\partial} \oplus k = \mathrm{d}\mathbb{P}_k(\partial\hat{K})$, which ends the proof. $\square$

**Proposition 10** (Decomposition of divergence free elements). We denote by $\mathcal{L}^{f,i}$ the Legendre polynomial of degree $i$ on the face $f$ of $\hat{K}$, normalised such that

$$ \int_{\partial \hat{K}} (\mathcal{L}^{f,i})^2 = 1. $$

Suppose that $\mathbf{u} \in \mathrm{d}\mathbb{P}_k(\hat{K})$ is divergence free. Then $\mathbf{u}$ can be uniquely decomposed as

$$ \mathbf{u} = \bar{\mathbf{v}}_{\mathbf{u}}^{0} + \bar{\mathbf{v}}_{\mathbf{u}}^{1} + \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \bar{\mathbf{v}}_{\mathbf{u}}^{f,i}, \quad (14) $$

where
Proof. We first prove that $\bar{\mathbf{v}}_u^1$, if it exists, is orthogonal to $\Phi_k$. We denote by $\mathbf{u}_\Phi$ an element of $\Phi_k$. Then a $\psi^\Phi$ exists such that $\mathbf{u}_\Phi = \nabla^\perp \Psi^\Phi$. As $\mathbf{u}_\Phi \in \Phi_k$, $\Psi^\Phi$ is constant on $\partial \hat{K}$. Then

$$
\begin{align*}
\int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \mathbf{u}_{\Phi} &= \int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \nabla^{\perp} \Psi^{\Phi} \\
&= \int_{\hat{K}} (\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \cdot \nabla \Psi^{\Phi} \\
&= \int_{\hat{K}} \nabla \cdot ( (\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \Psi^{\Phi} ) \\
&= \int_{\partial \hat{K}} \mathrm{Tr} \left( (\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \Psi^{\Phi} \right) \\
\int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \mathbf{u}_{\Phi} &= 0,
\end{align*}
$$

because both $(\bar{\mathbf{v}}_u^1)^\perp$ and $\Psi^\Phi$ are constant on $\partial \hat{K}$. We have then proven that $\bar{\mathbf{v}}_u^1 \perp \Phi_k$.
We can now prove the uniqueness of the decomposition. We suppose that

$$0 = \bar{\mathbf{v}}_{\mathbf{u}}^{0} + \bar{\mathbf{v}}_{\mathbf{u}}^{1} + \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \bar{\mathbf{v}}_{\mathbf{u}}^{f, i},$$

where the different components ensure the properties of the proposition. By definition, $\bar{\mathbf{v}}_u^0$ is orthogonal to the $\bar{\mathbf{v}}_u^{f,i}$, and we proved that it is also orthogonal to $\bar{\mathbf{v}}_u^1$, so that it vanishes. We take the trace of the remaining part, which gives

$$0 = \mathrm{Tr}(\bar{\mathbf{v}}_{\mathbf{u}}^{1}) + \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \mathrm{Tr}(\bar{\mathbf{v}}_{\mathbf{u}}^{f, i}).$$

Taking the scalar product by any $\mathcal{L}^{g,j}$ for $j \ge 1$ gives

$$\int_{\partial \hat{K}} \mathrm{Tr} \left( \bar{\mathbf{v}}_{\mathbf{u}}^{g,j} \right) \mathcal{L}^{g,j} = 0.$$

As $\mathrm{Tr}(\bar{\mathbf{v}}_u^{g,j})$ has a moment only on $\{g,j\}$, it is actually zero. As it is orthogonal to $\Phi_k$, we get $\bar{\mathbf{v}}_u^{g,j} = 0$ for all $g,j$. It remains then $\bar{\mathbf{v}}_u^1 = 0$, and so each component is zero, which proves the uniqueness of the decomposition.

We now prove the existence. We consider one $\mathcal{L}^{f,i}$ of $R_k(\partial \hat{K})$, for $i \ge 1$. $\mathcal{L}^{f,i}$ is orthogonal to $\mathbf{k}$, and so using using Proposition 8, $\mathcal{L}^{f,i}$ is in Range $C_k^\partial$, so that a $\psi^{f,i}$ exists such that $C_k^\partial(\psi^{f,i}) = \mathcal{L}^{f,i}$. Denoting by $\mathcal{P}$ the orthogonal projection on $\Phi_k$, we define $\mathbf{e}^{f,i} := \nabla^\perp (\psi^{f,i}) - \mathcal{P}(\nabla^\perp (\psi^{f,i}))$. Then $\mathbf{e}^{f,i}$ is orthogonal to $\Phi_k$, is divergence free, and is such that $\mathrm{Tr}(\mathbf{e}^{f,i}) = \mathcal{L}^{f,i}$. We consider one divergence free $\mathbf{u} \in d\mathbb{P}_k(\hat{K})$. We define

$$\lambda_{f,i} := \int_{\partial \hat{K}} \mathrm{Tr}(\mathbf{u}) \mathcal{L}^{f,i},$$

and set $\bar{\mathbf{v}}_u^{f,i} = \lambda_{f,i} \mathbf{e}^{f,i}$. We also set $\bar{\mathbf{v}}_u^0 = \mathcal{P}(\mathbf{u})$, and

$$\bar{\mathbf{v}}_{\mathbf{u}}^{1} := \mathbf{u} - \bar{\mathbf{v}}_{\mathbf{u}}^{0} - \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \bar{\mathbf{v}}_{\mathbf{u}}^{f, i}.$$

Then $\bar{\mathbf{v}}_u^0$ and the $\bar{\mathbf{v}}_u^{f,i}$ ensure all the properties required. It remains to prove that $\bar{\mathbf{v}}_u^1$ is constant. We know that $\bar{\mathbf{v}}_u^1$ is divergence free, orthogonal to $\Phi_k$, and that its trace is constant on each face, because all the components in $\mathcal{L}^{f,i}$ for $i \ge 1$ were removed. We denote by $\mathbf{k}_y$ the opposite of the trace on the side linking $[0,0]$ to $[0,1]$, and $\mathbf{k}_x$ the opposite of the trace on the side linking $[0,0]$ to $[1,0]$, and consider $\mathbf{k} = (\mathbf{k}_x, \mathbf{k}_y)$. As $\nabla \cdot \bar{\mathbf{v}}_u^1 = 0$, we have

$$
\begin{align*}
0 &= \int_{\hat{K}} \nabla \cdot \bar{\mathbf{v}}_{\mathbf{u}}^{1} \\
&= \int_{\partial \hat{K}} \mathrm{Tr} (\bar{\mathbf{v}}_{\mathbf{u}}^{1}) \\
&= \int_{(0,0)}^{(1,0)} \mathrm{Tr} (\bar{\mathbf{v}}_{\mathbf{u}}^{1}) + \int_{(1,0)}^{(1,1)} \mathrm{Tr} (\bar{\mathbf{v}}_{\mathbf{u}}^{1}) + \int_{(1,1)}^{(0,0)} \mathrm{Tr} (\bar{\mathbf{v}}_{\mathbf{u}}^{1}), \\
&= -\mathbf{k}_x - \mathbf{k}_y + \int_{(1,0)}^{(1,1)} \mathrm{Tr} (\bar{\mathbf{v}}_{\mathbf{u}}^{1}) \\
0 &= -\mathbf{k}_x - \mathbf{k}_y + \sqrt{2} \mathrm{Tr} (\bar{\mathbf{v}}_{\mathbf{u}}^{1}) |_{[(1,0),(1,1)]}
\end{align*}
$$
• $\bar{\mathbf{v}}_u^0$ is in the set

$$ \Phi_k = \{ \mathbf{u} \in d\mathbb{P}_k \quad \nabla \cdot \mathbf{u} = 0 \quad \mathrm{Tr}(\mathbf{u}) = 0 \}. $$

• $\bar{\mathbf{v}}_u^1$ is constant.

• $\bar{\mathbf{v}}_u^{f,i}$ is orthogonal to $\Phi_k$ and such that

$$ - \nabla \cdot \bar{\mathbf{v}}_u^{f,i} = 0. $$

$$ - \mathrm{Tr}(\bar{\mathbf{v}}_u^{f,i}) \text{ is orthogonal to all } \mathcal{L}^{g,j} \text{ for } \{g,j\} \neq \{f,i\}. $$

Note that the sum over the integers $i$ in the sum of $\bar{\mathbf{v}}_u^{f,i}$ in (14) begins at 1 and not 0 for excluding the ones that would have a constant trace. Note also that the set $\Phi_k$ is the same as in the decomposition used in [16, Proposition 3.1] for the derivation of classical conformal divergence free finite elements of Raviart-Thomas [32, 33] or Brezzi-Douglas-Marini [15] types (these are respectively referred as RT and BDM in [6]).

*Proof.* We first prove that $\bar{\mathbf{v}}_u^1$, if it exists, is orthogonal to $\Phi_k$. We denote by $\mathbf{u}_\Phi$ an element of $\Phi_k$. Then a $\psi^\Phi$ exists such that $\mathbf{u}_\Phi = \nabla^\perp \Psi^\Phi$. As $\mathbf{u}_\Phi \in \Phi_k$, $\Psi^\Phi$ is constant on $\partial \hat{K}$. Then

$$ \begin{aligned} \int_{\hat{K}} \bar{\mathbf{v}}_u^1 \cdot \mathbf{u}_\Phi &= \int_{\hat{K}} \bar{\mathbf{v}}_u^1 \cdot \nabla^\perp \Psi^\Phi \\ &= \int_{\hat{K}} (\bar{\mathbf{v}}_u^1)^\perp \cdot \nabla \Psi^\Phi \\ &= \int_{\hat{K}} \nabla \cdot ( (\bar{\mathbf{v}}_u^1)^\perp \Psi^\Phi ) \\ &= \int_{\partial \hat{K}} \mathrm{Tr} ( (\bar{\mathbf{v}}_u^1)^\perp \Psi^\Phi ) \\ \int_{\hat{K}} \bar{\mathbf{v}}_u^1 \cdot \mathbf{u}_\Phi &= 0, \end{aligned} $$

because both $(\bar{\mathbf{v}}_u^1)^\perp$ and $\Psi^\Phi$ are constant on $\partial \hat{K}$. We have then proven that $\bar{\mathbf{v}}_u^1 \perp \Phi_k$.

We can now prove the uniqueness of the decomposition. We suppose that

$$ 0 = \bar{\mathbf{v}}_u^0 + \bar{\mathbf{v}}_u^1 + \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^k \bar{\mathbf{v}}_u^{f,i}, $$

where the different components ensure the properties of the proposition. By definition, $\bar{\mathbf{v}}_u^0$ is orthogonal to the $\bar{\mathbf{v}}_u^{f,i}$, and we proved that it is also orthogonal to $\bar{\mathbf{v}}_u^1$, so that it vanishes. We take the trace of the remaining part, which gives

$$ 0 = \mathrm{Tr}(\bar{\mathbf{v}}_u^1) + \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^k \mathrm{Tr}(\bar{\mathbf{v}}_u^{f,i}). $$

Taking the scalar product by any $\mathcal{L}^{g,j}$ for $j \ge 1$ gives

$$ \int_{\partial \hat{K}} \mathrm{Tr}(\bar{\mathbf{v}}_u^{g,j}) \mathcal{L}^{g,j} = 0. $$

As $\mathrm{Tr}(\bar{\mathbf{v}}_u^{g,j})$ has a moment only on face $g$ for degree $j$, it is actually zero. As it is orthogonal to $\Phi_k$, we get $\bar{\mathbf{v}}_u^{g,j} = 0$ for all $g, j$. It remains then $\bar{\mathbf{v}}_u^1 = 0$, and so each component is zero, which proves the uniqueness of the decomposition.

We now prove the existence. We consider one $\mathcal{L}^{f,i}$ of $d\mathbb{P}_k(\partial \hat{K})$, for $i \ge 1$. $\mathcal{L}^{f,i}$ is orthogonal to $\mathcal{K}$, and so using using Proposition 9, $\mathcal{L}^{f,i}$ is in Range $C_k^\partial$, so that a $\psi^{f,i}$ exists such that $C_k^\partial(\psi^{f,i}) = \mathcal{L}^{f,i}$. Denoting by $\mathcal{P}$ the orthogonal projection on $\Phi_k$, we define $\mathbf{e}^{f,i} := \nabla^\perp (\psi^{f,i}) - \mathcal{P}(\nabla^\perp (\psi^{f,i}))$. Then $\mathbf{e}^{f,i}$ is orthogonal to $\Phi_k$, is divergence free, and is such that $\mathrm{Tr}(\mathbf{e}^{f,i}) = \mathcal{L}^{f,i}$. We consider one divergence free $\mathbf{u} \in d\mathbb{P}_k(\hat{K})$. We define

$$ \lambda_{f,i} := \int_{\partial \hat{K}} \mathrm{Tr}(\mathbf{u}) \mathcal{L}^{f,i}, $$

and set $\bar{\mathbf{v}}_u^{f,i} = \lambda_{f,i} \mathbf{e}^{f,i}$. We also set $\bar{\mathbf{v}}_u^0 = \mathcal{P}(\mathbf{u})$, and

$$ \bar{\mathbf{v}}_u^1 := \mathbf{u} - \bar{\mathbf{v}}_u^0 - \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^k \bar{\mathbf{v}}_u^{f,i}. $$
which means that the trace on the side linking [1, 0] to [1, 1] is equal to $\frac{1}{\sqrt{2}}(\mathbf{k}_x + \mathbf{k}_y)$, which is also equal to the trace of $\mathbf{k}$. This means $\bar{\mathbf{v}}_u^1 - \mathbf{k}$ is divergence free, and that $\mathrm{Tr}(\bar{\mathbf{v}}_u^1 - \mathbf{k}) = 0$, and so $\bar{\mathbf{v}}_u^1 - \mathbf{k} \in \Phi_k$.

As $\mathbf{k}$ is orthogonal to $\Phi_k$, and so is $\bar{\mathbf{v}}_u^1$, we conclude that $\bar{\mathbf{v}}_u^1 - \mathbf{k}$ is also orthogonal to $\Phi_k$. This gives $\bar{\mathbf{v}}_u^1 - \mathbf{k} = 0$, and $\bar{\mathbf{v}}_u^1$ is therefore constant. $\square$

It is important to note that the decomposition of Proposition 9 was proven on the reference element. However, all the properties of the different spaces are invariant by linear transformations. This means that the decomposition of Proposition 9 holds actually on any straight triangular cell of the mesh.

## 4.4 $\ker(\nabla_{\mathcal{G}'\cdot})$ and $\mathrm{Range}(\nabla^{\perp})$

**Proposition 10.** If $\mathbb{K}$ denotes the set of uniform vectors, then

$$ \ker(\nabla_{\mathcal{G}'\cdot}) = \mathrm{Range}(\nabla^{\perp}) \oplus \mathbb{K}. $$

*Proof.* We begin by proving that $\mathrm{Range}(\nabla^{\perp}) \subset \ker(\nabla_{\mathcal{G}'\cdot})$. We consider an element $\mathbf{u}$ of $\mathrm{Range}(\nabla^{\perp})$. Then a $\Psi \in \mathbb{P}_{k+1}$ exists such that $\mathbf{u} = \nabla^{\perp}\psi$. Then

$$ \forall c \in \mathcal{C} \quad \nabla \cdot (\nabla^{\perp} \psi) = 0. $$

Also, as $\psi$ is continuous, the jump of $\nabla^{\perp}\psi$ across faces vanishes. We have then proven that $\mathrm{Range}(\nabla^{\perp}) \subset \ker(\nabla_{\mathcal{G}'\cdot})$.

We now prove that $\mathrm{Range}(\nabla^{\perp}) \perp \mathbb{K}$. We denote by $\mathbf{k} = (\mathbf{k}_x, \mathbf{k}_y)^T \in \mathbb{K}$, and by $\psi$ an element of $\mathbb{P}_k$. Then

$$ \mathbf{k} \cdot \nabla^{\perp} \psi = \begin{pmatrix} \mathbf{k}_x \\ \mathbf{k}_y \end{pmatrix} \cdot \begin{pmatrix} -\partial_y \psi \\ \partial_x \psi \end{pmatrix} = -\mathbf{k}_x \partial_y \psi + \mathbf{k}_y \partial_x \psi = \begin{pmatrix} \mathbf{k}_y \\ -\mathbf{k}_x \end{pmatrix} \cdot \nabla \psi. $$

We denote by $\mathbf{k}^{\perp} = (\mathbf{k}_y, -\mathbf{k}_x)^T$. Then as $\mathbf{k}^{\perp}$ is uniform, we have on all cells:

$$ \nabla \cdot (\mathbf{k}^{\perp} \psi) = \mathbf{k}^{\perp} \cdot \nabla \psi. $$

Then

$$ \begin{aligned} \sum_{c \in \mathcal{C}} \int_c \mathbf{k}^{\perp} \cdot \nabla \psi &= \sum_{c \in \mathcal{C}} \int_c \nabla \cdot (\mathbf{k}^{\perp} \psi) \\ &= \sum_{c \in \mathcal{C}} \int_c \psi \mathbf{k}^{\perp} \cdot \mathbf{n}_{out} \\ &= \sum_{c \in \mathcal{C}} \sum_{f \in \mathcal{F}(c)} \int_f \psi \mathbf{k}^{\perp} \cdot \mathbf{n}_{out} \\ &= -\sum_f \int_f [\psi \mathbf{k}^{\perp} \cdot \mathbf{n}_f] \\ \sum_{c \in \mathcal{C}} \int_c \mathbf{k}^{\perp} \cdot \nabla \psi &= 0, \end{aligned} $$

because both $\mathbf{k}^{\perp}$ and $\psi$ are continuous across the faces. We have then proven that

$$ \mathbb{K} \perp \mathrm{Range}(\nabla^{\perp}). $$

We also remark that $\mathbb{K} \subset \ker(\nabla_{\mathcal{G}'\cdot})$, because $\nabla_{\mathcal{G}'\cdot}$ is a derivation operator. For the moment, we have proven that

$$ \mathrm{Range}(\nabla^{\perp}) \oplus \mathbb{K} \subset \ker \nabla_{\mathcal{G}'\cdot}. $$

Suppose now that an element $\mathbf{u} \in \mathbf{d}\mathbb{P}_k$ is such that its divergence is 0, namely

$$ \left\{ \begin{array}{ll} \forall c \in \mathcal{C} & \nabla \cdot \mathbf{u} = 0 \\ \forall f \in \mathcal{F} & [\mathbf{u} \cdot \mathbf{n}_f] = 0. \end{array} \right. $$

As $\nabla \cdot \mathbf{u} = 0$ on all the cells, $\mathbf{u}$ can be decomposed as in Proposition 9; this decomposition involves three types of components:
Then $\bar{\mathbf{v}}_u^0$ and the $\bar{\mathbf{v}}_u^{f,i}$ ensure all the properties required. It remains to prove that $\bar{\mathbf{v}}_u^1$ is constant. We know that $\bar{\mathbf{v}}_u^1$ is divergence free, orthogonal to $\Phi_k$, and that its trace is constant on each face, because all the components in $\mathcal{L}^{f,i}$ for $i \ge 1$ were removed. We denote by $\mathbf{k}_y$ the opposite of the trace on the side linking $[0,0]$ to $[0,1]$, and $\mathbf{k}_x$ the opposite of the trace on the side linking $[0,0]$ to $[1,0]$, and consider $\mathbf{k} = (\mathbf{k}_x, \mathbf{k}_y)$. As $\nabla \cdot \bar{\mathbf{v}}_u^1 = 0$, we have

$$
\begin{align*}
0 &= \int_{\bar{K}} \nabla \cdot \bar{\mathbf{v}}_u^1 \\
&= \int_{\partial \bar{K}} \mathbf{Tr}(\bar{\mathbf{v}}_u^1) \\
&= \int_{(0,0)}^{(1,0)} \mathbf{Tr}(\bar{\mathbf{v}}_u^1) + \int_{(1,0)}^{(1,1)} \mathbf{Tr}(\bar{\mathbf{v}}_u^1) + \int_{(1,1)}^{(0,0)} \mathbf{Tr}(\bar{\mathbf{v}}_u^1), \\
&= -\mathbf{k}_x - \mathbf{k}_y + \int_{(1,0)}^{(1,1)} \mathbf{Tr}(\bar{\mathbf{v}}_u^1) \\
0 &= -\mathbf{k}_x - \mathbf{k}_y + \sqrt{2} \mathbf{Tr}(\bar{\mathbf{v}}_u^1)_{|[(1,0),(1,1)]}
\end{align*}
$$

which means that the trace on the side linking $[1,0]$ to $[1,1]$ is equal to $\frac{1}{\sqrt{2}}(\mathbf{k}_x + \mathbf{k}_y)$, which is also equal to the trace of $\mathbf{k}$. This means $\bar{\mathbf{v}}_u^1 - \mathbf{k}$ is divergence free, and that $\mathbf{Tr}(\bar{\mathbf{v}}_u^1 - \mathbf{k}) = 0$, and so $\bar{\mathbf{v}}_u^1 - \mathbf{k} \in \Phi_k$.

As $\mathbf{k}$ is orthogonal to $\Phi_k$, and so is $\bar{\mathbf{v}}_u^1$, we conclude that $\bar{\mathbf{v}}_u^1 - \mathbf{k}$ is also orthogonal to $\Phi_k$. This gives $\bar{\mathbf{v}}_u^1 - \mathbf{k} = 0$, and $\bar{\mathbf{v}}_u^1$ is therefore constant. $\square$

It is important to note that the decomposition of Proposition 10 was proven on the reference element. However, all the properties of the different spaces are invariant by linear transformations. This means that the decomposition of Proposition 10 holds actually on any straight triangular cell of the mesh.

## 4.4 ker $(\nabla_{\mathcal{G}'})$ and Range $(\nabla^{\perp})$

**Proposition 11.** If $\mathbb{K}$ denotes the set of uniform vectors, then

$$ \ker (\nabla_{\mathcal{G}'}) = \text{Range} (\nabla^{\perp}) \oplus \mathbb{K}. $$

*Proof.* We begin by proving that Range $(\nabla^{\perp}) \subset \ker (\nabla_{\mathcal{G}'})$. We consider an element $\mathbf{u}$ of Range $(\nabla^{\perp})$. Then a $\Psi \in \mathbb{P}_{k+1}$ exists such that $\mathbf{u} = \nabla^{\perp}\psi$. Then

$$ \forall c \in \mathcal{C} \quad \nabla \cdot (\nabla^{\perp} \psi) = 0. $$

Also, as $\psi$ is continuous, the jump of $\nabla^{\perp}\psi$ across faces vanishes. We have then proven that Range $(\nabla^{\perp}) \subset \ker (\nabla_{\mathcal{G}'})$.

We now prove that Range $(\nabla^{\perp}) \perp \mathbb{K}$. We denote by $\mathbf{k} = (\mathbf{k}_x, \mathbf{k}_y)^T \in \mathbb{K}$, and by $\psi$ an element of $\mathbb{P}_k$. Then

$$ \mathbf{k} \cdot \nabla^{\perp} \psi = \begin{pmatrix} \mathbf{k}_x \\ \mathbf{k}_y \end{pmatrix} \cdot \begin{pmatrix} -\partial_y \psi \\ \partial_x \psi \end{pmatrix} = -\mathbf{k}_x \partial_y \psi + \mathbf{k}_y \partial_x \psi = \begin{pmatrix} \mathbf{k}_y \\ -\mathbf{k}_x \end{pmatrix} \cdot \nabla \psi. $$

We denote by $\mathbf{k}^{\perp} = (\mathbf{k}_y, -\mathbf{k}_x)^T$. Then as $\mathbf{k}^{\perp}$ is uniform, we have on all cells:

$$ \nabla \cdot (\mathbf{k}^{\perp} \psi) = \mathbf{k}^{\perp} \cdot \nabla \psi. $$

Then

$$
\begin{align*}
\sum_{c \in \mathcal{C}} \int_c \mathbf{k}^{\perp} \cdot \nabla \psi &= \sum_{c \in \mathcal{C}} \int_c \nabla \cdot (\mathbf{k}^{\perp} \psi) \\
&= \sum_{c \in \mathcal{C}} \int_{\partial c} \psi \mathbf{k}^{\perp} \cdot \mathbf{n}_{out} \\
&= \sum_{c \in \mathcal{C}} \sum_{f \in \mathcal{F}(c)} \int_f \psi \mathbf{k}^{\perp} \cdot \mathbf{n}_{out} \\
&= -\sum_f \int_f [\psi \mathbf{k}^{\perp} \cdot \mathbf{n}_f] \\
&\qquad \sum_{c \in \mathcal{C}} \int_c \mathbf{k}^{\perp} \cdot \nabla \psi = 0,
\end{align*}
$$
• The ones of $\Phi_k$, which have a trace equal to 0. This set is of dimension $\frac{k(k-1)}{2}$ on each cell.

• The constant component; this set is of dimension 2 on each cell.

• The components $\bar{v}_i^f$ for $1 \le i \le k$ and for all faces. This set is of dimension 3k on each cell.

Let us see now what is the effect of the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ on these different components. We first remark that the traces of the different components are orthogonal two at a time, which means that we can consider the effect of $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ component by component:

• The ones of $\Phi_k$ are not affected by the zero jump constraint, because their trace is already equal to 0. This induces $N \frac{k(k-1)}{2}$ components in $d\mathbb{P}_k$.

• The piecewise constant components, with the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is a set that was already identified in previous publications [15, 3, 23], and this space is $\nabla^{\perp} \mathbb{P}_1 \oplus \mathbf{K}$, which is of dimension $1 + \frac{N}{2}$.

• Concerning the components $\bar{v}_i^f$ for $1 \le i \le k$, the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is inducing $k \# \mathcal{F}$ free constraints on a space of dimension $3kN$. This induces a space of dimension

$$3kN - k \# \mathcal{F} = 3kN - k \frac{3N}{2} = \frac{3kN}{2}.$$

Adding the dimension of these different sets gives the dimension of $\ker \nabla_{\mathcal{G}'}$:

$$ \begin{aligned} \dim (\ker \nabla_{\mathcal{G}'}) &= N \frac{k(k-1)}{2} + 1 + \frac{N}{2} + \frac{3kN}{2} \\ &= \frac{N}{2} (k^2 - k + 1 + 3k) + 1 \\ &= \frac{N(k+1)^2}{2} + 1, \end{aligned} $$

which is exactly equal to $\mathrm{rank}(\nabla^{\perp}) + \dim \mathbf{K}$. We have then proven that $\mathrm{Range}(\nabla^{\perp}) \oplus \mathbf{K} \subset \ker \nabla_{\mathcal{G}'}$ and that the dimensions are equal, so that $\mathrm{Range}(\nabla^{\perp}) \oplus \mathbf{K} = \ker \nabla_{\mathcal{G}'}$. $\square$

## 4.5 Study of $\nabla_{\mathcal{G}'}$

The kernel of $\nabla_{\mathcal{G}'}$ was already characterised in Proposition 10. We now characterise its range.

**Proposition 11** (Range of $\nabla_{\mathcal{G}'}$). We have

$$d\mathbb{P}_{k-1}(C) \times d\mathbb{P}_k(\mathcal{F}) = \mathrm{Range}(\nabla_{\mathcal{G}'} \oplus \mathbf{K}),$$

where the sum is orthogonal for the scalar product defined in (7).

*Proof.* Following Proposition 10 and Proposition 7 we have

$$\dim (\ker \nabla_{\mathcal{G}'}) = \mathrm{rank} (\nabla^{\perp}) + 2 = \frac{N(k+1)^2}{2} - 1 + 2 = \frac{N(k+1)^2}{2} + 1.$$

Using the rank nullity theorem gives

$$\dim (d\mathbb{P}_k) = \mathrm{rank} (\nabla_{\mathcal{G}'}) + \dim (\ker \nabla_{\mathcal{G}'}).$$

Using Proposition 6 leads to

$$ \begin{aligned} \mathrm{rank}(\nabla_{\mathcal{G}'}) &= \dim d\mathbb{P}_k - \dim (\ker \nabla_{\mathcal{G}'}) \\ &= N(k+1)(k+2) - \left(\frac{N(k+1)^2}{2} + 1\right) \\ &= \frac{N(k+1)(2(k+2) - (k+1))}{2} - 1 \\ &= \frac{N(k+1)(2k+4-k-1)}{2} - 1 \\ \mathrm{rank}(\nabla_{\mathcal{G}'}) &= \frac{N(k+1)(k+3)}{2} - 1. \end{aligned} $$
because both $\mathbf{k}^{\perp}$ and $\psi$ are continuous across the faces. We have then proven that

$$ \mathbf{K} \perp \text{Range}(\nabla^{\perp}). $$

We also remark that $\mathbf{K} \subset \ker(\nabla_{\mathcal{G}'\cdot})$, because $\nabla_{\mathcal{G}'\cdot}$ is a derivation operator. For the moment, we have proven that

$$ \text{Range}(\nabla^{\perp}) \oplus \mathbf{K} \subset \ker \nabla_{\mathcal{G}'\cdot}. $$

Suppose now that an element $\mathbf{u} \in \mathbf{d}\mathbf{P}_k$ is such that its divergence is 0, namely

$$ \left\{ \begin{array}{ll} \forall c \in \mathcal{C} & \nabla \cdot \mathbf{u} = 0 \\ \forall f \in \mathcal{F} & [\mathbf{u} \cdot \mathbf{n}_f] = 0. \end{array} \right. $$

As $\nabla \cdot \mathbf{u} = 0$ on all the cells, $\mathbf{u}$ can be decomposed as in Proposition 10; this decomposition involves three types of components:

* The ones of $\Phi_k$, which have a trace equal to 0. This set is of dimension $\frac{k(k-1)}{2}$ on each cell.

* The constant component; this set is of dimension 2 on each cell.

* The components $\bar{\mathbf{v}}_i^f$ for $1 \le i \le k$ and for all faces. This set is of dimension $3k$ on each cell.

Let us see now what is the effect of the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ on these different components. We first remark that the traces of the different components are orthogonal two at a time, which means that we can consider the effect of $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ component by component:

* The ones of $\Phi_k$ are not affected by the zero jump constraint, because their trace is already equal to 0. This induces $N \frac{k(k-1)}{2}$ components in $\mathbf{d}\mathbf{P}_k$.

* The piecewise constant components, with the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is a set that was already identified in previous publications [18, 3, 26], and this space is $\nabla^{\perp}\mathbf{P}_1 \oplus \mathbf{K}$, which is of dimension $1 + \frac{N}{2}$.

* Concerning the components $\bar{\mathbf{v}}_i^f$ for $1 \le i \le k$, the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is inducing $k\#\mathcal{F}$ free constraints on a space of dimension $3kN$. This induces a space of dimension

$$ 3kN - k\#\mathcal{F} = 3kN - k \frac{3N}{2} = \frac{3kN}{2}. $$

Adding the dimension of these different sets gives the dimension of $\ker \nabla_{\mathcal{G}'\cdot}$:

$$ \begin{aligned} \dim (\ker \nabla_{\mathcal{G}'\cdot}) &= N \frac{k(k-1)}{2} + 1 + \frac{N}{2} + \frac{3kN}{2} \\ &= \frac{N}{2} (k^2 - k + 1 + 3k) + 1 \\ &= \frac{N(k+1)^2}{2} + 1, \end{aligned} $$

which is exactly equal to $\text{rank}(\nabla^{\perp}) + \dim \mathbf{K}$. We have then proven that $\text{Range}(\nabla^{\perp}) \oplus \mathbf{K} \subset \ker \nabla_{\mathcal{G}'\cdot}$ and that the dimensions are equal, so that $\text{Range}(\nabla^{\perp}) \oplus \mathbf{K} = \ker \nabla_{\mathcal{G}'\cdot}$. $\square$

## 4.5 Study of $\nabla_{\mathcal{G}'\cdot}$

The kernel of $\nabla_{\mathcal{G}'\cdot}$ was already characterised in Proposition 11. We now characterise its range.

**Proposition 12** (Range of $\nabla_{\mathcal{G}'\cdot}$). We have

$$ \mathbf{d}\mathbf{P}_{k-1}(\mathcal{C}) \times \mathbf{d}\mathbf{P}_k(\mathcal{F}) = \text{Range}(\nabla_{\mathcal{G}'\cdot}) \oplus \mathbf{K}, $$

where the sum is orthogonal for the scalar product defined in (7).
We prove now that $\mathbb{K}$ is orthogonal to Range $(\nabla_{\mathcal{D}'})$. We denote by $k$ an element of $d\mathbb{P}_{k-1}(\mathcal{C}) \times d\mathbb{P}_k(\mathcal{F})$, which has the same value on all the cells and faces. We also denote by $k$ this value. We denote by $\mathbf{u}$ an element of $d\mathbb{P}_k(\mathcal{C})$. Then

$$
\begin{align*}
\langle \nabla_{\mathcal{D}'} \cdot \mathbf{u} | k \rangle_{[d\mathbb{P}_{k-1}(\mathcal{C}) \times d\mathbb{P}_k(\mathcal{F})]} &= \sum_{c \in \mathcal{C}} \int_c k \nabla \cdot \mathbf{u} + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\
&= \sum_{c \in \mathcal{C}} \int_c \nabla \cdot (k\mathbf{u}) + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\
&= \sum_{c \in \mathcal{C}} \int_{\partial c} k\mathbf{u} \cdot \mathbf{n}_{out} + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\
&= \sum_{c \in \mathcal{C}} \sum_{f \in \mathcal{F}(c)} \int_f k\mathbf{u} \cdot \mathbf{n}_{out} + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\
&= \sum_{f \in \mathcal{F}} \int_f (k\mathbf{u}_L - k\mathbf{u}_R) \cdot \mathbf{n}_f + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\
&= - \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\
\langle \nabla_{\mathcal{D}'} \cdot \mathbf{u} | k \rangle_{[d\mathbb{P}_{k-1}(\mathcal{C}) \times d\mathbb{P}_k(\mathcal{F})]} &= 0.
\end{align*}
$$

We thus have proven that Range $(\nabla_{\mathcal{D}'}) \perp \mathbb{K}$. As rank $(\nabla_{\mathcal{D}'}) = \dim (d\mathbb{P}_{k-1}(\mathcal{C}) \times d\mathbb{P}_k(\mathcal{F})) - 1$, this actually means that

$$d\mathbb{P}_{k-1}(\mathcal{C}) \times d\mathbb{P}_k(\mathcal{F}) = \text{Range}(\nabla_{\mathcal{D}'}) \oplus \mathbb{K},$$

which ends the proof.

<figure><img src="image_9.png" /></figure>

## 4.6 Summary on the de-Rham complex

Gathering all the results of this section, the following proposition was proven

**Proposition 12.** *The discrete diagram*

$$\mathbb{P}_{k+1} \xrightarrow{\nabla^{\perp}} \mathbf{dP}_k(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{D}'}} \mathbb{dP}_k(\mathcal{F}) \times \mathbb{dP}_{k-1}(\mathcal{C}),$$

where $\nabla_{\mathcal{D}'}$ is the $\nabla$ in the sense of distributions, ensures the harmonic gap property given in Definition 1. Moreover

$$\left\{ \begin{array}{l} \mathbb{P}_{k+1}/\mathbb{K} = \ker (\nabla^{\perp}) \\ (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_{k-1}(\mathcal{C})) / k = \text{Range}(\nabla_{\mathcal{D}'}) \end{array} \right..$$

By changing the representation of the linear forms, the following proposition is also obtained:

**Proposition 13.** *The discrete diagram*

$$\mathbb{P}_{k+1} \xrightarrow{\nabla} \mathbf{dP}_k(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} \mathbb{dP}_k(\mathcal{F}) \times \mathbb{dP}_{k-1}(\mathcal{C}),$$

where $\nabla_{\mathcal{D}'}^{\perp}$ is $\nabla^{\perp}$ in the sense of distributions, ensures the harmonic gap property given in Definition 1. Moreover

$$\left\{ \begin{array}{l} \mathbb{P}_{k+1}/\mathbb{K} = \ker (\nabla) \\ (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_{k-1}(\mathcal{C})) / k = \text{Range}(\nabla_{\mathcal{D}'}^{\perp}) \end{array} \right..$$

The location of the degrees of freedom for the two discrete de-Rham complexes found for triangles are summarised in Figure 5. Note that compared with the finite element spaces $d\mathbf{R}T_{k+1}^{\Delta}$ and $d\mathbf{N}_{k+1}^{\Delta}$ discussed in section 3, the space $d\mathbb{P}_k$ represents a significant improvement regarding the number of degrees of freedom, as

$$\dim d\mathbf{R}T_{k+1}^{\Delta} - d\mathbb{P}_k = d\mathbf{N}_{k+1}^{\Delta} - d\mathbb{P}_k = k + 1.$$

Also, Raviart-Thomas and Nédélec finite element basis are known to be difficult to generate on simplices, whereas the generation of a basis for $d\mathbb{P}_k$ is straightforward. Therefore, using the basis $d\mathbb{P}_k$ instead of $d\mathbf{R}T_{k+1}^{\Delta}$ or $d\mathbf{N}_{k+1}^{\Delta}$ for discontinuous approximations seems to be very beneficial.
Proof. Following Proposition 11 and Proposition 8 we have

$$ \dim (\ker \nabla_{\mathcal{D}'\cdot}) = \mathrm{rank} (\nabla^{\perp}) + 2 = \frac{N(k+1)^2}{2} - 1 + 2 = \frac{N(k+1)^2}{2} + 1. $$

Using the rank nullity theorem gives

$$ \dim (d\mathbb{P}_k) = \mathrm{rank} (\nabla_{\mathcal{D}'\cdot}) + \dim (\ker \nabla_{\mathcal{D}'\cdot}). $$

Using Proposition 7 leads to

$$ \begin{aligned} \mathrm{rank} (\nabla_{\mathcal{D}'\cdot}) &= \dim d\mathbb{P}_k - \dim (\ker \nabla_{\mathcal{D}'\cdot}) \\ &= N(k+1)(k+2) - \left( \frac{N(k+1)^2}{2} + 1 \right) \\ &= \frac{N(k+1)(2(k+2) - (k+1))}{2} - 1 \\ &= \frac{N(k+1)(2k+4-k-1)}{2} - 1 \\ \mathrm{rank} (\nabla_{\mathcal{D}'\cdot}) &= \frac{N(k+1)(k+3)}{2} - 1. \end{aligned} $$

We prove now that $K$ is orthogonal to Range $(\nabla_{\mathcal{D}'\cdot})$. We denote by $k$ an element of $d\mathbb{P}_{k-1}(C) \times d\mathbb{P}_k(\mathcal{F})$, which has the same value on all the cells and faces. We also denote by $k$ this value. We denote by $\mathbf{u}$ an element of $d\mathbb{P}_k(C)$. Then

$$ \begin{aligned} \langle \nabla_{\mathcal{D}'\cdot} \cdot \mathbf{u} | k \rangle_{[d\mathbb{P}_{k-1}(C) \times d\mathbb{P}_k(\mathcal{F})]} &= \sum_{c \in \mathcal{C}} \int_c k \nabla \cdot \mathbf{u} + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\ &= \sum_{c \in \mathcal{C}} \int_c \nabla \cdot (k\mathbf{u}) + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\ &= \sum_{c \in \mathcal{C}} \int_{\partial c} k\mathbf{u} \cdot \mathbf{n}_{out} + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\ &= \sum_{c \in \mathcal{C}} \sum_{f \in \mathcal{F}(c)} \int_f k\mathbf{u} \cdot \mathbf{n}_{out} + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\ &= \sum_{f \in \mathcal{F}} \int_f (k\mathbf{u}_L - k\mathbf{u}_R) \cdot \mathbf{n}_f + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \\ &= -\sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket + \sum_{f \in \mathcal{F}} \int_f k \llbracket \mathbf{u} \cdot \mathbf{n}_f \rrbracket \end{aligned} \quad \langle \nabla_{\mathcal{D}'\cdot} \cdot \mathbf{u} | k \rangle_{[d\mathbb{P}_{k-1}(C) \times d\mathbb{P}_k(\mathcal{F})]} = 0. $$

We thus have proven that Range $(\nabla_{\mathcal{D}'\cdot}) \perp K$. As $\mathrm{rank} (\nabla_{\mathcal{D}'\cdot}) = \dim (d\mathbb{P}_{k-1}(C) \times d\mathbb{P}_k(\mathcal{F})) - 1$, this actually means that

$$ d\mathbb{P}_{k-1}(C) \times d\mathbb{P}_k(\mathcal{F}) = \mathrm{Range} (\nabla_{\mathcal{D}'\cdot}) \oplus K, $$

which ends the proof. $\square$

## 4.6 Summary on the de-Rham complex

Gathering all the results of this section, the following proposition was proven

**Proposition 13.** *The discrete diagram*

$$ \mathbb{P}_{k+1} \xrightarrow{\nabla^{\perp}} \mathbf{dP}_k(C) \xrightarrow{\nabla_{\mathcal{D}'\cdot}} d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_{k-1}(C), $$

where $\nabla_{\mathcal{D}'\cdot}$ is the $\nabla$ in the sense of distributions, ensures the Proposition 1. Moreover

$$ \left\{ \begin{array}{l} \mathbb{P}_{k+1}/K = \ker (\nabla^{\perp}) \\ (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_{k-1}(C))/k = \mathrm{Range} (\nabla_{\mathcal{D}'\cdot}). \end{array} \right. $$

By changing the representation of the linear forms, the following proposition is also obtained:
<figure><img src="image_10.png" /><figcaption>Figure 5: Representation of the finite element spaces involved in the grad/curl and curl/div de-Rham complex for triangular meshes for $k = 0, 1, 2$ and $3$. Same code for colors are used as in Figure 1.</figcaption></figure>
**Proposition 14.** *The discrete diagram*

$$ \mathbb{P}_{k+1} \xrightarrow{\nabla} \mathbf{d}\mathbb{P}_k(C) \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} \mathbf{d}\mathbb{P}_k(\mathcal{F}) \times \mathbf{d}\mathbb{P}_{k-1}(C), $$

where $\nabla_{\mathcal{D}'}^{\perp}$ is $\nabla^{\perp}$ in the sense of distributions, ensures the Proposition 1. Moreover

$$ \left\{ \begin{array}{l} \mathbb{P}_{k+1}/\mathbb{K} = \ker(\nabla) \\ (\mathbf{d}\mathbb{P}_k(\mathcal{F}) \times \mathbf{d}\mathbb{P}_{k-1}(C))/k = \text{Range}(\nabla_{\mathcal{D}'}^{\perp}). \end{array} \right. $$

The location of the degrees of freedom for the two discrete de-Rham complexes found for triangles are summarised in Figure 5. Note that compared with the finite element spaces $\mathbf{d}\mathbb{R}\mathbf{T}_{k+1}^{\Delta}$ and $\mathbf{d}\mathbb{N}_{k+1}^{\Delta}$ discussed in section 3, the space $\mathbf{d}\mathbb{P}_k$ represents a significant improvement regarding the number of degrees of freedom, as

$$ \dim \mathbf{d}\mathbb{R}\mathbf{T}_{k+1}^{\Delta} - \mathbf{d}\mathbb{P}_k = \mathbf{d}\mathbb{N}_{k+1}^{\Delta} - \mathbf{d}\mathbb{P}_k = k+1. $$

Also, Raviart-Thomas and Nédélec finite element basis are known to be difficult to generate on simplices, whereas the generation of a basis for $\mathbf{d}\mathbb{P}_k$ is straightforward. Therefore, using the basis $\mathbf{d}\mathbb{P}_k$ instead of $\mathbf{d}\mathbb{R}\mathbf{T}_{k+1}$ or $\mathbf{d}\mathbb{N}_{k+1}$ for discontinuous approximations seems to be very beneficial.

# 5 The case of Cartesian meshes

## 5.1 Why the Cartesian case is more complicated

Inspired by the $\mathbf{d}\mathbb{P}_0$ triangular case, we consider the following discrete divergence

$$ \begin{array}{l} \nabla_{\mathcal{D}'} : \mathbf{d}\mathbb{Q}_0 \mapsto \mathbf{d}\mathbb{P}_0(\mathcal{F}) \\ \mathbf{u} \mapsto a \text{ such that } a_f := [\mathbf{u} \cdot \mathbf{n}_f]. \end{array} \qquad (15) $$

We directly see that $\dim \mathbf{d}\mathbb{Q}_0 = 2\#\mathcal{C} = 2N$ and $\dim \mathbf{d}\mathbb{P}_0(\mathcal{F}) = \#\mathcal{F} = 2N$. If $N_x$ is the number of cells in the $x$ direction and $N_y$ is the number of cells in the $y$ direction, then $N = N_x N_y$. Also, it is easy to see that $\ker \nabla_{\mathcal{D}'}$ is composed of $N_x + N_y$ components. This means that

$$ \mathrm{rank} \nabla_{\mathcal{D}'} = 2N - N_x - N_y, $$

whereas we expect the $\nabla_{\mathcal{D}'}$ to be of rank $2N - 1$. Also, the kernel of the divergence is much smaller than what is expected for correctly approximating continuous divergence free vectors.

A second problem that we see is that when deriving an element of $\mathbb{Q}_k$, it does not give an element of $\mathbb{Q}_{k-1}$. Therefore, the finite element space that should be put before the discrete $\nabla^{\perp}$ operator is difficult to determine. If we put $\mathbb{Q}_0$, then the derivative will be 0 and the range of the discrete $\nabla^{\perp}$ will be reduced to 0. If $\mathbb{Q}_1$ is used, then the $\nabla^{\perp}$ will be in a space larger than $\mathbb{Q}_0$.

## 5.2 Finite element space definition

We need to define some finite element spaces for discretising the different spaces, A, B and C that are involved in the de-Rham diagram (4). Based on what was done in the triangular case (12), we propose to start by the continuous finite element space $\mathbb{Q}_{k+1}$ for the space A. For the space B, the initial plan is to take $\mathbf{d}\mathbb{Q}_k$, but we need to enrich it with the curl of $\mathbb{Q}_{k+1}$. For simplifying, we relax the continuity conditions on $\nabla^{\perp}\mathbb{Q}_{k+1}$, and add to $\mathbf{d}\mathbb{Q}_k$ all the piecewise discontinuous polynomials that have the same degree as the one of $\nabla^{\perp}\mathbb{Q}_{k+1}$. This leads to choose $\widehat{\mathbf{d}\mathbb{Q}}_k^{\text{curl}}(\mathcal{C})$ for discretizing the space B. Last, taking the divergence in the sense of distributions for $\widehat{\mathbf{d}\mathbb{Q}}_k^{\text{curl}}(\mathcal{C})$ naturally maps to $\mathbf{d}\mathbb{Q}_k(\mathcal{F}) \times \widehat{\mathbf{d}\mathbb{Q}}_{k-1}(\mathcal{C})$. This is why the following diagram is considered

$$ \mathbb{Q}_{k+1} \xrightarrow{\nabla^{\perp}} \widehat{\mathbf{d}\mathbb{Q}}_k^{\text{curl}}(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} \mathbf{d}\mathbb{Q}_k(\mathcal{F}) \times \widehat{\mathbf{d}\mathbb{Q}}_{k-1}(\mathcal{C}). \qquad (16) $$
# 5 The case of Cartesian meshes

## 5.1 Why the Cartesian case is more complicated

Inspired by the $d\mathbb{P}_0$ triangular case, we consider the following discrete divergence

$$ \nabla_{\mathcal{D}'} : \begin{array}{rcl} \mathbf{d}\mathbb{Q}_0 & \mapsto & d\mathbb{P}_0(\mathcal{F}) \\ \mathbf{u} & \mapsto & a \text{ such that } a_f := [\mathbf{u} \cdot \mathbf{n}_f]. \end{array} \qquad (15) $$

We directly see that $\dim \mathbf{d}\mathbb{Q}_0 = 2\#\mathcal{C} = 2N$ and $\dim d\mathbb{P}_0(\mathcal{F}) = \#\mathcal{F} = 2N$. If $N_x$ is the number of cells in the $x$ direction and $N_y$ is the number of cells in the $y$ direction, then $N = N_x N_y$. Also, it is easy to see that $\ker \nabla_{\mathcal{D}'}$ is composed of $N_x + N_y$ components. This means that

$$ \mathrm{rank} \nabla_{\mathcal{D}'} = 2N - N_x - N_y, $$

whereas we expect the $\nabla_{\mathcal{D}'}$ to be of rank $2N - 1$. Also, the kernel of the divergence is much smaller than what is expected for correctly approximating continuous divergence free vectors.

A second problem that we see is that when deriving an element of $\mathbb{Q}_k$, it does not give an element of $\mathbb{Q}_{k-1}$. Therefore, the finite element space that should be put before the discrete $\nabla^\perp$ operator is difficult to determine. If we put $\mathbb{Q}_0$, then the derivative will be 0 and the range of the discrete $\nabla^\perp$ will be reduced to 0. If $\mathbb{Q}_1$ is used, then the $\nabla^\perp$ will be in a space larger than $\mathbb{Q}_0$.

## 5.2 Finite element space definition

We need to define some finite element spaces for discretising the different spaces, A, B and C that are involved in the de-Rham diagram (4). Based on what was done in the triangular case (12), we propose to start by the continuous finite element space $\mathbb{Q}_{k+1}$ for the space A. For the space B, the initial plan is to take $\mathbf{d}\mathbb{Q}_k$, but we need to enrich it with the curl of $\mathbb{Q}_{k+1}$. For simplifying, we relax the continuity conditions on $\nabla^\perp \mathbb{Q}_{k+1}$, and add to $\mathbf{d}\mathbb{Q}_k$ all the piecewise discontinuous polynomials that have the same degree as the one of $\nabla^\perp \mathbb{Q}_{k+1}$. This leads to choose $\widehat{\mathbf{d}\mathbb{Q}}_k^{\mathrm{curl}}(\mathcal{C})$ for discretizing the space B. Last, taking the divergence in the sense of distributions for $\widehat{\mathbf{d}\mathbb{Q}}_k^{\mathrm{curl}}(\mathcal{C})$ naturally maps to $d\mathbb{Q}_k(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(\mathcal{C})$. This is why the following diagram is considered

$$ \mathbb{Q}_{k+1} \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} \widehat{\mathbf{d}\mathbb{Q}}_k^{\mathrm{curl}}(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{D}'}^{\perp}} d\mathbb{Q}_k(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(\mathcal{C}). \qquad (16) $$

## 5.3 Dimension of the finite element spaces

We first compute the dimension of each of the finite element spaces involved in (16), induced by the number of faces, points and cells that was computed in Proposition 2.

**Proposition 14** (Dimension of the finite element spaces).

$$ \left\{ \begin{array}{l} \dim \mathbb{Q}_{k+1} = N(k+1)^2 \\ \dim \widehat{\mathbf{d}\mathbb{Q}}_k^{\mathrm{curl}}(\mathcal{C}) = N(2(k+1)^2 + 2k+1) \\ \dim \left( d\mathbb{Q}_k(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(\mathcal{C}) \right) = N(k^2 + 4k + 2) \end{array} \right. $$

*Proof.* The dimension of $\mathbb{Q}_{k+1}$ was already proven in Proposition 3. We are now interested in $\widehat{\mathbf{d}\mathbb{Q}}_k^{\mathrm{curl}}(\mathcal{C})$. Let us recall that

$$ \widehat{\mathbf{d}\mathbb{Q}}_k^{\mathrm{curl}}(\mathcal{C}) = [(d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k-1,k+1}) \times (d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k+1,k-1})] \oplus \mathrm{Vec} \left( \begin{array}{c} x^k y^{k+1} \\ x^{k+1} y^k \end{array} \right). $$

We first focus on the dimension of $\mathbb{Q}_{k,k} + \mathbb{Q}_{k-1,k+1}$. The elements of $\mathbb{Q}_{k-1,k+1}$ are all in $\mathbb{Q}_{k,k}$, except for the case in which the degree in $y$ is $k+1$. A basis of these polynomials that are in $\mathbb{Q}_{k-1,k+1}$ but not in $\mathbb{Q}_{k,k}$ is therefore given by $x^i y^{k+1}$ for $0 \le i \le k-1$, so that

$$ \dim (\mathbb{Q}_{k,k} + \mathbb{Q}_{k-1,k+1}) = (k+1)^2 + k. $$
<figure><img src="image_9.png" /><figcaption>Figure 5: Representation of the finite element spaces involved in the gradcurl and curl/div de-
Rham complex for triangular meshes for k = 0, 1, 2 and 3. Same code for colors are used as in
Figure 1.</figcaption></figure>
The space $\mathbb{Q}_{k,k} + \mathbb{Q}_{k+1,k-1}$ has the same dimension, so that

$$ \dim \widehat{\mathbb{d}\mathbb{Q}}_k^{\text{curl}}(\mathcal{C}) = N (\dim (\mathbb{Q}_{k,k} + \mathbb{Q}_{k-1,k+1}) + \dim (\mathbb{Q}_{k,k} + \mathbb{Q}_{k+1,k-1}) + 1) \\ = N (2 ((k+1)^2 + k) + 1) $$

$$ \dim \widehat{\mathbb{d}\mathbb{Q}}_k^{\text{curl}}(\mathcal{C}) = N (2(k+1)^2 + 2k + 1). $$

It remains to compute the dimension of $\mathbb{d}\mathbb{Q}_k(\mathcal{F}) \times \widehat{\mathbb{d}\mathbb{Q}}_{k-1}(\mathcal{C})$. We recall that

$$ \widehat{\mathbb{d}\mathbb{Q}}_{k-1}(\mathcal{C}) = \mathbb{d}\mathbb{Q}_{k-1}(\mathcal{C}) + \mathbb{d}\mathbb{Q}_{k,k-1}(\mathcal{C}) + \mathbb{d}\mathbb{Q}_{k-1,k}(\mathcal{C}). $$

This leads to

$$ \dim \widehat{\mathbb{d}\mathbb{Q}}_{k-1}(\mathcal{C}) = N(k^2 + 2k) \\ = Nk(k + 2). $$

$$ \dim (\mathbb{d}\mathbb{Q}_k(\mathcal{F}) \times \widehat{\mathbb{d}\mathbb{Q}}_{k-1}(\mathcal{C})) = \dim (\mathbb{d}\mathbb{Q}_k(\mathcal{F})) + \dim \widehat{\mathbb{d}\mathbb{Q}}_{k-1}(\mathcal{C}) \\ = \#\mathcal{F}(k+1) + Nk(k+2) \\ = 2N(k+1) + Nk(k+2) $$

$$ \dim (\mathbb{d}\mathbb{Q}_k(\mathcal{F}) \times \widehat{\mathbb{d}\mathbb{Q}}_{k-1}(\mathcal{C})) = N(k^2 + 4k + 2). $$

□

## 5.4 Study of $\nabla_{\mathcal{G}}^{\perp}$

We are now interested in the study of the $\nabla_{\mathcal{G}}^{\perp}$ operator.

**Proposition 15** ($\nabla_{\mathcal{G}}^{\perp}$ in the quadrangular case). We have

$$ \left\{ \begin{array}{l} \dim (\ker \nabla_{\mathcal{G}}^{\perp}) = 1 \\ \text{rank} (\nabla_{\mathcal{G}}^{\perp}) = N(k+1)^2 - 1. \end{array} \right. $$

The proof is exactly the same as Proposition 7.

## 5.5 Discrete divergence free polynomials on the reference cell

We consider the following application

$$ \mathbf{C}_k^{\partial} : \psi \in \mathbb{d}\mathbb{Q}_{k+1}(\hat{K}) \mapsto \text{Tr}(\nabla^{\perp}\psi) \in R_k(\partial\hat{K}) \quad (17) $$

**Proposition 16.** We denote by $k$ the constant elements of $R_k(\partial\hat{K})$. Then

$$ R_k(\partial\hat{K}) = \text{Range } \mathbf{C}_k^{\partial} \oplus k, $$

where the sum is orthogonal.

*Proof.* The proof that the sum of Range $\mathbf{C}_k^{\partial}$ and $k$ is direct and orthogonal follows exactly the same lines as the proof for the triangular case of Proposition 8.

We are now interested in the study of the kernel of $\mathbf{C}_k^{\partial}$. Suppose that an element $\psi$ is such that $\mathbf{C}_k^{\partial}(\psi) = 0$. We consider the classical Lagrange basis of $\mathbb{d}\mathbb{Q}_{k+1}(\hat{K})$. Then $\psi$ is such that its value on the boundary of $\hat{K}$ is constant, and may take any value on the degrees of freedom matching with the interior nodes. This means that

$$ \dim (\ker \mathbf{C}_k^{\partial}) = 1 + k^2. $$

Using the rank-nullity theorem gives

$$ \begin{aligned} \text{rank} (\mathbf{C}_k^{\partial}) &= \dim \mathbb{d}\mathbb{Q}_{k+1} - (1 + k^2) \\ &= (k + 2)^2 - 1 - k^2 \\ &= k^2 + 4k + 4 - 1 - k^2 \end{aligned} $$

$$ \text{rank} (\mathbf{C}_k^{\partial}) = 4k + 3. $$

We also know that $\dim R_k(\partial\hat{K}) = 4(k+1)$. We have then Range $\mathbf{C}_k^{\partial} \oplus k \subset R_k(\partial\hat{K})$, and $\dim (\text{Range } \mathbf{C}_k^{\partial} \oplus k) = \dim R_k(\partial\hat{K})$, so that Range $\mathbf{C}_k^{\partial} \oplus k = R_k(\partial\hat{K})$, which ends the proof.

□
## 5.3 Dimension of the finite element spaces

We first compute the dimension of each of the finite element spaces involved in (16), induced by the number of faces, points and cells that was computed in Proposition 3.

**Proposition 15** (Dimension of the finite element spaces).

$$ egin{cases} \dim \mathbb{Q}_{k+1} = N(k+1)^2 \\ \dim \widehat{d\mathbb{Q}}_{k}^{\text{curl}}(C) = N(2(k+1)^2 + 2k+1) \\ \dim (d\mathbb{Q}_k(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(C)) = N(k^2 + 4k + 2) \end{cases} $$

*Proof.* The dimension of $\mathbb{Q}_{k+1}$ was already proven in Proposition 4. We are now interested in $\widehat{d\mathbb{Q}}_{k}^{\text{curl}}(C)$. Let us recall that

$$ \widehat{d\mathbb{Q}}_{k}^{\text{curl}}(C) = [(d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k-1,k+1}) \times (d\mathbb{Q}_{k,k} + d\mathbb{Q}_{k+1,k-1})] \oplus \text{Vec} \begin{pmatrix} x^k y^{k+1} \\ x^{k+1} y^k \end{pmatrix}. $$

We first focus on the dimension of $\mathbb{Q}_{k,k} + \mathbb{Q}_{k-1,k+1}$. The elements of $\mathbb{Q}_{k-1,k+1}$ are all in $\mathbb{Q}_{k,k}$, except for the case in which the degree in $y$ is $k+1$. A basis of these polynomials that are in $\mathbb{Q}_{k-1,k+1}$ but not in $\mathbb{Q}_{k,k}$ is therefore given by $x^i y^{k+1}$ for $0 \le i \le k-1$, so that

$$ \dim (\mathbb{Q}_{k,k} + \mathbb{Q}_{k-1,k+1}) = (k+1)^2 + k. $$

The space $\mathbb{Q}_{k,k} + \mathbb{Q}_{k+1,k-1}$ has the same dimension, so that

$$ \begin{aligned} \dim \widehat{d\mathbb{Q}}_{k}^{\text{curl}}(C) &= N(\dim (\mathbb{Q}_{k,k} + \mathbb{Q}_{k-1,k+1}) + \dim (\mathbb{Q}_{k,k} + \mathbb{Q}_{k+1,k-1}) + 1) \\ &= N(2((k+1)^2 + k) + 1) \end{aligned} $$

$$ \dim \widehat{d\mathbb{Q}}_{k}^{\text{curl}}(C) = N(2(k+1)^2 + 2k+1). $$

It remains to compute the dimension of $d\mathbb{Q}_k(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(C)$. We recall that

$$ \widehat{d\mathbb{Q}}_{k-1}(C) = d\mathbb{Q}_{k-1}(C) + d\mathbb{Q}_{k,k-1}(C) + d\mathbb{Q}_{k-1,k}(C). $$

This leads to

$$ \begin{aligned} \dim \widehat{d\mathbb{Q}}_{k-1}(C) &= N(k^2 + 2k) \\ &= Nk(k+2). \end{aligned} $$

We finally find

$$ \begin{aligned} \dim (d\mathbb{Q}_k(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(C)) &= \dim (d\mathbb{Q}_k(\mathcal{F})) + \dim \widehat{d\mathbb{Q}}_{k-1}(C) \\ &= \#\mathcal{F}(k+1) + Nk(k+2) \\ &= 2N(k+1) + Nk(k+2) \end{aligned} $$

$$ \dim (d\mathbb{Q}_k(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(C)) = N(k^2 + 4k + 2). $$

□

## 5.4 Study of $\nabla_{\mathcal{G}}^{\perp}$

We are now interested in the study of the $\nabla_{\mathcal{G}}^{\perp}$ operator.

**Proposition 16** ($\nabla_{\mathcal{G}}^{\perp}$ in the quadrangular case). We have

$$ \left\{ \begin{array}{l} \dim (\ker \nabla_{\mathcal{G}}^{\perp}) = 1 \\ \operatorname{rank} (\nabla_{\mathcal{G}}^{\perp}) = N(k+1)^2 - 1. \end{array} \right. $$

The proof is exactly the same as Proposition 8.
**Proposition 17** (Decomposition of divergence free elements). We denote by $\mathcal{L}^{f,i}$ the Legendre polynomial of degree $i$ on the face $f$ of $\hat{K}$, normalised such that

$$ \int_{\partial \hat{K}} (\mathcal{L}^{f,i})^2 = 1. $$

Suppose that $\mathbf{u} \in \hat{\mathbf{dQ}}_k^{\text{curl}}(\hat{K})$ is divergence free. Then $\mathbf{u}$ can be uniquely decomposed as

$$ \mathbf{u} = \bar{\mathbf{v}}_{\mathbf{u}}^{0} + \bar{\mathbf{v}}_{\mathbf{u}}^{1} + \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \bar{\mathbf{v}}_{\mathbf{u}}^{f,i}, $$

where

* $\mathbf{v}_{\mathbf{u}}^0$ is in the set

$$ \Phi_k = \{ \mathbf{u} \in \hat{\mathbf{dQ}}_k^{\text{curl}} \quad \nabla \cdot \mathbf{u} = 0 \quad \text{Tr}(\mathbf{u}) = 0 \}. $$

* $\bar{\mathbf{v}}_{\mathbf{u}}^1$ is in $\text{Vec}((\begin{pmatrix} 1 \\ 0 \end{pmatrix}, (\begin{pmatrix} 0 \\ 1 \end{pmatrix}, (\begin{pmatrix} 1-2x \\ 2y-1 \end{pmatrix})).$

* $\bar{\mathbf{v}}_{\mathbf{u}}^{f,i}$ is orthogonal to $\Phi_k$ and such that

$$ - \nabla \cdot \bar{\mathbf{v}}_{\mathbf{u}}^{f,i} = 0. $$

$$ - \text{Tr}(\bar{\mathbf{v}}_{\mathbf{u}}^{f,i}) \text{ is orthogonal to all } \mathcal{L}^{g,j} \text{ for } \{g, j\} \neq \{f, i\}. $$

*Proof.* We first prove that $\bar{\mathbf{v}}_{\mathbf{u}}^1 \perp \Phi_k$. We denote by $\mathbf{u}_{\Phi}$ an element of $\Phi_k$. As in the proof of Proposition 9, a $\psi^{\Phi}$ exists such that $\mathbf{u}_{\Phi} = \nabla^{\perp} \Psi^{\Phi}$, and the trace of $\Psi^{\Phi}$ is constant. Then

$$ \begin{aligned} \int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \mathbf{u}_{\Phi} &= \int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \nabla^{\perp} \Psi^{\Phi} \\ &= \int_{\hat{K}} (\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \cdot \nabla \Psi^{\Phi} \\ &= \int_{\hat{K}} \nabla \cdot ((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \Psi^{\Phi}) - \int_{\hat{K}} \nabla \cdot ((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp}) \Psi^{\Phi} \\ &= \int_{\partial \hat{K}} \text{Tr}((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \Psi^{\Phi}) \\ \int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \mathbf{u}_{\Phi} &= \text{Tr}(\Psi^{\Phi}) \int_{\partial \hat{K}} \text{Tr}((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp}), \end{aligned} $$

because $\nabla \cdot ((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp}) = 0$, and because $\text{Tr}(\Psi^{\Phi})$ is constant. Computing the integral on $\partial \hat{K}$ of all the components of $(\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp}$ (namely $(1,0)^T$, $(0,1)^T$ and $(1-2x, 2y-1)^T$) leads to

$$ \int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \mathbf{u}_{\Phi} = 0. $$

We have then proven that $\bar{\mathbf{v}}_{\mathbf{u}}^1 \perp \Phi_k$. The proof of uniqueness follows exactly the same lines as the proof of Proposition 9.

The strategy for the existence begins similarly but then differs. We consider one $\mathcal{L}^{f,i}$ of $R_k(\partial \hat{K})$, for $i \ge 1$. Then, as was done in the proof of Proposition 9, a $\psi^{f,i}$ exists such that $C_k^{\partial}(\psi^{f,i}) = \mathcal{L}^{f,i}$. Denoting by $\mathcal{P}$ the orthogonal projection on $\Phi_k$, we define $\mathbf{e}^{f,i} := \nabla^{\perp}(\psi^{f,i}) - \mathcal{P}(\nabla^{\perp}(\psi^{f,i}))$. Then $\mathbf{e}^{f,i}$ is orthogonal to $\Phi_k$, is divergence free, and is such that $\text{Tr}(\mathbf{e}^{f,i}) = \mathcal{L}^{f,i}$. We consider now one divergence free $\mathbf{u} \in \hat{\mathbf{dQ}}_k^{\text{curl}}(\hat{K})$. We define

$$ \lambda_{f,i} := \int_{\partial \hat{K}} \text{Tr}(\mathbf{u}) \mathcal{L}^{f,i}, $$

and set $\bar{\mathbf{v}}_{\mathbf{u}}^{f,i} = \lambda_{f,i} \mathbf{e}^{f,i}$. We also set $\bar{\mathbf{v}}_{\mathbf{u}}^0 = \mathcal{P}(\mathbf{u})$, and

$$ \bar{\mathbf{v}}_{\mathbf{u}}^{1} := \mathbf{u} - \bar{\mathbf{v}}_{\mathbf{u}}^{0} - \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \bar{\mathbf{v}}_{\mathbf{u}}^{f,i}. $$
## 5.5 Discrete divergence free polynomials on the reference cell

We consider the following application

$$ C_k^\partial : \psi \in d\mathbb{Q}_{k+1}(\hat{K}) \mapsto \mathrm{Tr}(\nabla^\perp \psi) \in d\mathbb{P}_k(\partial \hat{K}) \quad (17) $$

**Proposition 17.** We denote by $k$ the constant elements of $d\mathbb{P}_k(\partial \hat{K})$. Then

$$ d\mathbb{P}_k(\partial \hat{K}) = \mathrm{Range} \, C_k^\partial \oplus k, $$

where the sum is orthogonal.

*Proof.* The proof that the sum of Range $C_k^\partial$ and $k$ is direct and orthogonal follows exactly the same lines as the proof for the triangular case of Proposition 9.

We are now interested in the study of the kernel of $C_k^\partial$. Suppose that an element $\psi$ is such that $C_k^\partial(\psi) = 0$. We consider the classical Lagrange basis of $d\mathbb{Q}_{k+1}(\hat{K})$. Then $\psi$ is such that its value on the boundary of $\hat{K}$ is constant, and may take any value on the degrees of freedom matching with the interior nodes. This means that

$$ \dim(\ker C_k^\partial) = 1 + k^2. $$

Using the rank-nullity theorem gives

$$ \begin{aligned} \mathrm{rank}(C_k^\partial) &= \dim d\mathbb{Q}_{k+1} - (1 + k^2) \\ &= (k + 2)^2 - 1 - k^2 \\ &= k^2 + 4k + 4 - 1 - k^2 \end{aligned} $$

$$ \mathrm{rank}(C_k^\partial) = 4k + 3. $$

We also know that $\dim d\mathbb{P}_k(\partial \hat{K}) = 4(k+1)$. We have then $\mathrm{Range} \, C_k^\partial \oplus k \subset d\mathbb{P}_k(\partial \hat{K})$, and $\dim(\mathrm{Range} \, C_k^\partial \oplus k) = \dim d\mathbb{P}_k(\partial \hat{K})$, so that $\mathrm{Range} \, C_k^\partial \oplus k = d\mathbb{P}_k(\partial \hat{K})$, which ends the proof. □

**Proposition 18 (Decomposition of divergence free elements).** We denote by $\mathcal{L}^{f,i}$ the Legendre polynomial of degree $i$ on the face $f$ of $\hat{K}$, normalised such that

$$ \int_{\partial \hat{K}} (\mathcal{L}^{f,i})^2 = 1. $$

Suppose that $\mathbf{u} \in \hat{\mathbb{d}}\mathbb{Q}_k^{\mathrm{curl}}(\hat{K})$ is divergence free. Then $\mathbf{u}$ can be uniquely decomposed as

$$ \mathbf{u} = \bar{\mathbf{v}}_{\mathbf{u}}^{0} + \bar{\mathbf{v}}_{\mathbf{u}}^{1} + \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \bar{\mathbf{v}}_{\mathbf{u}}^{f,i}, $$

where

* $\bar{\mathbf{v}}_{\mathbf{u}}^{0}$ is in the set

$$ \Phi_k = \{ \mathbf{u} \in \hat{\mathbb{d}}\mathbb{Q}_k^{\mathrm{curl}} \quad \nabla \cdot \mathbf{u} = 0 \quad \mathrm{Tr}(\mathbf{u}) = 0 \}. $$

* $\bar{\mathbf{v}}_{\mathbf{u}}^{1}$ is in $\mathrm{Vec}((\begin{pmatrix} 1 \\ 0 \end{pmatrix}, (\begin{pmatrix} 0 \\ 1 \end{pmatrix}, (\begin{pmatrix} 1 - 2x \\ 2y - 1 \end{pmatrix})).

* $\bar{\mathbf{v}}_{\mathbf{u}}^{f,i}$ is orthogonal to $\Phi_k$ and such that

$$ - \nabla \cdot \bar{\mathbf{v}}_{\mathbf{u}}^{f,i} = 0. $$

$$ - \mathrm{Tr}(\bar{\mathbf{v}}_{\mathbf{u}}^{f,i}) \text{ is orthogonal to all } \mathcal{L}^{g,j} \text{ for } \{g, j\} \neq \{f, i\}. $$
<figure><img src="image_11.png" /><figcaption>Figure 6: Value of the trace for $\bar{\mathbf{v}}_u^1$ and the basis used for representing it.</figcaption></figure>

Then $\bar{\mathbf{v}}_u^0$ and the $\bar{\mathbf{v}}_u^{f,i}$ ensure all the properties required. It remains to prove that $\bar{\mathbf{v}}_u^1$ is in the expected space. We know that $\bar{\mathbf{v}}_u^1$ is divergence free, orthogonal to $\Phi_k$, and that its trace is constant on each face, because all the components in $\mathcal{L}^{f,i}$ for $i \ge 1$ were removed. We denote by $b_0, b_1, b_2$ and $b_3$ the values of the trace on the boundary of $\bar{\mathbf{v}}_u^1$ (see Figure 6). We then define

$$ \mathbf{k} := \frac{b_1 - b_3}{2} \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \frac{b_2 - b_0}{2} \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \frac{b_0 + b_2}{2} \begin{pmatrix} 1 - 2x \\ 2y - 1 \end{pmatrix}, $$

then $\mathrm{Tr}(\mathbf{k}) = \mathrm{Tr}(\bar{\mathbf{v}}_u^1)$ (note that as $\nabla \cdot \bar{\mathbf{v}}_u^1 = 0$, we have $b_0 + b_1 + b_2 + b_3 = 0$). Also, $\mathbf{k}$ is divergence free and orthogonal to $\Phi_k$. This means that $\bar{\mathbf{v}}_u^1 - \mathbf{k}$ is divergence free and orthogonal to $\Phi_k$. But as $\mathrm{Tr}(\bar{\mathbf{v}}_u^1 - \mathbf{k}) = 0$, $\bar{\mathbf{v}}_u^1 - \mathbf{k}$ is also in $\Phi_k$, and so $\bar{\mathbf{v}}_u^1 - \mathbf{k} = 0$, which proves that $\bar{\mathbf{v}}_u^1$ is in $\mathrm{Vec}(\begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 - 2x \\ 2y - 1 \end{pmatrix})$. This ends the proof. $\square$

It is important to note that the decomposition of Proposition 17 was proven on the reference element, but holds also on all the cells of a Cartesian meshes, the basis function $(\begin{pmatrix} 1 - 2x \\ 2y - 1 \end{pmatrix})$ being replaced by $(\begin{pmatrix} \frac{2}{L_y^x} & (m_x - x) \\ \frac{2}{L_y} & (y - m_y) \end{pmatrix})$, where $(m_x, m_y)$ is the centre of the cell and $L_x$ and $L_y$ are the size of the cell in the directions $x$, and $y$.

## 5.6 ker $(\nabla_{\mathcal{G}'\cdot})$ and Range $(\nabla^{\perp})$

**Proposition 18.** We denote by $\mathbf{K}$ the set of uniform vectors. Then

$$ \ker (\nabla_{\mathcal{G}'\cdot}) = \mathrm{Range} (\nabla^{\perp}) \oplus \mathbf{K}. $$

*Proof.* The proof of Range $(\nabla^{\perp})$ and $\mathbf{K}$ being orthogonal and the two being subvectorial spaces of $\ker (\nabla_{\mathcal{G}'\cdot})$ is exactly the same as in the proof of Proposition 10.

Concerning the dimension of this space, we proceed as in the proof of Proposition 10, but need to rewrite it because the dimensions are different. Suppose now that an element $\mathbf{u} \in \widehat{\mathbf{d}\mathbf{Q}}_k^{\mathrm{curl}}$ is such that its divergence $\nabla_{\mathcal{G}'\cdot}$ is 0, namely

$$ \left\{ \begin{array}{ll} \forall c \in \mathcal{C} & \nabla \cdot \mathbf{u} = 0 \\ \forall f \in \mathcal{F} & [\mathbf{u} \cdot \mathbf{n}_f] = 0. \end{array} \right. $$

As $\nabla \cdot \mathbf{u} = 0$ on all the cells, $\mathbf{u}$ can be decomposed as in Proposition 17; this decomposition involves three types of components:

* The ones of $\Phi_k$, which have a trace equal to 0. This set is of dimension $k^2 + 1$ on each cell.
<figure><img src="image_10.png" /><figcaption>Figure 6: Value of the trace for $\bar{\mathbf{v}}_u^1$ and the basis used for representing it.</figcaption></figure>

*Proof.* We first prove that $\bar{\mathbf{v}}_u^1 \perp \Phi_k$. We denote by $\mathbf{u}_\Phi$ an element of $\Phi_k$. As in the proof of Proposition 10, a $\psi^\Phi$ exists such that $\mathbf{u}_\Phi = \nabla^\perp \Psi^\Phi$, and the trace of $\Psi^\Phi$ is constant. Then

$$ 
\begin{align*}
\int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \mathbf{u}_{\Phi} &= \int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \nabla^{\perp} \Psi^{\Phi} \\
&= \int_{\hat{K}} (\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \cdot \nabla \Psi^{\Phi} \\
&= \int_{\hat{K}} \nabla \cdot ((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \Psi^{\Phi}) - \int_{\hat{K}} \nabla \cdot ((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp}) \Psi^{\Phi} \\
&= \int_{\partial \hat{K}} \mathbf{Tr}((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp} \Psi^{\Phi}) \\
\int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \mathbf{u}_{\Phi} &= \mathbf{Tr}(\Psi^{\Phi}) \int_{\partial \hat{K}} \mathbf{Tr}((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp}),
\end{align*}
$$

because $\nabla \cdot ((\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp}) = 0$, and because $\mathbf{Tr}(\Psi^{\Phi})$ is constant. Computing the integral on $\partial \hat{K}$ of all the components of $(\bar{\mathbf{v}}_{\mathbf{u}}^{1})^{\perp}$ (namely $(1,0)^T$, $(0,1)^T$ and $(1-2x, 2y-1)^T$) leads to

$$ \int_{\hat{K}} \bar{\mathbf{v}}_{\mathbf{u}}^{1} \cdot \mathbf{u}_{\Phi} = 0. $$

We have then proven that $\bar{\mathbf{v}}_u^1 \perp \Phi_k$. The proof of uniqueness follows exactly the same lines as the proof of Proposition 10.

The strategy for the existence begins similarly but then differs. We consider one $\mathcal{L}^{f,i}$ of $d\mathbb{P}_k(\partial \hat{K})$, for $i \ge 1$. Then, as was done in the proof of Proposition 10, a $\psi^{f,i}$ exists such that $C_k^\partial(\psi^{f,i}) = \mathcal{L}^{f,i}$. Denoting by $\mathcal{P}$ the orthogonal projection on $\Phi_k$, we define $\mathbf{e}^{f,i} := \nabla^\perp (\psi^{f,i}) - \mathcal{P}(\nabla^\perp (\psi^{f,i}))$. Then $\mathbf{e}^{f,i}$ is orthogonal to $\Phi_k$, is divergence free, and is such that $\mathbf{Tr}(\mathbf{e}^{f,i}) = \mathcal{L}^{f,i}$. We consider now one divergence free $\mathbf{u} \in \hat{d\mathbb{Q}}_k^{\text{curl}}(\hat{K})$. We define

$$ \lambda_{f,i} := \int_{\partial \hat{K}} \mathbf{Tr}(\mathbf{u}) \mathcal{L}^{f,i}, $$

and set $\bar{\mathbf{v}}_u^{f,i} = \lambda_{f,i} \mathbf{e}^{f,i}$. We also set $\bar{\mathbf{v}}_u^0 = \mathcal{P}(\mathbf{u})$, and

$$ \bar{\mathbf{v}}_{\mathbf{u}}^{1} := \mathbf{u} - \bar{\mathbf{v}}_{\mathbf{u}}^{0} - \sum_{f \in \mathcal{F}(\hat{K})} \sum_{i=1}^{k} \bar{\mathbf{v}}_{\mathbf{u}}^{f,i}. $$

Then $\bar{\mathbf{v}}_u^0$ and the $\bar{\mathbf{v}}_u^{f,i}$ ensure all the properties required. It remains to prove that $\bar{\mathbf{v}}_u^1$ is in the expected space. We know that $\bar{\mathbf{v}}_u^1$ is divergence free, orthogonal to $\Phi_k$, and that its trace is constant on each face, because all the components in $\mathcal{L}^{f,i}$ for $i \ge 1$ were removed. We denote by $b_0, b_1, b_2$ and $b_3$ the values of the trace on the boundary of $\bar{\mathbf{v}}_u^1$ (see Figure 6). We then define

$$ \mathbf{k} := \frac{b_1 - b_3}{2} \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \frac{b_2 - b_0}{2} \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \frac{b_0 + b_2}{2} \begin{pmatrix} 1 - 2x \\ 2y - 1 \end{pmatrix}, $$
• One component in the dimension 3 vectorial space.

$$ \hat{\mathbf{K}} = \mathrm{Vec} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} \frac{2}{L_{i,j}^x} & (m_{i,j}^x - x) \\ \frac{2}{L_{i,j}^y} & (y - m_{i,j}^y) \end{pmatrix} \right). $$

• The components $\bar{\mathbf{v}}_i^f$ for $1 \le i \le k$ and for all faces. This set is of dimension $4k$ on each cell.

Let us see now what is the effect of the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ on these different components. We first remark that the traces of the different components are orthogonal two at a time, which means that we can consider the effect of $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ component by component:

• The ones of $\Phi_k$ are not affected by the zero jump constraint, because their trace is already equal to 0. This induces $N k^2$ components in $\hat{d\mathbb{Q}}_k^{\mathrm{curl}}$.

• The components of $\hat{\mathbf{K}}$, with the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is a set that is identified in Proposition 22, proven in Appendix A: it is of dimension $N + 1$.

• Concerning the components $\bar{\mathbf{v}}_i^f$ for $1 \le i \le k$, the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is inducing $k \# \mathcal{F}$ free constraints on a space of dimension $4kN$. This induces a space of dimension

$$ 4kN - k \# \mathcal{F} = 4kN - k2N = 2kN. $$

Adding the dimension of these different sets gives the dimension of $\ker \nabla_{\mathcal{D}'}$:

$$ \begin{aligned} \dim (\ker \nabla_{\mathcal{D}'}) &= Nk^2 + N + 1 + 2kN \\ &= N(k^2 + 2k + 1) + 1 \\ &= N(k + 1)^2 + 1, \end{aligned} $$

which is exactly equal to $\mathrm{rank}(\nabla^{\perp}) + \dim \mathbf{K}$. We have then proven that $\nabla^{\perp} \oplus \mathbf{K} \subset \ker \nabla_{\mathcal{D}'}$ and that the dimensions are equal, so that $\mathrm{Range}(\nabla^{\perp}) \oplus \mathbf{K} = \ker \nabla_{\mathcal{D}'}$.

□

## 5.7 Study of $\nabla_{\mathcal{D}'}$

The kernel of $\nabla_{\mathcal{D}'}$ was already characterised in Proposition 18. We now characterise its range.

**Proposition 19 (Range of $\nabla_{\mathcal{D}'}$).** We have

$$ \hat{d\mathbb{Q}}_{k-1}(\mathcal{C}) \times d\mathbb{Q}_k(\mathcal{F}) = \mathrm{Range}(\nabla_{\mathcal{D}'} \oplus \mathbf{K}), $$

where the sum is orthogonal for the scalar product defined in (7).

The proof is exactly the same as the proof of Proposition 11.

## 5.8 Summary on the de-Rham complex

Gathering all the results of this section, the following proposition was proven

**Proposition 20.** The discrete diagram

$$ \mathbb{P}_{k+1} \xrightarrow{\nabla^{\perp}} \hat{d\mathbb{Q}}_k^{\mathrm{curl}}(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{D}'}} d\mathbb{Q}_k(\mathcal{F}) \times \hat{d\mathbb{Q}}_{k-1}(\mathcal{C}), $$

where $\nabla_{\mathcal{D}'}$ is $\nabla$ in the sense of distributions, ensures the harmonic gap property given in Definition 1. Moreover

$$ \left\{ \begin{array}{l} \mathbb{Q}_{k+1}/\mathbf{K} = \ker (\nabla_{\mathcal{D}'}^{\perp}) \\ \left( d\mathbb{Q}_k(\mathcal{F}) \times \hat{d\mathbb{Q}}_{k-1}(\mathcal{C}) \right) /k = \mathrm{Range}(\nabla_{\mathcal{D}'}). \end{array} \right. $$

By changing the representation of the linear forms, which is equivalent to rotating of $\pi/2$ the vector spaces, the following proposition is also obtained:
then $\mathrm{Tr}(\mathbf{k}) = \mathrm{Tr}(\bar{\mathbf{v}}_u^1)$ (note that as $\nabla \cdot \bar{\mathbf{v}}_u^1 = 0$, we have $b_0 + b_1 + b_2 + b_3 = 0$). Also, $\mathbf{k}$ is divergence free and orthogonal to $\Phi_k$. This means that $\bar{\mathbf{v}}_u^1 - \mathbf{k}$ is divergence free and orthogonal to $\Phi_k$. But as $\mathrm{Tr}(\bar{\mathbf{v}}_u^1 - \mathbf{k}) = 0$, $\bar{\mathbf{v}}_u^1 - \mathbf{k}$ is also in $\Phi_k$, and so $\bar{\mathbf{v}}_u^1 - \mathbf{k} = 0$, which proves that $\bar{\mathbf{v}}_u^1$ is in $\mathrm{Vec}(\begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 - 2x \\ 2y - 1 \end{pmatrix})$. This ends the proof. □

It is important to note that the decomposition of Proposition 18 was proven on the reference element, but holds also on all the cells of a Cartesian meshes, the basis function $\begin{pmatrix} 1 - 2x \\ 2y - 1 \end{pmatrix}$ being replaced by $\begin{pmatrix} \frac{2}{L_x} & (m_x - x) \\ \frac{2}{L_y} & (y - m_y) \end{pmatrix}$, where $(m_x, m_y)$ is the centre of the cell and $L_x$ and $L_y$ are the size of the cell in the directions $x$, and $y$.

## 5.6 $\ker(\nabla_{\mathcal{G}'\cdot})$ and Range $(\nabla^{\perp})$

**Proposition 19.** We denote by $\mathbb{K}$ the set of uniform vectors. Then

$$ \ker(\nabla_{\mathcal{G}'\cdot}) = \mathrm{Range}(\nabla^{\perp}) \oplus \mathbb{K}. $$

*Proof.* The proof of Range $(\nabla^{\perp})$ and $\mathbb{K}$ being orthogonal and the two being subvectorial spaces of $\ker(\nabla_{\mathcal{G}'\cdot})$ is exactly the same as in the proof of Proposition 11.

Concerning the dimension of this space, we proceed as in the proof of Proposition 11, but need to rewrite it because the dimensions are different. Suppose now that an element $\mathbf{u} \in \hat{\mathbf{dQ}}_k^{\mathrm{curl}}$ is such that its divergence $\nabla_{\mathcal{G}'\cdot}$ is 0, namely

$$ \begin{cases} \forall c \in \mathcal{C} & \nabla \cdot \mathbf{u} = 0 \\ \forall f \in \mathcal{F} & [\mathbf{u} \cdot \mathbf{n}_f] = 0. \end{cases} $$

As $\nabla \cdot \mathbf{u} = 0$ on all the cells, $\mathbf{u}$ can be decomposed as in Proposition 18; this decomposition involves three types of components:

* The ones of $\Phi_k$, which have a trace equal to 0. This set is of dimension $k^2 + 1$ on each cell.

* One component in the dimension 3 vectorial space.

$$ \hat{\mathbb{K}} = \mathrm{Vec} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} \frac{2}{L_{i,j}^x} & (m_{i,j}^x - x) \\ \frac{2}{L_{i,j}^y} & (y - m_{i,j}^y) \end{pmatrix} \right). $$

* The components $\bar{\mathbf{v}}_i^f$ for $1 \le i \le k$ and for all faces. This set is of dimension $4k$ on each cell.

Let us see now what is the effect of the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ on these different components. We first remark that the traces of the different components are orthogonal two at a time, which means that we can consider the effect of $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ component by component:

* The ones of $\Phi_k$ are not affected by the zero jump constraint, because their trace is already equal to 0. This induces $Nk^2$ components in $\hat{\mathbf{dQ}}_k^{\mathrm{curl}}$.

* The components of $\hat{\mathbb{K}}$, with the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is a set that is identified in Proposition 23, proven in Appendix A: it is of dimension $N + 1$.

* Concerning the components $\bar{\mathbf{v}}_i^f$ for $1 \le i \le k$, the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is inducing $k\#\mathcal{F}$ free constraints on a space of dimension $4kN$. This induces a space of dimension

$$ 4kN - k\#\mathcal{F} = 4kN - k2N = 2kN. $$
**Proposition 21.** *The discrete diagram*

$$ \mathbb{Q}_{k+1} \xrightarrow{\nabla} \widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{div}}(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{D}}^{\perp}} \mathbb{d}\mathbb{Q}_{k}(\mathcal{F}) \times \widehat{\mathbb{d}\mathbb{Q}}_{k-1}(\mathcal{C}), $$

where $\nabla_{\mathcal{D}}^{\perp}$ is $\nabla^{\perp}$ in the sense of distributions, ensures the harmonic gap property given in Definition 1. Moreover

$$ \left\{ \begin{array}{l} \mathbb{Q}_{k+1}/\mathbb{K} = \ker(\nabla) \\ \left( \mathbb{d}\mathbb{Q}_{k}(\mathcal{F}) \times \widehat{\mathbb{d}\mathbb{Q}}_{k-1}(\mathcal{C}) \right) / \mathbb{k} = \text{Range}(\nabla_{\mathcal{D}}^{\perp}). \end{array} \right. $$

We now compare the number of degrees of freedom of the basis $\widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{div}}$ and $\widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{curl}}$ with the quadrangular basis discussed in section 3 on each cell:

$$ \dim \widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{div}} - \dim \mathbb{dRT}_{k+1}^{\square} = \dim \widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{curl}} - \dim \mathbb{dN}_{k+1}^{\square} = 2(k+2)(k+1) - (2(k+1)^2 + 2k+1) = 1. $$

Therefore the difference of number of degrees of freedom is negligible. Considering that the vector basis $\widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{div}}$ and $\widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{curl}}$ do not have a Lagrange basis (this is why no representation of these basis was proposed), it seems that using these new basis has few benefits with respect to $\mathbb{dRT}_{k+1}$ and $\mathbb{dN}_{k+1}$.

# 6 Conclusion

In this article, two-dimensional discrete de-Rham structures in which the vector space is a discontinuous approximation space were discussed. We first recalled that by relaxing the normal or tangential continuity properties of the classical conformal space, a set of discontinuous approximation space can be designed as in [25]. These discontinuous spaces, $\mathbb{dN}$ and $\mathbb{dRT}$, are discontinuous versions of the Nédélec N and Raviart-Thomas RT approximation spaces.

Then the de-Rham structure of the natural discontinuous vectorial space $\mathbb{dP}_k$ on triangles, used for example for the discontinuous Galerkin method was investigated. We proved that for straight triangular meshes, a discrete de-Rham complex can be built for which the harmonic gap property is ensured for any order of approximation.

Based on the finite element spaces and discrete $\nabla^{\perp}/(\nabla_{\mathcal{D}}^{\perp})$ or $\nabla/(\nabla_{\mathcal{D}}^{\perp})$ that were used for the triangular case, the same harmonic gap property was proven for discontinuous spaces of vectors. However, the space of vectors is not the classical $\mathbb{d}\mathbb{Q}_k$ approximation space which is usually used in the discontinuous Galerkin method, but rather an enriched version, $\widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{curl}}$ and $\widehat{\mathbb{d}\mathbb{Q}}_{k}^{\text{div}}$, depending on the diagram considered. Note that no diagram that would be based on the so-called serendipity elements was addressed, but this could be a way to derive velocity approximation spaces that can be put in a de-Rham diagram with fewer degrees of freedom (still at the price of getting a nonoptimal order of accuracy on general quads).

It is important to note that only the harmonic gap property was addressed in this article. No approximation property was addressed, which is still far from the framework that was developed in [1] for conforming finite element approximation. Still, the harmonic gap property is an algebraic property that we believe to be useful in the context of the derivation of curl preserving numerical schemes for hyperbolic systems. Some curl-preserving schemes that were developed in the finite volume scheme context [16, 15, 23] rely on the existence of the following discrete decomposition [3] ($\mathbb{CR}$ denotes the Crouzeix-Raviart finite element space [14]), which reads on periodic triangular meshes as

$$ \mathbb{dP}_0/\mathbb{R}^2 = \nabla^{\perp}\mathbb{P}_1 \oplus \nabla\mathbb{CR}, \qquad (18) $$

and on the preservation of the solenoidal component, this property being strongly linked with the correct low Mach number behaviour on triangular and tetrahedral meshes [19, 24]. The harmonic gap property discussed in this article directly induces the following Hodge-Helmholtz decomposition [1]:

$$ \mathbf{B}/\mathbb{R}^2 = \text{Range}(\nabla^{\perp}) \oplus \text{Range}(\nabla_{\mathcal{D}}^{\perp}, \mathcal{S}), \qquad (19) $$
Adding the dimension of these different sets gives the dimension of $\ker \nabla_{\mathcal{G}'}$:

$$ 
\begin{aligned}
\dim (\ker \nabla_{\mathcal{G}'}) &= Nk^2 + N + 1 + 2kN \\
&= N(k^2 + 2k + 1) + 1 \\
&= N(k + 1)^2 + 1,
\end{aligned}
$$

which is exactly equal to $\text{rank}(\nabla^{\perp}) + \dim \mathbb{K}$. We have then proven that $\nabla^{\perp} \oplus \mathbb{K} \subset \ker \nabla_{\mathcal{G}'}$ and that the dimensions are equal, so that $\text{Range}(\nabla^{\perp}) \oplus \mathbb{K} = \ker \nabla_{\mathcal{G}'}$.

□

## 5.7 Study of $\nabla_{\mathcal{G}'}$.

The kernel of $\nabla_{\mathcal{G}'}$ was already characterised in Proposition 19. We now characterise its range.

**Proposition 20 (Range of $\nabla_{\mathcal{G}'}$).** We have

$$ \widehat{d\mathbb{Q}}_{k-1}(\mathcal{C}) \times \widehat{d\mathbb{Q}}_{k}(\mathcal{F}) = \text{Range}(\nabla_{\mathcal{G}'} \oplus \mathbb{K}), $$

where the sum is orthogonal for the scalar product defined in (7).

The proof is exactly the same as the proof of Proposition 12.

## 5.8 Summary on the de-Rham complex

Gathering all the results of this section, the following proposition was proven

**Proposition 21.** The discrete diagram

$$ \mathbb{P}_{k+1} \xrightarrow{\nabla^{\perp}} \widehat{d\mathbb{Q}}_{k}^{\text{curl}}(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{G}'}} \widehat{d\mathbb{Q}}_{k}(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(\mathcal{C}), $$

where $\nabla_{\mathcal{G}'}$ is $\nabla$ in the sense of distributions, ensures the Proposition 1. Moreover

$$ \left\{ \begin{array}{l} \mathbb{Q}_{k+1}/\mathbb{K} = \ker (\nabla_{\mathcal{G}'}) \\ \left( \widehat{d\mathbb{Q}}_{k}(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(\mathcal{C}) \right) /k = \text{Range}(\nabla_{\mathcal{G}'}). \end{array} \right. $$

By changing the representation of the linear forms, which is equivalent to rotating of $\pi/2$ the vector spaces, the following proposition is also obtained:

**Proposition 22.** The discrete diagram

$$ \mathbb{Q}_{k+1} \xrightarrow{\nabla} \widehat{d\mathbb{Q}}_{k}^{\text{div}}(\mathcal{C}) \xrightarrow{\nabla_{\mathcal{G}'}^{\perp}} \widehat{d\mathbb{Q}}_{k}(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(\mathcal{C}), $$

where $\nabla_{\mathcal{G}'}^{\perp}$ is $\nabla^{\perp}$ in the sense of distributions, ensures the Proposition 1. Moreover

$$ \left\{ \begin{array}{l} \mathbb{Q}_{k+1}/\mathbb{K} = \ker (\nabla) \\ \left( \widehat{d\mathbb{Q}}_{k}(\mathcal{F}) \times \widehat{d\mathbb{Q}}_{k-1}(\mathcal{C}) \right) /k = \text{Range}(\nabla_{\mathcal{G}'}^{\perp}). \end{array} \right. $$

We now compare the number of degrees of freedom of the basis $\widehat{d\mathbb{Q}}_{k}^{\text{div}}$ and $\widehat{d\mathbb{Q}}_{k}^{\text{curl}}$ with the quadrangular basis discussed in section 3 on each cell:

$$ \dim \widehat{d\mathbb{Q}}_{k}^{\text{div}} - \dim \widehat{d\mathbb{R}}_{k+1}^{\text{div}} = \dim \widehat{d\mathbb{Q}}_{k}^{\text{curl}} - \dim \widehat{d\mathbb{N}}_{k+1}^{\text{div}} = 2(k+2)(k+1) - (2(k+1)^2 + 2k+1) = 1. $$

Therefore the difference of number of degrees of freedom is negligible. Considering that the vector basis $\widehat{d\mathbb{Q}}_{k}^{\text{div}}$ and $\widehat{d\mathbb{Q}}_{k}^{\text{curl}}$ do not have a Lagrange basis (this is why no representation of these basis was proposed), it seems that using these new basis has few benefits with respect to $\widehat{d\mathbb{R}}_{k+1}$ and $\widehat{d\mathbb{N}}_{k+1}$.
where the ⋆ denotes the adjoint operator. For example, the diagram of Proposition 12 induces the following discrete Hodge-Helmholtz decomposition

$$d\mathbb{P}_k(C)/\mathbb{R}^2 = \text{Range}(\nabla^{\perp} \mathbb{P}_{k+1}) \oplus \text{Range}(\nabla_{\mathcal{D}'}^{\perp,*} (d\mathbb{P}_k(\mathcal{F}) \times d\mathbb{P}_{k-1}(C))), \quad (20)$$

which can be seen as the high order extension of (18). We are currently investigating the use of (19) for the derivation of high order discontinuous Galerkin schemes that preserve a curl. Investigations include also discontinuous Galerkin schemes that preserve a divergence.

## References

[1] Douglas Norman Arnold. Finite element exterior calculus. SIAM, 2018.

[2] Douglas Norman Arnold, Daniele Boffi, and Richard Steven Falk. Quadrilateral H (div) finite elements. SIAM Journal on Numerical Analysis, 42(6):2429–2451, 2005.

[3] Douglas Norman Arnold and Richard Steven Falk. A uniformly accurate finite element method for the Reissner-Mindlin plate. SIAM J. Numer. Anal., 26:1276–1290, 1989.

[4] Douglas Norman Arnold, Richard Steven Falk, and Ragnar Winther. Finite element exterior calculus, homological techniques, and applications. Acta numerica, 15:1–155, 2006.

[5] Douglas Norman Arnold, Richard Steven Falk, and Ragnar Winther. Finite element exterior calculus: from Hodge theory to numerical stability. Bulletin of the American mathematical society, 47(2):281–354, 2010.

[6] Douglas Norman Arnold and Anders Bernhard Logg. Periodic table of the finite elements. SIAM News, 47(9):212, 2014.

[7] Jérôme Bonelle. Compatible discrete operator schemes on polyhedral meshes for elliptic and Stokes equations. Ph. D. Thesis, 2014.

[8] Jérôme Bonelle and Alexandre Ern. Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes. ESAIM: Mathematical Modelling and Numerical Analysis, 48(2):553–581, 2014.

[9] Jérôme Bonelle and Alexandre Ern. Analysis of compatible discrete operator schemes for the Stokes equations on polyhedral meshes. IMA Journal of numerical analysis, 35(4):1672–1697, 2015.

[10] Alain Bossavit. Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism. IEE Proceedings A (Physical Science, Measurement and Instrumentation, Management and Education, Reviews), 135(8):493–500, 1988.

[11] Alain Bossavit. Computational electromagnetism: variational formulations, complementarity, edge elements. Academic Press, 1998.

[12] Franco Brezzi, Jim Douglas, and Luisa Donatella Marini. Two families of mixed finite elements for second order elliptic problems. Numerische Mathematik, 47(2):217–235, 1985.

[13] Franco Brezzi and Michel Fortin. Mixed and hybrid finite element methods, volume 15. Springer Science & Business Media, 2012.

[14] Michel Crouzeix and Pierre-Arnaud Raviart. Conforming and nonconforming finite element methods for solving the stationary Stokes equations I. Revue française d'automatique informatique recherche opérationnelle. Mathématique, 7(R3):33–75, 1973.

[15] Stéphane Dellacherie, Jonathan Jung, Pascal Omnes, and Pierre-Arnaud Raviart. Construction of modified Godunov type schemes accurate at any Mach number for the compressible Euler system. Mathematical Models and Methods in Applied Sciences, November 2016.

[16] Stéphane Dellacherie, Pascal Omnes, and Felix Rieper. The influence of cell geometry on the Godunov scheme applied to the linear wave equation. Journal of Computational Physics, 229(14):5315–5338, 2010.
# 6 Conclusion

In this article, two-dimensional discrete de-Rham structures in which the vector space is a discontinuous approximation space were discussed. We first recalled that by relaxing the normal or tangential continuity properties of the classical conformal space, a set of discontinuous approximation space can be designed as in [29]. These discontinuous spaces, dN and dRT, are discontinuous versions of the Nédélec N and Raviart-Thomas RT approximation spaces.

Then the de-Rham structure of the natural discontinuous vectorial space dP_k on triangles, used for example for the discontinuous Galerkin method was investigated. We proved that for straight triangular meshes, a discrete de-Rham complex can be built for which the Proposition 1 is ensured for any order of approximation.

Based on the finite element spaces and discrete $\nabla^{\perp}/(\nabla_{\mathcal{D}'}')$ or $\nabla/(\nabla_{\mathcal{D}'}')$ that were used for the triangular case, Proposition 1 was proven for discontinuous spaces of vectors. However, the space of vectors is not the classical dQ_k approximation space which is usually used in the discontinuous Galerkin method, but rather an enriched version, $\hat{dQ}_k^{\text{curl}}$ and $\hat{dQ}_k^{\text{div}}$, depending on the diagram considered. Note that no diagram that would be based on the so-called serendipity elements was addressed, but this could be a way to derive velocity approximation spaces that can be put in a de-Rham diagram with fewer degrees of freedom (still at the price of getting a nonoptimal order of accuracy on general quads).

It is important to note that only Proposition 1 was addressed in this article. The bounded cochain projection property was not addressed, which is still far from the framework that was developed in [1] for conforming finite element approximation. Still, Proposition 1 is an algebraic property that we believe to be useful in the context of the derivation of curl preserving numerical schemes for hyperbolic systems. Some curl-preserving schemes that were developed in the finite volume scheme context [19, 18, 26] rely on the existence of the following discrete decomposition [3] (CR denotes the Crouzeix-Raviart finite element space [17]), which reads on periodic triangular meshes as

$$dP_0/\mathbb{R}^2 = \nabla^{\perp} P_1 \oplus \nabla CR, \qquad (18)$$

and on the preservation of the solenoidal component, this property being strongly linked with the correct low Mach number behaviour on triangular and tetrahedral meshes [22, 27]. The Proposition 1 discussed in this article directly induces the following Hodge-Helmholtz decomposition [1]:

$$B/\mathbb{R}^2 = \text{Range}(\nabla^{\perp}) \oplus \text{Range}(\nabla_{\mathcal{D}'}^{\perp,*}), \qquad (19)$$

where the ⋆ denotes the adjoint operator. For example, the diagram of Proposition 13 induces the following discrete Hodge-Helmholtz decomposition

$$dP_k(C)/\mathbb{R}^2 = \text{Range}(\nabla^{\perp} P_{k+1}) \oplus \text{Range}(\nabla_{\mathcal{D}'}^{\perp,*} (dP_k(F) \times dP_{k-1}(C))), \qquad (20)$$

which can be seen as the high order extension of (18). In a submitted article [31], we explain how to preserve a curl or a divergence constraint of a hyperbolic system of conservation law with the approximation spaces that were discussed in this article.

# References

[1] Douglas Norman Arnold. Finite element exterior calculus. SIAM, 2018.

[2] Douglas Norman Arnold, Daniele Boffi, and Richard Steven Falk. Quadrilateral H (div) finite elements. SIAM Journal on Numerical Analysis, 42(6):2429–2451, 2005.

[3] Douglas Norman Arnold and Richard Steven Falk. A uniformly accurate finite element method for the Reissner-Mindlin plate. SIAM J. Numer. Anal., 26:1276–1290, 1989.

[4] Douglas Norman Arnold, Richard Steven Falk, and Ragnar Winther. Finite element exterior calculus, homological techniques, and applications. Acta numerica, 15:1–155, 2006.

[5] Douglas Norman Arnold, Richard Steven Falk, and Ragnar Winther. Finite element exterior calculus: from Hodge theory to numerical stability. Bulletin of the American mathematical society, 47(2):281–354, 2010.
[17] Daniele Antonio Di Pietro and Jérôme Droniou. The Hybrid High-Order method for polytopal meshes. *Number 19 in Modeling, Simulation and Application*, 2020.

[18] Alexandre Ern and Jean-Luc Guermond. *Finite elements I: Approximation and interpolation*, volume 72. Springer Nature, 2021.

[19] Hervé Guillard. On the behavior of upwind schemes in the low Mach number limit. IV: P0 approximation on triangular and tetrahedral cells. *Computers & Fluids*, 38(10):1969–1972, 2009.

[20] Ralf Hiptmair. Discrete Hodge operators. *Numerische Mathematik*, 90:265–289, 2001.

[21] Ralf Hiptmair. Finite elements in computational electromagnetism. *Acta Numerica*, 11:237–339, 2002.

[22] Qingguo Hong, Yuwen Li, and Jinchao Xu. An extended Galerkin analysis in finite element exterior calculus. *Mathematics of Computation*, 91(335):1077–1106, 2022.

[23] Jonathan Jung and Vincent Perrier. Steady low Mach number flows: identification of the spurious mode and filtering method. *Journal of Computational Physics*, 468:111462, 2022.

[24] Jonathan Jung and Vincent Perrier. Behavior of the discontinuous Galerkin method for compressible flows at low Mach number on triangles and tetrahedrons. *SIAM Journal on Scientific Computing*, 46(1):A452–A482, 2024.

[25] Martin Werner Licht. Complexes of discrete distributional differential forms and their homology theory. *Foundations of Computational Mathematics*, 17(4):1085–1122, 2017.

[26] Riccardo Milani, Jérôme Bonelle, and Alexandre Ern. Artificial compressibility methods for the incompressible Navier-Stokes equations using lowest-order face-based schemes on polytopal meshes. *Computational Methods in Applied Mathematics*, 22(1):133–154, 2022.

[27] Pierre-Arnaud Raviart and Jean-Marie Thomas. A mixed finite element method for 2-nd order elliptic problems. In *Mathematical aspects of finite element methods*, pages 292–315. Springer, 1977.

[28] Pierre-Arnaud Raviart and Jean-Marie Thomas. Primal hybrid finite element methods for 2nd order elliptic equations. *Mathematics of Computation*, 31(138):391–413, 1977.

## A Proof for low order quads

**Proposition 22.** Suppose that a periodic Cartesian mesh is composed of N cells, and that on each cell (of mid point $(m_{i,j}^x, m_{i,j}^y)$ and of length $L_{i,j}^x$ and $L_{i,j}^y$), a vector **u** is in

$$ \mathrm{Vec} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} \frac{2}{L_{i,j}^x} (m_{i,j}^x - x) \\ \frac{2}{L_{i,j}^y} (y - m_{i,j}^y) \end{pmatrix} \right), $$

then the vectorial space such that

$$ \forall f \quad [\mathbf{u} \cdot \mathbf{n}_f] = 0, \qquad (21) $$

is of dimension $N+1$.

*Proof.* We denote by $N_x$ (resp. $N_y$) the number of cells in the $x$ (resp. $y$) direction. For each cell $i,j$, we denote by $\alpha_{i,j}$, $\beta_{i,j}$ and $\gamma_{i,j}$ the coefficients such that

$$ \mathbf{u}_{|c_{i,j}} = \alpha_{i,j} \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \beta_{i,j} \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \gamma_{i,j} \begin{pmatrix} \frac{2}{L_{i,j}^x} (m_{i,j}^x - x) \\ \frac{2}{L_{i,j}^y} (y - m_{i,j}^y) \end{pmatrix}. $$

Then the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is equivalent to the following equations

$$ \forall j \quad \begin{cases} \forall 0 \le i \le N_x - 2 \quad \alpha_{i+1,j} - \alpha_{i,j} = -(\gamma_{i+1,j} + \gamma_{i,j}) \\ \alpha_{0,j} - \alpha_{N_x-1,j} = -(\gamma_{0,j} + \gamma_{N_x-1,j}) \end{cases} \qquad (22) $$
[6] Douglas Norman Arnold and Anders Bernhard Logg. Periodic table of the finite elements. *SIAM News*, 47(9):212, 2014.

[7] Dinshaw S. Balsara. Divergence-free adaptive mesh refinement for magnetohydrodynamics. *Journal of Computational Physics*, 174(2):614–648, 2001.

[8] Dinshaw S. Balsara. Second-order-accurate schemes for magnetohydrodynamics with divergence-free reconstruction. *The Astrophysical Journal Supplement Series*, 151(1):149, 2004.

[9] Dinshaw S. Balsara and Daniel Shields Spicer. A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations. *Journal of Computational Physics*, 149(2):270–292, 1999.

[10] Jérôme Bonelle. Compatible discrete operator schemes on polyhedral meshes for elliptic and Stokes equations. *Ph. D. Thesis*, 2014.

[11] Jérôme Bonelle and Alexandre Ern. Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes. *ESAIM: Mathematical Modelling and Numerical Analysis*, 48(2):553–581, 2014.

[12] Jérôme Bonelle and Alexandre Ern. Analysis of compatible discrete operator schemes for the Stokes equations on polyhedral meshes. *IMA Journal of numerical analysis*, 35(4):1672–1697, 2015.

[13] Alain Bossavit. Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism. *IEE Proceedings A (Physical Science, Measurement and Instrumentation, Management and Education, Reviews)*, 135(8):493–500, 1988.

[14] Alain Bossavit. *Computational electromagnetism: variational formulations, complementarity, edge elements*. Academic Press, 1998.

[15] Franco Brezzi, Jim Douglas, and Luisa Donatella Marini. Two families of mixed finite elements for second order elliptic problems. *Numerische Mathematik*, 47(2):217–235, 1985.

[16] Franco Brezzi and Michel Fortin. *Mixed and hybrid finite element methods*, volume 15. Springer Science & Business Media, 2012.

[17] Michel Crouzeix and Pierre-Arnaud Raviart. Conforming and nonconforming finite element methods for solving the stationary Stokes equations I. *Revue française d'automatique informatique recherche opérationnelle. Mathématique*, 7(R3):33–75, 1973.

[18] Stéphane Dellacherie, Jonathan Jung, Pascal Omnes, and Pierre-Arnaud Raviart. Construction of modified Godunov type schemes accurate at any Mach number for the compressible Euler system. *Mathematical Models and Methods in Applied Sciences*, November 2016.

[19] Stéphane Dellacherie, Pascal Omnes, and Felix Rieper. The influence of cell geometry on the Godunov scheme applied to the linear wave equation. *Journal of Computational Physics*, 229(14):5315–5338, 2010.

[20] Daniele Antonio Di Pietro and Jérôme Droniou. *The Hybrid High-Order method for polytopal meshes*. Number 19 in Modeling, Simulation and Applications series. Springer International Publishing, March 2020.

[21] Alexandre Ern and Jean-Luc Guermond. *Finite elements I: Approximation and interpolation*, volume 72. Springer Nature, 2021.

[22] Hervé Guillard. On the behavior of upwind schemes in the low Mach number limit. IV: P0 approximation on triangular and tetrahedral cells. *Computers & Fluids*, 38(10):1969–1972, 2009.

[23] Ralf Hiptmair. Discrete Hodge operators. *Numerische Mathematik*, 90:265–289, 2001.

[24] Ralf Hiptmair. Finite elements in computational electromagnetism. *Acta Numerica*, 11:237–339, 2002.

[25] Qingguo Hong, Yuwen Li, and Jinchao Xu. An extended Galerkin analysis in finite element exterior calculus. *Mathematics of Computation*, 91(335):1077–1106, 2022.
$$ \forall i \quad \left\{ \begin{array}{l} \forall 0 \le j \le N_y - 2 \qquad \beta_{i,j+1} - \beta_{i,j} = -(\gamma_{i,j+1} + \gamma_{i,j}) \\ \beta_{i,0} - \beta_{i,N_y-1} = -(\gamma_{i,0} + \gamma_{i,N_y-1}) \end{array} \right. \qquad (23) $$

We first consider equation (22). If this equation is seen for each j as a system in $\alpha_{*,j}$, then this system is singular, with a kernel directed by (1, 1, ..., 1), and the right hand side is compatible with this kernel if and only if

$$ \forall j \quad \sum_{i} \gamma_{i,j} = 0. \qquad (24) $$

If this last constraint holds, then if one $\alpha_{*,j}$ is known (for example $\alpha_{0,j}$), then the $\alpha_{i,j}$ are known for all i. This makes $N_y$ parameters for recovering the $\alpha_{i,j}$ once the $\gamma_{i,j}$ are known.

In the same manner, by considering (23), we can prove that if and only if

$$ \forall i \quad \sum_{j} \gamma_{i,j} = 0, \qquad (25) $$

the system in $\beta$ has a $N_x$ parameter family solution determined by the $\beta_{i,0}$.

It remains to study the system on the $\gamma$ coefficients defined by (24), (25). For this system, we consider the coefficients $\gamma_{i,j}$ for $i > 0$ and for $j > 0$ as parameters (this makes $N - (N_x + N_y - 1) = N - N_x - N_y + 1$ parameters). Then the $\gamma_{0,j}$ for $j \ge 1$ are determined by (24), and the $\gamma_{i,0}$ for $i \ge 1$ are determined by (25):

$$ \begin{array}{ll} \forall j \ge 1 & \gamma_{0,j} = -\sum_{i \ge 1} \gamma_{i,j} \\ \forall i \ge 1 & \gamma_{i,0} = -\sum_{j \ge 1} \gamma_{i,j} \end{array} \qquad (26) $$

It remains to determine $\gamma_{0,0}$, which is a priori given by two equations

$$ \gamma_{0,0} = -\sum_{i \ge 1} \gamma_{i,0} $$

$$ \gamma_{0,0} = -\sum_{j \ge 1} \gamma_{0,j} $$

However, if (26) is considered, the two equations give the same value, namely

$$ \gamma_{0,0} = -\sum_{i \ge 1} \sum_{j \ge 1} \gamma_{i,j}. $$

We therefore have been able to express all the $\gamma_{0,j}$ and all the $\gamma_{i,0}$ provided all the $\gamma_{i,j}$ for $i \ge 1$ and $j \ge 1$ are known.

Finally, a basis of the vectorial space determined by (21) was found. Its parameters are

* the $\gamma_{i,j}$ for $i \ge 1$ and $j \ge 1$. These are $N - N_x - N_y + 1$ parameters.

* the $\alpha_{0,j}$ for $j \ge 0$. These are $N_y$ parameters.

* the $\beta_{i,0}$ for $i \ge 0$. These are $N_x$ parameters.

This makes a total of $N + 1$ parameters.

□
[26] Jonathan Jung and Vincent Perrier. Steady low Mach number flows: identification of the spurious mode and filtering method. *Journal of Computational Physics*, 468:111462, 2022.

[27] Jonathan Jung and Vincent Perrier. Behavior of the discontinuous Galerkin method for compressible flows at low Mach number on triangles and tetrahedrons. *SIAM Journal on Scientific Computing*, 46(1):A452–A482, 2024.

[28] Jonathan Jung and Vincent Perrier. A curl preserving finite volume scheme by space velocity enrichment. application to the low mach number accuracy problem. *Journal of Computational Physics*, 515:113252, 2024.

[29] Martin Werner Licht. Complexes of discrete distributional differential forms and their homology theory. *Foundations of Computational Mathematics*, 17(4):1085–1122, 2017.

[30] Riccardo Milani, Jérôme Bonelle, and Alexandre Ern. Artificial compressibility methods for the incompressible Navier–Stokes equations using lowest-order face-based schemes on polytopal meshes. *Computational Methods in Applied Mathematics*, 22(1):133–154, 2022.

[31] Vincent Perrier. Development of discontinuous galerkin methods for hyperbolic systems that preserve a curl or a divergence constraint, 2024.

[32] Pierre-Arnaud Raviart and Jean-Marie Thomas. A mixed finite element method for 2-nd order elliptic problems. In *Mathematical aspects of finite element methods*, pages 292–315. Springer, 1977.

[33] Pierre-Arnaud Raviart and Jean-Marie Thomas. Primal hybrid finite element methods for 2nd order elliptic equations. *Mathematics of Computation*, 31(138):391–413, 1977.

[34] Maurizio Tavelli and Michael Dumbser. A pressure-based semi-implicit space–time discontinuous Galerkin method on staggered unstructured meshes for the solution of the compressible Navier–Stokes equations at all Mach numbers. *Journal of Computational Physics*, 341:341–376, 2017.

## A Proof for low order quads

**Proposition 23.** Suppose that a periodic Cartesian mesh is composed of $N$ cells, and that on each cell (of mid point $(m_{i,j}^x, m_{i,j}^y)$ and of length $L_{i,j}^x$ and $L_{i,j}^y$), a vector $\mathbf{u}$ is in

$$ \mathrm{Vec} \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \begin{pmatrix} \frac{2}{L_{i,j}^x} (m_{i,j}^x - x) \\ \frac{2}{L_{i,j}^y} (y - m_{i,j}^y) \end{pmatrix} \right), $$

then the vectorial space such that

$$ \forall f \quad [\mathbf{u} \cdot \mathbf{n}_f] = 0, \qquad (21) $$

is of dimension $N + 1$.

*Proof.* We denote by $N_x$ (resp. $N_y$) the number of cells in the $x$ (resp. $y$) direction. For each cell $i, j$, we denote by $\alpha_{i,j}$, $\beta_{i,j}$ and $\gamma_{i,j}$ the coefficients such that

$$ \mathbf{u}_{|c_{i,j}} = \alpha_{i,j} \begin{pmatrix} 1 \\ 0 \end{pmatrix} + \beta_{i,j} \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \gamma_{i,j} \begin{pmatrix} \frac{2}{L_{i,j}^x} (m_{i,j}^x - x) \\ \frac{2}{L_{i,j}^y} (y - m_{i,j}^y) \end{pmatrix}. $$

Then the constraint $[\mathbf{u} \cdot \mathbf{n}_f] = 0$ is equivalent to the following equations

$$ \forall j \quad \left\{ \begin{array}{l} \forall 0 \le i \le N_x - 2 \quad \alpha_{i+1,j} - \alpha_{i,j} = -(\gamma_{i+1,j} + \gamma_{i,j}) \\ \alpha_{0,j} - \alpha_{N_x-1,j} = -(\gamma_{0,j} + \gamma_{N_x-1,j}) \end{array} \right. \qquad (22) $$

$$ \forall i \quad \left\{ \begin{array}{l} \forall 0 \le j \le N_y - 2 \quad \beta_{i,j+1} - \beta_{i,j} = -(\gamma_{i,j+1} + \gamma_{i,j}) \\ \beta_{i,0} - \beta_{i,N_y-1} = -(\gamma_{i,0} + \gamma_{i,N_y-1}) \end{array} \right. \qquad (23) $$
We first consider equation (22). If this equation is seen for each j as a system in $\alpha_{*,j}$, then this system is singular, with a kernel directed by $(1, 1, \dots, 1)$, and the right hand side is compatible with this kernel if and only if

$$ \forall j \quad \sum_i \gamma_{i,j} = 0. \qquad (24) $$

If this last constraint holds, then if one $\alpha_{*,j}$ is known (for example $\alpha_{0,j}$), then the $\alpha_{i,j}$ are known for all $i$. This makes $N_y$ parameters for recovering the $\alpha_{i,j}$ once the $\gamma_{i,j}$ are known.

In the same manner, by considering (23), we can prove that if and only if

$$ \forall i \quad \sum_j \gamma_{i,j} = 0, \qquad (25) $$

the system in $\beta$ has a $N_x$ parameter family solution determined by the $\beta_{i,0}$.

It remains to study the system on the $\gamma$ coefficients defined by (24), (25). For this system, we consider the coefficients $\gamma_{i,j}$ for $i > 0$ and for $j > 0$ as parameters (this makes $N - (N_x + N_y - 1) = N - N_x - N_y + 1$ parameters). Then the $\gamma_{0,j}$ for $j \ge 1$ are determined by (24), and the $\gamma_{i,0}$ for $i \ge 1$ are determined by (25):

$$ \begin{aligned} \forall j \ge 1 \quad & \gamma_{0,j} = -\sum_{i \ge 1} \gamma_{i,j} \\ \forall i \ge 1 \quad & \gamma_{i,0} = -\sum_{j \ge 1} \gamma_{i,j} \end{aligned} \qquad (26) $$

It remains to determine $\gamma_{0,0}$, which is a priori given by two equations

$$ \begin{aligned} \gamma_{0,0} &= -\sum_{i \ge 1} \gamma_{i,0} \\ \gamma_{0,0} &= -\sum_{j \ge 1} \gamma_{0,j} \end{aligned} $$

However, if (26) is considered, the two equations give the same value, namely

$$ \gamma_{0,0} = -\sum_{i \ge 1} \sum_{j \ge 1} \gamma_{i,j}. $$

We therefore have been able to express all the $\gamma_{0,j}$ and all the $\gamma_{i,0}$ provided all the $\gamma_{i,j}$ for $i \ge 1$ and $j \ge 1$ are known.

Finally, a basis of the vectorial space determined by (21) was found. Its parameters are

* the $\gamma_{i,j}$ for $i \ge 1$ and $j \ge 1$. These are $N - N_x - N_y + 1$ parameters.

* the $\alpha_{0,j}$ for $j \ge 0$. These are $N_y$ parameters.

* the $\beta_{i,0}$ for $i \ge 0$. These are $N_x$ parameters.

This makes a total of $N + 1$ parameters.

□