Lastly, focusing on the degrees of freedom rather than on each element, one
can rearrange the order of terms so that the above relation turns into

$$ \sum_{\mathcal{S} \in \Sigma^h} \left( \sum_{K \in \tau^h, \mathcal{S} \in \Sigma_{W(K)}} \Phi_{\mathcal{S}}^K(Y^h) + \sum_{\Gamma \in \partial \tau^h, \mathcal{S} \in \Sigma_{W(\Gamma)}} \Phi_{\mathcal{S}}^{\Gamma}(Y^h) \right) = 0. \quad (2.97) $$

The partial residuals being independent from each other, it reduces to asking

$$ \sum_{K \in \tau^h, \mathcal{S} \in \Sigma_{W(K)}} \Phi_{\mathcal{S}}^K(Y^h) + \sum_{\Gamma \in \partial \tau^h, \mathcal{S} \in \Sigma_{W(\Gamma)}} \Phi_{\mathcal{S}}^{\Gamma}(Y^h) = 0, \quad \forall \mathcal{S} \in \Sigma^h, \quad (2.98) $$

allowing in turn to approximate the dynamics (2.61) as simply as

$$ \frac{dY_{\mathcal{S}}^K}{dt} = - \left( \sum_{K \in \tau^h, \mathcal{S} \in \Sigma_{W(K)}} \Phi_{\mathcal{S}}^K(Y^h) + \sum_{\Gamma \in \partial \tau^h, \mathcal{S} \in \Sigma_{W(\Gamma)}} \Phi_{\mathcal{S}}^{\Gamma}(Y^h) \right), \quad \forall \mathcal{S} \in \Sigma^h. \quad (2.99) $$

The values $Y_{\mathcal{S}}^h, \mathcal{S} \in \Sigma^h$ of the discrete solution $Y^h$ in $(\tau^h, W^h, \Sigma^h)$ can therefore be updated without having to know explicitly the basis functions that are associated to the set of degrees of freedom $\Sigma_{W(K)}$ in each space $W(K)$.

<figure><img src="image_29.png" /><figcaption>**Figure 2.24.** Example of a Residual Distribution scheme construction based on an initial primal formulation (from left to right), assuming a globally continuous discretisation setting and a characterisation made entirely through a set of point-values degrees of freedom. The control volumes are greyed, the degrees of freedom are represented with black circles.</figcaption></figure>

Any Residual Distribution scheme therefore considers a dual perspective, using implicitly control volumes $C_{\mathcal{S}}$ that encompass the support of their associated degree of freedom $\mathcal{S} \in \Sigma^h$, whose volume correspond to the sum of the $\beta_{\mathcal{S}}^K|K|$ area portions of each primal cell $K$ contributing to the local dynamic (see Figure 2.25). Hence, in the following, the term *dual mesh* is to be understood in the variational sense that is associated to $W^h$, not geometrically.

<figure><img src="image_30.png" /><figcaption>**Figure 2.25.** Comparison between the geometrical dual (left) and the variational (right) duals associated to the discretisation setting of the above figure. The green points represent the vertices of the control volumes, the black ones the degrees of freedom.</figcaption></figure>
**Practical definition of a Residual Distribution scheme** Given the nature of the update procedure (2.99), defining a Residual Distribution scheme reduces in practice to approximating the total residuals (2.91) and determining a distribution scheme (2.94) guaranteeing (2.93). If those two choices can be done independently, a single comprehensive discretisation can be performed by designing directly the expression of the partial residuals under the more general condition

$$ \sum_{\mathcal{G} \in \Sigma_{W(K)}} \Phi_{\mathcal{G}}^{K}(Y^h) = \Phi^{K}(Y^h) \quad \text{and} \quad \sum_{\mathcal{G} \in \Sigma_{W(\Gamma)}} \Phi_{\mathcal{G}}^{\Gamma}(Y^h) = \Phi^{\Gamma}(Y^h), \quad (2.100) $$

replacing (2.93). Indeed, the residual distribution (2.97) only involves the partial residuals, and the constraint (2.92) would be ensured by (2.100), making the scheme naturally conservative regardless the construction process.

Additional properties can then be enforced by applying further constraints on the distribution coefficients, or by choosing wisely the approximate dynamics. Limiting and stabilising techniques can also be employed [64, 65], either modifying the distribution weights or adding extra stabilisation terms that preserve the conservation constraints. The Residual Distribution framework is therefore extremely flexible and allows us to recast a wide variety of numerical schemes. We there exemplify this feature by reformulating a subclass of the Discontinuous Galerkin schemes described in the above Section 2.4.1, that will be of use when deriving an entropy conservative Flux Reconstruction scheme in the later Section 2.5.

### 2.4.3.1 Recasting a Discontinuous Galerkin scheme

Let us consider again the Discontinuous Galerkin framework presented in Section 2.4.1 and show that any scheme can be reformulated in the Residual Distribution framework whenever the characterisation of the test and trial spaces are chosen identical. To begin with, we start by recalling that the general formulation of a Discontinuous Galerkin method writes

$$ \begin{aligned} & -\sum_{K \in \tau^h} \left( \int_K \nabla \varphi \cdot \mathcal{F}(Y^K) \, dx + \int_{\partial K} \varphi F_n(Y^K, Y^{K-}) \, dx \right) \\ & + \sum_{\Gamma \in \tau^h} \int_{\Gamma} \varphi (\mathcal{F}(Y^{\Gamma}) \cdot n - F_b(Y^{\Gamma}, g)) \, dx = 0, \quad \forall \varphi \in V^h, \quad (2.68) \end{aligned} $$

where $F_n$ and $F_b$ are numerical fluxes consistent with $\mathcal{F}$, and where the approximate solution $Y^h$ is looked for in the trial space

$$ W^h = \{ w \in L^2(\Omega, \mathbb{R}^q) \text{ s.t. } \forall K \in \tau^h, w|_K \in W(K) := (\mathbb{P}_k(K))^d \}. $$

Choosing then the test space $V^h$ as

$$ V^h = \{ v \in L^2(\Omega, \mathbb{R}^q) \text{ s.t. } \forall K \in \tau^h, v|_K \in V(K) := \mathbb{P}_k(K) \}, $$
one can select a same characterisation for the test function and for each component of the solution, all being polynomials of degree $k \in \mathbb{N}$. Thus, defining a set of degrees of freedom $\sum_{\mathbb{P}_k(K)}$ whose values entirely determine functions in $\mathbb{P}_k(K)$ and denoting by $\{\varphi_\sigma\}_{\sigma \in \sum_{\mathbb{P}_k(K)}}$ its associated dual basis in $\mathbb{P}_k(K)$, one can decompose the expression of each term as

$$Y^K = \sum_{\sigma \in \sum_{\mathbb{P}_k(K)}} Y_\sigma^K \varphi_\sigma \quad \text{and} \quad \varphi = \sum_{\sigma \in \sum_{\mathbb{P}_k(K)}} \alpha_\sigma \varphi_\sigma, \qquad (2.101)$$

where $Y_\sigma^K \in \mathbb{R}^q$ stands for the discrete solution values and $\alpha_\sigma \in \mathbb{R}$ for the decomposition coefficients of the test function $\varphi$ on the associated basis of $\mathbb{P}_k(K)$. Thus, as the expression (2.68) is linear in $\varphi$, finding a solution $Y^h$ reduces to finding the discrete values $\{Y_\sigma^K, \sigma \in \sum_{\mathbb{P}_k(K)}, K \in \tau^h\}$ whose associated reconstruction $Y^h$ fulfils

$$\begin{aligned} & -\sum_{K \in \tau^h} \left( \int_K \nabla \varphi_\sigma \cdot \mathcal{F}(Y^K) \,dx + \int_{\partial K} \varphi_\sigma F_n(Y^K, Y^{K-}) \,dx \right) \qquad (2.102) \\ & + \sum_{\Gamma \in \tau^h} \int_{\Gamma} \varphi_\sigma (\mathcal{F}(Y^\Gamma) \cdot n - F_b(Y^\Gamma, g)) \,dx = 0, \quad \forall \sigma \in \sum_{\mathbb{P}_k(K)}. \end{aligned}$$

Comparing then the relations (2.102) with the expression (2.96) of the Residual Distribution framework and using the fact that the degrees of freedom are identical for each component of the solution in $W^h$ and for any test function of $V^h$, we can observe that considering partial residuals of the form

$$\Phi_{\sigma}^{DG, K}(Y^h) = \int_K \nabla \varphi_{\sigma} \cdot \mathcal{F}(Y^h) \,dx + \int_{\partial K} \varphi_{\sigma} F_n(Y^K, Y^{K-}) \,dx \qquad (2.103)$$

$$\Phi_{\sigma}^{DG, \Gamma}(Y^h) = \int_{\Gamma} \varphi_{\sigma} (\mathcal{F}(Y^{\Gamma}) \cdot n - F_b(Y^{\Gamma}, g)) \,dx \qquad (2.104)$$

define a Discontinuous Galerkin scheme in the Residual Distribution framework whenever the basis $\{\varphi_\sigma\}_{\sigma \in \sum_{\mathbb{P}_k(K)}}$ is nodal. Indeed, so defined, the functions $\{\varphi_\sigma\}_{\sigma \in \sum_{\mathbb{P}_k(K)}}$ would form a Lagrangian basis of $\mathbb{P}_k(K)$, ensuring through the relation $\sum_{\sigma \in \sum_{\mathbb{P}_k(K)}} \varphi_\sigma = 1$ that the total residuals

$$\Phi^{DG, K}(Y^h) = \sum_{\sigma \in \sum_{\mathbb{P}_k(K)}} \int_K \nabla \varphi_\sigma \cdot \mathcal{F}(Y^h) \,dx + \int_{\partial K} \varphi_\sigma F_n(Y^h, Y^{h, n}) \,dx \qquad (2.105)$$

$$\Phi^{DG, \Gamma}(Y^h) = \sum_{\sigma \in \sum_{\mathbb{P}_k(\Gamma)}} \int_{\Gamma} \varphi_{\sigma} (\mathcal{F}(Y^h) \cdot n - F_b(Y^h, g)) \,dx, \qquad (2.106)$$

fulfil the conservation relations (2.100). Furthermore, the degrees of freedom being then of pointwise nature, the evolution of the discrete solution is immediately intelligible through the update relations (2.99), and the dual mesh is corresponds exactly to the illustration provided in Figure 2.25. The practical definition of the scheme then only relies on the definition of the numerical fluxes $F_b$ and $F_n$ and of the quadrature rules.
### 2.4.3.2 Residual Distribution scheme properties

Having described the general framework of Residual Distribution schemes and observed the flexibility of its definition by reformulating a subclass of Discontinuous Galerkin scheme, let us state conditions on the residuals so that the designed Residual Distribution scheme provides reliable and accurate approximate solutions. For interested readers, proofs of the presented results can be found in [66, 67, 68] and [69].

First of all, let us recall that the conservation property of a Residual Distribution scheme can be ensured directly at the level of the partial residuals, as already observed in the above developments, requiring the relation (2.100).

**Proposition 2.4.1 Conservative Residual Distribution schemes**

Any Residual Distribution scheme defined by (2.98) is conservative whenever its associated partial residuals fulfil the relations

$$ \sum_{\sigma \in \Sigma_{W(K)}} \Phi_{\sigma}^{K}(Y^{h}) = \Phi^{K}(Y^{h}) \quad \text{and} \quad \sum_{\sigma \in \Sigma_{W(\Gamma)}} \Phi_{\sigma}^{\Gamma}(Y^{h}) = \Phi^{\Gamma}(Y^{h}) \quad (2.100) $$

for consistent approximations

$$ \left\{ \begin{array}{l} \Phi^{K}(Y^{h}) \approx \int_{K} \nabla \cdot \mathcal{F}(Y^{h}) \, dx \\ \\ \Phi^{\Gamma}(Y^{h}) \approx - \int_{\Gamma} \mathcal{F}(Y^{h}) \cdot n + F_{b}(Y^{h}, g) \, dx \end{array} \right. \quad (2.107) $$

of the total residuals.

Once a conservative Residual Distribution scheme has been defined, one can guarantee the convergence of the approximate solution towards a weak solution of the steady analogue of the problem (2.1) by reformulating of the classical Lax-Wendroff theorem in the Residual Distribution framework, whose application only requires the satisfaction of an *a-posteriori* estimate on the partial residuals. The well-definition of any Residual Distribution scheme is therefore ensured by the following theorem.

**Theorem 2.4.1 Lax-Wendroff theorem for Residual Distribution**

Let us consider a sequence of shape-regular meshes $(\tau^h)_{h \to 0}$ on which is defined a Residual Distribution scheme (2.98). Assuming that for any mesh $\tau^h$ and any cell $K \in \tau^h$

* the partial residuals satisfy the conservation relations (2.100),

* the resulting approximate total residuals are consistent with (2.91),

* there exists a constant $M \in \mathbb{R}_+$ so that the *a-posteriori* error estimate

$$ \| \Phi_{\mathcal{G}}^{K}(Y^{K}) \| \le c \sum_{\mathcal{G}', \mathcal{G}' \in \Sigma_{W(K)}} \| Y_{\mathcal{G}}^{K} - Y_{\mathcal{G}'}^{K} \| \quad (2.108) $$

holds for any approximate solution $Y^h \in W^h$ whose discrete counterpart satisfies $\|Y^h\|_{\infty} \le M$ and where the constant $c$ depends only on $M$ and $\tau^h$,
and more generally that

* there exists a global constant $C$ such that for any mesh $\tau^h$ of the sequence, the associated discrete solution satisfies $||Y^h||_\infty \le C$,

* there exists a function $Y \in L^2(\Omega, \mathbb{R}^q)$ such that up to a subsequence,

$$ ||Y^h - Y||_{L^2(\Omega)} \to 0, $$

then the limit $Y = \lim_{h \to 0} Y^h$ exists and is a weak solution to the steady analogue of the problem (2.1).

While the approximate solution to the steady analogue of the problem (2.1) converges towards a weak solution given the hypothesis of Theorem 2.4.1, further conditions on the residuals should be imposed to ensure its physical relevance. Considering here entropy conservation as a selection criterion, we state in Proposition 2.4.2 sufficient conditions for the solution to be entropy conservative with respect to a given entropy function and fluxes pair.

**Proposition 2.4.2** Convergence towards weak entropy solutions

Let us consider $(S, \mathcal{G})$ an entropy pair for the steady analogue of the system (2.1a), $G_n$ a numerical entropy flux whose normal component on each element boundary is consistent with $\mathcal{G} \cdot n$, and $U := \nabla_Y S$ the entropy variable vector associated with the entropy flux $S$. Then, given a Residual Distribution scheme expressed by (2.98) and for which

* the assumptions of Theorem 2.4.1 hold,

* the partial residuals $\{\Phi_{\mathcal{G}}^K\}_{\mathcal{G} \in \Sigma_{W(K)}}$ and $\{\Phi_{\mathcal{G}}^\Gamma\}_{\mathcal{G} \in \Sigma_{W(\Gamma)}}$ further satisfy

$$ \left\{ \sum_{\mathcal{G} \in \Sigma_{W(K)}} \langle U_{\mathcal{G}}, \Phi_{\mathcal{G}}^K \rangle \ge \int_{\partial K} G_n(U^K, U^{K-}) \, d\gamma, \quad \forall K \in \tau^h, \quad (2.109a) \right. $$

$$ \left. \sum_{\sigma \in \Sigma_{W(\Gamma)}} \langle U_{\mathcal{G}}, \Phi_{\mathcal{G}}^{\Gamma} \rangle \ge \int_{\Gamma} G_n(U^{\Gamma}, U_b) - \mathcal{G}(U^{\Gamma}) \cdot n \, d\gamma, \quad \forall \Gamma \in \partial \tau^h, \quad (2.109b) \right. $$

for any shape-regular mesh $\tau^h$ and a given quadrature rule denoted $\int$,

the limit weak solution $Y = \lim_{h \to 0} Y^h$ is entropy stable, i.e. for any $\varphi \in C^1(\Omega)$ such that $\varphi > 0$, it holds

$$ - \int_{\tau^h} \nabla \varphi \cdot \mathcal{G}(Y) \, dx + \int_{\partial \tau^h} \varphi \mathcal{G}(g) \cdot n \, d\gamma \le 0, $$

where $\int$ denotes the scheme's quadrature rule. Moreover, when the equality holds, the scheme and limit weak solution are further entropy conservative.

Having defined the Residual Distribution framework and observed the flexibility it confers to the scheme design, we can now derive residuals expressions that fulfil all the assumptions of Proposition 2.4.1, Proposition 2.4.2, and Theorem 2.4.1, ensuring the resulting scheme to be at least entropy conservative.
### 2.4.4 Methods strengths and drawbacks

To conclude this chapter, let us put in perspective the example of numerical methods and frameworks described in the above three sections with their practical usability and relevance with respect to the nature of the problem that we wish to tackle, described in Section 2.1.

**Discontinuous Galerkin** The first approach presented in Section 2.4.1 consists of designing Discontinuous Galerkin methods that consider a finite element perspective. Their design is thus extremely compact and benefits from an extensive flexibility as both approximations of the flux functions and the solution are considered as cross-cell discontinuous. However, as a side-effect, the practical use of any Discontinuous Galerkin scheme requires a large number of degrees of freedom, further rapidly growing along with any mesh refinement. Moreover, the method is *non-conforming*, the space $V^h$ not being contained in $C_0^\infty(\Omega)$, and is known to be unstable by itself [53, 70]. It thus requires the use of further stabilisation terms [71], and even limiting techniques [72] across the discontinuities that may develop and propagate in the approximate solution.

**Flux Reconstruction** As flexible as the Discontinuous Galerkin framework though being entirely discrete, the Flux Reconstruction approach reduces the impact that may have the discontinuity of the flux to the dynamic approximation by introducing a correction term. This correction being determined on a boundary stencil, the continuity of the comprehensive flux approximation is then retrieved compactly and only requires a very small computational effort. However, the solution and initial flux points being located within the cell, the update procedure remains sensitive to the shape of the lifting functions that propagates the correction amount from the boundary to the solution points. As a consequence, ensuring the stability of the scheme may be delicate and requires to consider the comprehensive discretisation.

**Residual Distribution framework** Very general and encompassing particularly both the Flux Reconstruction and Discontinuous Galerkin frameworks under some weak assumptions, Residual Distribution methods act on the mesh variational dual and can consider either a direct or a variational-based formulation. Both the discrete solution and dynamics are then described through the values of degrees of freedom rather than from coefficients associated to given basis functions. Hence, one does not require the explicit knowledge of the basis functions that are dual to the selected degrees of freedom characterisation for updating the solution, and the framework is not bound to any specific geometry. Especially, any hybrid mesh containing
possibly non-convex polygonal cells can be considered, preventing the need
of designing conformal meshes. However, as in the Flux Reconstruction and
Discontinuous Galerkin framework, additional stabilisation terms may be re-
quired and are thus subject to constraints supported on control volumes that
are made of several primal cells. In addition, the convergence towards a physi-
cally relevant solution is still not ensured by definition of the framework itself,
requiring close investigations and requires for example the addition of further
entropy stabilisation terms.

We now propose to go in this direction and to derive an entropy conservative
scheme that can be expressed either in the Residual Distribution or in the
Flux Reconstruction framework.
## 2.5 An entropy conservative scheme

Making use of the Residual Distribution framework flexibility and results on entropy stability and conservation, we first develop a simple criterion that is sufficient for the reformulated scheme to be at least entropy stable. We then go further and propose a construction of a Residual Distribution scheme yielding through a reciprocal recast an entropy conservative Flux Reconstruction scheme that can be used on any polygonal mesh, thus extending the results of [56, 73, 74, 75, 76] that could be applied only to simplicial meshes or only granting linear stability.

We thus start by reformulating the Flux Reconstruction scheme derived in Section 2.4.2 for its expression to fit the Residual Distribution framework, and ask the inequality relations (2.109) to hold for the resulting residuals. Those relations already providing a criterion for the scheme to be entropy stable, we then focus on deriving a specific residuals expression that ensures the equality case, and thus the entropy conservativeness of the resulting scheme. Lastly, by remembering that the only freedom left in the design of the initial Flux Reconstruction scheme lies in the expression of the correction functions and having in perspective the obtained entropy conservative residuals, we derive a system whose solutions provide correction functions ensuring that the resulting Flux Reconstruction schemes is entropy conservative directly in its original formulation.

### 2.5.1 Reformulation of a Flux Reconstruction scheme

To begin with, let us recast the Flux Reconstruction scheme derived in the two last paragraphs of Section 2.4.2 into the Residual Distribution framework. We use here the initial formulation (2.88) from which the discrete operator $D_F^h$ can be expressed at the solution points $\kappa_Y$ of any element $K \in \tau^h$ as

$$D_F^h(x_m) := -\nabla \cdot (F^{K,D})(x_m) - \sum_{j \in [1, N_C]} (\alpha_j \nabla \cdot h_j(x_m)), \quad \forall x_m \in \Sigma_Y, \quad (2.110)$$

any Flux Reconstruction scheme having an entirely pointwise definition. Furthermore, let us recall that by choosing identical flux and solution points, i.e. by defining $\kappa_F = \kappa_Y$, and up to assuming

$$h_j(x_l) \cdot n = 1 \quad \forall j \in [1, N_C], \forall x_l \in \kappa_C, \quad (2.78)$$

the scheme (2.88) is conservative at any discrete level that quantifies the cell-wise local conservations through quadrature rules whose associated points coincide with the flux points within the cell and with correction points on the boundary. In other words, defining the quadrature points $\kappa_G = \kappa_F$ for an inner
quadrature and $\varkappa_G = \varkappa_C$ for a boundary quadrature, it comes

$$ \sum_{g \in [1, |\varkappa_G|]} \omega_g \left( \nabla \cdot (F_g^{K, D}) + \sum_{j \in [1, N_C]} (\alpha_j \nabla \cdot h_j(x_g)) \right) = 0, \quad \forall K \in \tau^h, \quad (2.79) $$

expression from which a Residual Distribution formulation can be derived. Indeed, the solution $Y^K$ is reconstructed in each element $K \in \tau^h$ as a polynomial of order $k \in \mathbb{N}$ whose characterisation involves exactly $|\varkappa_G| = \dim \mathbb{P}_k(K)$ point-values degrees of freedom. One can therefore define partial residuals that describe the scheme (2.88) by splitting the expression (2.79) on the nodal basis $\varphi = \{\varphi_\sigma\}_{\sigma \in [1, |\mathbb{P}_k(K)|]}$ of $\mathbb{P}_k(\mathbb{R})$ that is associated to the discrete solution points $\varkappa_G$. We then obtain for any cell $K \in \tau^h$ the relation

$$ \sum_{\sigma \in \Sigma_{\mathbb{P}_k(K)}} \left( \sum_{g \in [1, |\varkappa_G|]} \omega_g \left( -\varphi_\sigma \nabla \cdot (F_g^{K, D}) + \varphi_\sigma \sum_{j \in [1, N_C]} (\alpha_j \nabla \cdot h_j(x_g)) \right) \right) = 0. \quad (2.111) $$

Thus, choosing for the Residual Distribution scheme a quadrature rule $\int$ that is expressly based on the points $\varkappa_G$ for each cell $K$, we can then define the residuals

$$ \Phi_{\sigma}^{K,FR}(Y^h) := \int_K \varphi_{\sigma} \nabla \cdot F^{K,D}(Y^K) + \varphi_{\sigma} \sum_{j \in [1, N_C]} (\alpha_j \nabla \cdot h_j) dx. \quad (2.112) $$

An analogous reasoning can be performed on the boundary, yielding to defining the boundary residuals as

$$ \Phi_{\sigma}^{\Gamma,FR}(Y^h) = \int_{\Gamma} \varphi_{\sigma} (\mathcal{F}^I(Y^h, g) - F^h(Y^h) \cdot n) \, d\gamma, \quad (2.113) $$

where $\gamma$ represents the path skimming the faces. Hence, a Residual Distribution formulation of the Flux Reconstruction scheme (2.88) reads

$$ \sum_{K \subset \tau^h, \sigma \in K} \Phi_{\sigma}^{K,FR}(Y^h) + \sum_{\Gamma \subset \partial \tau^h, \sigma \in \Gamma} \Phi_{\sigma}^{\Gamma,FR}(Y^h) = 0, \quad (2.114) $$

and is entirely determined by the definition of the lifting functions $\{h_j\}_{j \in [1, N_C]}$. However, even though this last formulation indeed expresses a valid reformulation of the Flux Reconstruction scheme, the discrete conservation principle obtained in (2.84) does not necessarily transfer to the discrete setting considered in (2.100).

**Note** The below development are only focused on adjusting the inner residuals for the scheme to fulfil the inequality (2.109a), and are enough to retrieve an entropy stable Flux Reconstruction scheme. Nevertheless, one should keep in mind that the definition of boundary residuals has a strong impact on the scheme relevance and properties. For interested readers, the detailed analysis [68] can therefore only be recommended.

<figure><img src="image_31.png" /></figure>
## 2.5.2 Discrete conservation of the recast scheme

The derivation of a Flux Reconstruction scheme that is also conservative in its Residual Distribution formulation can be performed by designing a correction function whose expression naturally enforces the conservation of (2.112)–(2.114) in the sense of (2.100). To ease its selection, we propose here a selection criterion based on the residuals expressions. To this aim, we start by noticing that the formulation (2.112) shares an extensive similarity with the Discontinuous Galerkin scheme after having been reformulated as (2.102) in Section 2.4.3. Indeed, applying a Gauss formula on the last term of (2.112), it comes

$$ \Phi_{\sigma}^{K,FR}(Y^h) := -\int_K \nabla\varphi_{\sigma} \cdot F^{K,D}(Y^K) \,dx - \int_K \nabla\varphi_{\sigma} \cdot \sum_{j \in [1, N_C]} \alpha_j h_j \,dx \\ + \int_{\partial K} \varphi_{\sigma} \left( \sum_{j \in [1, N_C]} \alpha_j h_j \cdot n + F^{K,D}(Y^K) \cdot n \right) d\gamma, \quad (2.115) $$

Recalling then that $\alpha_j = (\mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n - F^{K,D}(Y_j^K) \cdot n)$, that the quadrature points of the boundary integral are taken as $\varkappa_C$, and that the assumption $h_j(x_l) \cdot n = \delta_{jl}$ holds for any $j \in [1, N_C]$ and any $x_l \in \varkappa_C$, it can be observed that $\sum_{j \in [1, N_C]} \alpha_j h_j(x_g) = 1$ at each boundary quadrature point $x_g \in \varkappa_G$, reducing the expression (2.115) to

$$ \Phi_{\sigma}^{K,FR}(Y^h) := -\int_K \nabla\varphi_{\sigma} \cdot F^{K,D}(Y^K) \,dx - \int_K \nabla\varphi_{\sigma} \cdot \sum_{j \in [1, N_C]} \alpha_j h_j \,dx \\ + \int_{\partial K} \varphi_{\sigma} \mathcal{F}^I(Y^K, Y^{K-}) \,d\gamma, \quad (2.116) $$

expression sharing similarities with the formulation (2.103). Therefore, up to choosing a quadrature rule for the Discontinuous Galerkin scheme that is identical to the one selected for defining the expression (2.112) of the Flux Reconstruction scheme (2.88), one can identify the correction term

$$ r_{\sigma} := -\int_K \nabla\varphi_{\sigma} \cdot \sum_{j \in [1, N_C]} \alpha_j h_j \,dx \quad (2.117) $$

and derive the relation

$$ \Phi_{\sigma}^{K,FR}(Y^h) = \Phi_{\sigma}^{K,DG}(Y^h) + r_{\sigma}. \quad (2.118) $$

Straightforward and not particularly useful in practice, this last formulation however provides a selection criterion over the correction function (2.77) for the Flux Reconstruction scheme to be discretely conservative in the Residual Distribution sense. Indeed, the reformulated Discontinuous Galerkin scheme (2.103) already being conservative, the relation (2.100) combined with the linearity of the integral operator implies that the total contribution of the correction terms to a single degrees of freedom should vanish, forcing the
terms $r_σ$ to only modify the distribution of the total residuals locally in their associated primal element. In other words, we require

$$ \sum_{\sigma \in \Sigma_{\mathbb{P}_k}(K)} r_{\sigma} = - \sum_{\sigma \in \Sigma_{\mathbb{P}_k}(K)} \left( \int_K \nabla \varphi_{\sigma} \cdot \sum_{j \in [1, N_C]} (\alpha_j h_j) dx \right) = 0. \quad (2.119) $$

Provided that both conditions (2.78) and (2.119) are fulfilled, the Flux Reconstruction scheme (2.88) can therefore be properly embedded in the Residual Distribution framework, the conservation relation (2.100) being guaranteed for the associated residuals (2.112). We can therefore pursue the analysis and derive further conditions on the correction term to grasp physically relevant solutions.

### 2.5.3 Entropy conservation of the recast scheme

As introduced in Section 2.1.5 and considered in Section 2.4.3 while stating the properties of a Residual Distribution scheme, ensuring that the weak solution towards which the approximated ones converge as $h \to 0$ is physically relevant can be done by asking the entropy conditions (2.109a) to be fulfilled for any $h \in \mathbb{R}_+$ and a given entropy pair $(S, G)$. Those conditions reading here under the assumption that the boundary quadrature points verify $\kappa_G = \kappa_C$, any approximate solution $Y^h$ should therefore fulfil

$$ \sum_{\mathcal{G} \in \Sigma_{W(K)}} \langle U_{\mathcal{G}}, \Phi_{\mathcal{G}}^{K} \rangle \geq \int_{\partial K} G_n(U^K, U^{K-}) d\gamma \quad \forall K \in \tau^h \quad (2.120) $$

for $U = \nabla_Y S(Y)$ the entropy variable that is associated to the entropy pair $(S, G)$ and where $G_n$ is a given numerical entropy flux consistent with $G \cdot n$.

**Correcting the Flux Reconstruction scheme** In practice, enforcing the relation (2.120) reduces to modifying the residuals $\Phi_{\sigma}^{K,FR}$ defined in (2.115) so that the assumptions of the Proposition 2.4.1, Proposition 2.4.2, and Theorem 2.4.1 are fulfilled. There, a simple approach is to enforce the property by construction, considering a modification of the residuals that is similar to the one already performed in (2.118) for turning a Discontinuous Galerkin scheme into an Flux Reconstruction one. Namely, we add a correction terms $\tau_{\sigma}$ to each of the partial Flux Reconstruction residuals $\Phi_{\sigma}^{K,FR}$, modifying their local distribution so as to enforce the entropy constraint (2.120) while preserving the conservation relation (2.100). Hence, as it reduces in the end to constructing an entropy stable Flux Reconstruction scheme out of an initial Discontinuous Galerkin formulation by modifying its residuals through two correction terms, we express the entropy stable residuals $\Phi_{\sigma}^{CS}$ as

$$ \Phi_{\sigma}^{CS} := \Phi_{\sigma}^{K,FR} + \tau_{\sigma} = \Phi_{\sigma}^{K,DG} + r_{\sigma} + \tau_{\sigma}, \quad (2.121) $$
where the correction terms $r_σ$ and $τ_σ$ modulates the initial Discontinuous Galerkin approximation of the flux dynamics. In particular, the term $r_σ$ retrieves the cross-cell flux continuity, is of the form (2.117) and is already constrained by (2.119). A-contrario, the expression of $τ_σ$ is responsible for the entropy stability of the comprehensive numerical scheme and can be designed freely, provided that its representation in $V^C(K)$ makes the comprehensive characterisation $Σ_{V(K)}$ unisolvent for $V(K)$, that the entropy condition (2.120) is enforced, and that the conservation relations (2.100) are preserved. Yet, as we are not interested here in either of the properties separately, those two correction terms will be considered as a whole in below developments, easing the presentation of the construction.

**Achieving entropy stability** Let us first present a condition on the correction term $r_σ + τ_σ$ ensuring that the resulting scheme is entropy stable. To this aim, we observe that enforcing the entropy stability condition (2.120) can be done by modifying the Discontinuous Galerkin residuals from terms that depend on the local entropy error. Indeed, as here the proposed construction of the numerical scheme (2.121) relies on a corrective approach, the entropy stability can be obtained a-posteriori. It is therefore enough to compute the total entropy error produced on each primal element of the Discontinuous Galerkin scheme and to compensate it locally by correcting accordingly the value of the partial residuals involved in the approximation of the local dynamics.

Even though a wide range of possibilities exist for distributing relevant correction amounts to the partial residuals, one should keep in mind that the comprehensive definition of the corrected residuals should preserve the Flux Reconstruction nature of the resulting numerical scheme. Hence, as a first consequence, the entropy error should be evaluated on each primal element through the quadrature rule of the initial Flux Reconstruction scheme, i.e. by defining

$$E := \left( \int_{\partial K} G_n(U^K, U^{K-}) d\gamma - \sum_{\sigma \in K} \langle U_\sigma^K, \Phi_\sigma^{K, DG} \rangle \right) \quad (2.122)$$

for $G_n$ a given numerical entropy flux consistent on each boundary with $G \cdot n$. Equally crucial to the well-definition of the scheme, a second consequence relates to the conservation relations. Indeed, while the correction terms should preserve the conservation relations (2.100), we already know that the Flux Reconstruction nature of the numerical scheme asks the term $r_σ$ to fulfil (2.119). An analogous condition has therefore to hold as well for $τ_σ$ alone, i.e.

$$\sum_{\sigma \in \Sigma_{\mathbb{P}_k}(K)} r_{\sigma} = 0 \quad \text{and} \quad \sum_{\sigma \in \Sigma_{\mathbb{P}_k}(K)} \tau_{\sigma} = 0. \quad (2.123)$$

Assuming then that the correction term $τ_σ$ admits a formulation that makes the set of degrees of freedom $Σ_{V(K)}$ characterising the comprehensive flux ap-
proximation unisolvent for the space $V(K)$, one can use the fact that the mapping $Y \to U(Y)$ is bijective to create correction terms that are constructed from the value of the non-corrected solution itself. One can then particularly ask the correction terms to fulfil the relation

$$ \sum_{\sigma \in \Sigma_{\mathbb{F}_k}(K)} \langle U_{\sigma}, r_{\sigma} + \tau_{\sigma} \rangle \geq E, \qquad (2.124) $$

ensuring the entropy stability of the resulting Flux Reconstruction scheme. Indeed, given the definition (2.122) it then comes

$$ \sum_{\sigma \in \Sigma_{\mathbb{F}_k}(K)} \langle U_{\sigma}, \Phi_{\sigma}^{K,CS} \rangle - \int_{\partial K} G_n(U^K, U^{K-}) d\gamma = \sum_{\sigma \in \Sigma_{\mathbb{F}_k}(K)} \langle U_{\sigma}, r_{\sigma} + \tau_{\sigma} \rangle - E \geq 0. $$

The equation (2.120) is then fulfilled and we naturally get entropy conservation for our Flux Reconstruction scheme (2.121).

Let us now put the constraints that are inherent to the Flux Reconstruction structure aside and focus momentarily on the Residual Distribution framework alone. We there derive in the next paragraph an expression of the correction term $r_{\sigma} + \tau_{\sigma}$ for which the equality holds, enforcing the entropy conservation of the resulting Residual Distribution scheme. We will then be able to fall back on the Flux Reconstruction framework in the second next paragraph by adding one further assumption on the shape of the correction term $\tau_{\sigma}$, and thus provide a definition of an entropy stable Flux Reconstruction scheme both in its classical expression (2.88) and in its Residual Distribution expression (2.121).

**Achieving entropy conservation at the residuals level** In light of the developments of the previous paragraphs, asking a recast Flux Reconstruction scheme to be entropy conservative requires to find an expression of the correction terms $r_{\sigma}$ and $\tau_{\sigma}$ that fulfil

$$ \left\{ \sum_{\sigma \in \Sigma_{W(K)}} (\Phi_{\sigma}^{K,CG} + r_{\sigma} + \tau_{\sigma}) = \int_{\partial K} F_n(Y^K, Y^{K-}) d\gamma. \quad (2.125a) \right. $$

$$ \left. \sum_{\sigma \in \Sigma_{W(K)}} \langle U_{\sigma}, \Phi_{\sigma}^{K,CG} + r_{\sigma} + \tau_{\sigma} \rangle = \int_{\partial K} G_n(U^K, U^{K-}) d\gamma. \quad (2.125b) \right. $$

There, using the bijective mapping $U: Y \to \nabla_Y S(Y)$, considering the term $r_{\sigma}$ in the shape (2.117), and recalling our above developments, it reduces to asking

$$ \left\{ \sum_{\sigma \in \Sigma_{W(K)}} r_{\sigma} = 0 \quad \text{and} \quad \sum_{\sigma \in \Sigma_{W(K)}} \tau_{\sigma} = 0, \quad (2.126) \right. $$

$$ \left. \sum_{\sigma \in \Sigma_{W(K)}} \langle U_{\sigma}, r_{\sigma} + \tau_{\sigma} \rangle = E \quad (2.127) \right. $$
for E the approximate entropy error term given in (2.122). Hence, considering the problem (2.126)-(2.127) by itself, one can design an entropy stable Residual Distribution scheme (2.121) by defining the following correction terms.

**Proposition 2.5.1** Correction terms ensuring entropy conservation

The Residual Distribution scheme (2.121) defined from the Discontinuous Galerkin residuals (2.103) to which is added the term $r_{\sigma} + \tau_{\sigma}$ defined as

$$r_{\sigma} + \tau_{\sigma} = \beta(U_{\sigma} - \bar{U}), \qquad (2.128)$$

for any $\sigma \in \Sigma_{W(K)}$, where

$$\beta = \frac{E}{\sum_{\sigma \in \Sigma_{W(K)}} (U_{\sigma} - \bar{U})^2} \quad \text{and} \quad \bar{U} = \frac{1}{|\Sigma_{W(K)}|} \sum_{\sigma \in \Sigma_{W(K)}} U_{\sigma}, \qquad (2.129)$$

is both conservative and entropy conservative.

**Proof.** Let us simply retrieve the relations (2.126) and (2.127) given the expression (2.128) of the correction terms.

▷ Considering first the conservation relations, one shall just recall that the condition on $r_{\sigma}$ is obtained by the assumption (2.119) on the correction function (2.77). The conservation relation on the term $\tau_{\sigma}$ is then obtained through the conservation of the comprehensive correction term, i.e.

$$\sum_{\sigma \in \Sigma_{W(K)}} \tau_{\sigma} = \sum_{\sigma \in \Sigma_{W(K)}} r_{\sigma} + \sum_{\sigma \in \Sigma_{W(K)}} \tau_{\sigma} = \sum_{\sigma \in \Sigma_{W(K)}} (r_{\sigma} + \tau_{\sigma}) \\ = \sum_{\sigma \in \Sigma_{W(K)}} \beta(U_{\sigma} - \bar{U}) = \sum_{\sigma \in \Sigma_{W(K)}} \beta U_{\sigma} - \sum_{\sigma \in \Sigma_{W(K)}} \beta U_{\sigma} = 0.$$

▷ The entropy conservation relation also follows from a simple computation. Indeed, writing

$$\sum_{\sigma \in \Sigma_{W(K)}} \langle U_{\sigma}, r_{\sigma} + \tau \rangle = \sum_{\sigma \in \Sigma_{W(K)}} \langle U_{\sigma}, \beta(U_{\sigma} - \bar{U}) \rangle = \sum_{\sigma \in \Sigma_{W(K)}} \alpha(\langle U_{\sigma}, U_{\sigma} \rangle - \langle U_{\sigma}, \bar{U}_{\sigma} \rangle) \\ = \sum_{\sigma \in \Sigma_{W(K)}} \beta(\langle U_{\sigma}, U_{\sigma} \rangle - 2\langle U_{\sigma}, \bar{U}_{\sigma} \rangle + \langle \bar{U}_{\sigma}, \bar{U}_{\sigma} \rangle) \\ = \beta \sum_{\sigma \in \Sigma_{W(K)}} (U_{\sigma} - \bar{U})^2 = E,$$

the relation (2.127) is automatically obtained.

So defined, the corrected scheme (2.121) can also be classified as a Residual Distribution one, and is both conservative and entropy conservative in the discrete sense. Furthermore, the entropy error is controlled by an estimate that only depends on the spatial dimension, initial Flux Reconstruction discretisa-
tion order, and the quadrature scheme (see [52] for a complete derivation of the estimate). Namely,

$$E = - \sum_{\sigma \in K} \langle U_{\sigma}, \Phi_{\sigma} \rangle + \int_{\partial K} G_n(U^h, U^{h,-}) d\gamma = \mathcal{O}(h^{d+k+1}). \quad (2.130)$$

However, it is not clear yet that the correction terms preserves the Flux Reconstruction structure of the recast scheme. To ascertain it, let us draw a further assumption on the shape of the correction term $\tau_{\sigma}$ and build a system having for unknown the correction function (2.77), whose any solution necessarily yields entropy conservation for any Flux Reconstruction scheme of the form (2.88) and can be characterised from a set of degrees of freedom $\Sigma_C$ that ensures the characterisation $\Sigma_{V(K)}$ of comprehensive flux approximation to be unisolvent for $V(K)$.

### 2.5.4 Flux Reconstruction nature of the resulting entropy stable Residual Distribution scheme

Now that we could identify in Section 2.5.3 an expression of the correction terms so that the resulting Residual Distribution scheme is entropy conservative, let us ascertain the Flux Reconstruction nature of the resulting scheme. Concluding the definition of the presented entropy conservative Flux Reconstruction scheme, a system having for unknown the correction function is in particular derived in the process, whose solvability motivates the entire Chapter 3 and whose solutions provide admissible correction functions that guarantee entropy conservation of the Flux Reconstruction scheme (2.121) in its classical definition (2.88).

**Flux Reconstruction nature of the residuals** A natural strategy to ensure that the scheme (2.121) enters the scope of the Flux Reconstruction framework is to prove the existence of a correction terms reformulation that guarantees to the comprehensive residuals a Flux Reconstruction structure.

Following this approach and recalling all the assumptions that we already ask on the Discontinuous Galerkin partial residuals $\Psi_{\sigma}^{K,CG}$, a sufficient way to ensure the existence of such a reformulation is to construct a set of functions $\psi_j$ contained in $V^C(K)$ so that the correction terms $r_{\sigma} + \tau_{\sigma}$ can read $\int_K \varphi_{\sigma} \sum_{j \in [1, N_C]} \psi_j dx$ for a the original number and locations of the initial correction points $\varkappa_C$. Indeed, so defined, the parallel with the formulation (2.116) is clearly apparent, and the expression of the unifying correction function can be expressed as simply as

$$\psi := \sum_{\sigma \in [1, N_C]} \alpha_j \psi_j. \quad (2.131)$$
**Proposition 2.5.2** **Existence of correction functions**

There exists a set of correction functions {$\psi_j$}$_{j \in [1, N_C]} \in V^C(K)$ so that

$$
r_{\sigma} + \tau_{\sigma} = \int_{K} \nabla \varphi_{\sigma} \sum_{j \in [1, N_C]} \alpha_j \psi_j \, dx, \qquad (2.132)
$$

**Proof.** Let us first assume that the entropy correction term can be expressed similarly as the correction term defining non-corrected Flux Reconstruction scheme (2.116), that is, reads

$$
\tau_{\sigma} = \int_{K} \nabla \varphi_{\sigma} \cdot \sum_{j \in [1, N_C]} \alpha_j \vartheta_j \, dx \qquad (2.133)
$$

for given lifting functions $\vartheta_j$, $j \in [1, N_C]$ that live in the Raviart – Thomas space $RT_k(K)$. The total correction term then writing

$$
r_{\sigma} + \tau_{\sigma} = \int_{K} \nabla \varphi_{\sigma} \sum_{j \in [1, N_C]} \alpha_j (h_j + \vartheta_j) \, dx, \quad (2.134)
$$

one can define the comprehensive correction function $\psi_j = h_j + \vartheta_j$. Recall-
ing now the constraints (2.128) that are associated to the term $\tau_\sigma$ to ensure
the entropy conservation, and the Admissibility conditions 2.4.1 applying
on both the lifting functions $h_j$ and $\vartheta_j$ to preserve the consistency of the
scheme, we therefore require the expression (2.131), representing the com-
prehensive correction function, to satisfy the Neumann problem

$$
\int_K \nabla \varphi_{\sigma} \sum_{j \in [1, N_C]} \alpha_j (\psi_j) \, dx = \beta (U_{\sigma} - \bar{U}) \quad (2.135)
$$

$$
\psi_j(x_l) \cdot n|_{\partial K} = \delta_{jl}, \quad \forall j \in [1, N_C], \forall x_l \in \mathcal{N}_C, \quad (2.78)
$$

$$
(\psi_j \cdot n)|_f = 0, \quad \forall j \in [1, N_C], \forall f \in \partial K, \text{ s.t. } x_j \notin f. \quad (2.82)
$$

where we further ask $(\psi_j^{K,C} \cdot n)|_{\partial K} \in \{w \in L^2(\partial K), \forall f \in \partial K, w|_f \in \mathbb{P}_k(f)\}$ to ensure the unisolvence of the boundary characterisation from the boundary constraints (2.78)–(2.82) alone.

▷ In the case where the cell K is either a simplex or has a quadrilateral shape, a solution is guaranteed by considering each lifting function $\psi_j$ in the Raviart – Thomas space of order $k \ge \min\{N_C - 1, 1\}$. Indeed, there, the conditions (2.78) and (2.82) are already fulfilled by constructing the function $\psi_j$ as the nodal normal basis function associated to the correction point $x_j \in \mathcal{N}_C$, and the constraint (2.135) can be enforced by making use of the internal basis functions of the Raviart – Thomas space.
▷ In the case where the cell K is a general polytope, one however requires the existence of a specific discretisation space where the condition conditions (2.78) and (2.82) can be obtained by construction and where it remains enough freedom to adapt the lifting functions to the condition (2.135). A class of admissible discretisation space exists, whose construction will be detailed in the entire Chapter 3 and applied to any polytope, not necessarily convex. We can by example consider the particular space (3.37), in which the lifting functions can be expressed from the set of basis functions that share a similar behaviour on the boundary with the ones of the Raviart – Thomas discretisation.

Hence, for any polytopal cell K, the problem (2.78)–(2.82)–(2.135) admits a solution and a reformulation of the residual as a Flux Reconstruction correction term is possible.

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

Highlighted in the above proof, the definition of the correction space $V^C(K)$ is therefore crucial to ensure that a given Flux Reconstruction scheme can be corrected to become entropy conservative. In particular, as the definition of the comprehensive flux discretisation framework $(K, V(K), \Sigma_{V(K)})$ is directly impacted by the correction space and characterisation method used to determine the correction function, one should make sure that the overall discretisation $\Sigma_{V(K)}$ is unisovlent in $V(K)$, i.e. that the combination of the degrees of freedom characterising the features of the solutions that are described from the possibly non-empty subspace $V^C(K) \cap V^D(K)$ are not linearly dependent nor contradictory. In other words, even if so far the proposed scheme has a Flux Reconstruction structure, it belongs to the Flux Reconstruction framework only if the chosen characterisation $\Sigma_C$ of the correction function makes the overall characterisation $\Sigma_{V(K)}$ unisolvent in $V(K)$. If this is granted in the simplicial and quadrilateral for any characterisation of correction functions living in the Raviart – Thomas space, verifying Proposition 2.4.1, one has however to keep it in mind when selecting the characterisation of the space $H_k(K)$ for cells having higher number of faces.

**Construction of the correction function** More than rooting the proof of Proposition 2.5.2, solving the system (2.78)–(2.82)–(2.135) allows us to construct the general correction function

$$ \psi = \sum_{j \in [1, N_C]} \alpha_j \psi_j, \qquad (2.136) $$

and to define an entropy conservative Flux Reconstruction scheme in its classical expression (2.88), the consistency of the original Flux Reconstruction scheme being naturally preserved by the definition given to the terms $\alpha_j$.
### 2.5.5 Summary of the construction

To conclude the development of an entropy stable Flux Reconstruction scheme, let us overview its construction and highlight the necessary assumptions for its practical definition.

**Theorem 2.5.1** *Entropy stable Flux Reconstruction scheme*

We consider here the approximation of the operator $D_F^h$ associated to the semi-discrete formulation (2.61) of the problem (2.1) by the particular the Flux Reconstruction scheme given by the formulation

$$D_F^h(x_m) = -\nabla \cdot (F^{K,D})(x_m) - \nabla \cdot (F^{K,C})(x_m), \quad x_m \in \mathcal{X}_Y,$$

where the flux points $\mathcal{X}_F$ coincide with the solution points $\mathcal{X}_Y$ and are located within the cell $K$, and where the correction points $\mathcal{X}_C$ are distributed on the cell's boundary. The initial flux approximation $F^{K,D}$ is then given by (2.72), and the correction function $F^{K,C}$ is then chosen upon constraints ensuring to the Flux Reconstruction scheme both a Residual Distribution formulation and entropy conservation.

▷ In particular, assuming that the correction function takes the form

$$F^{K,C} := \psi = \sum_{j \in [1, N_C]} (\alpha_j \nabla \cdot h_j), \quad (2.77)$$

then the correction procedure preserves the consistency of the initial Flux Reconstruction scheme by defining the constant $\alpha_j$ as $\alpha_j = (\mathcal{F}^I(Y_j^K, Y_j^{K-}) \cdot n - F^{K,D}(Y_j^K))$ and requiring the lifting functions $h_j, j \in [1, N_C]$ to fulfil

$$(\psi_j^{K,C} \cdot n)|_{\partial K} \in \{w \in L^2(\partial K), \forall f \in \partial K, w|_f \in \mathbb{P}_k(f)\} \quad (2.74)$$

$$\psi_j(x_l) \cdot n|_{\partial K} = \delta_{jl}, \quad \forall j \in [1, N_C], \forall x_l \in \mathcal{X}_C, \quad (2.78)$$

$$(\psi_j \cdot n)|_f = 0, \quad \forall j \in [1, N_C], \forall f \in \partial K, \text{s.t. } x_j \notin f. \quad (2.82)$$

▷ In addition, if the lifting functions are chosen so that the correction term

$$\psi_{\sigma} = \int_{K} \nabla \varphi_{\sigma} \cdot \psi \, dx \\ - \sum_{\sigma \in \Sigma_{\mathbb{P}_k(K)}} \psi_{\sigma} = 0, \quad (2.119)$$

for a given quadrature rule $\int$, any $\sigma \in \Sigma_{W(K)}$ degree of freedom, and where $\varphi_{\sigma}$ the corresponding dual basis function which are contributing to the definition of the element $(K, W(K), \Sigma_{W(K)})$ that is associated to the recast formulation (2.116) in $W(K)$, then the Flux Reconstruction scheme (2.88) is conservative and can be recast into the Residual Distribution framework, all the results of Section 2.4.3 therefore applying.
Lastly, if it further holds $\sum_{\sigma \in K} \langle U_\sigma, \psi_\sigma \rangle \ge E$, for E the entropy error (2.122) of the Discontinuous Galerkin scheme, then the Flux Reconstruction scheme (2.88) is additionally entropy stable.

To conclude this summary by an example, let us point out that as any correction function solving the discrete Neumann problem (2.78)–(2.82)–(2.135) satisfies all of the above criterion, any such solution necessarily yields an entropy conservative Flux Reconstruction scheme.

### 2.5.6 Conclusion and perspectives

Initially designed in the Residual Distribution framework, the entropy stable Flux Reconstruction scheme developed here is built in a corrective approach by constructing an entropy stabilisation term ensuring physically relevant and conservative solutions. Benefiting further from the shape flexibility of the Residual Distribution, the developed scheme can be used on any hybrid non-convex polygonal mesh, either in its residual expression (2.121) or in its direct formulation (2.88). Thus extending the results of [56, 73, 74, 75] to general polygonal grids, it is however subject to the existence of solution to the problem (2.78)–(2.82)–(2.135) for any polytope.

This problem requires in particular the existence of a finite dimensional variational space that can be partially characterised by $n$ sets of pointwise degrees of freedom that are supported on the face $f_j \in \partial K, j \in [1, n]$, and whose leftover freedom is enough for the correction function to be tuned against the inner entropy conservation constraint (2.128). In addition, as we allow the use of hybrid meshes, one should ensure that the nature of the boundary discretisation does not depend on the number of faces of the element. Hence, the classical Raviart–Thomas space used previously in [56], cannot be considered as such for elements having more than $4d$ faces or when considering the lowest flux discretisation order.

We therefore developed a new discretisation framework whose construction is detailed step by step in Chapter 3, and of which a more complete investigation can be found in [63]. There, the boundary characterisation is independent from the inner characterisation, providing the possibility to tune more finely and face-independently the coordinate-wise correction functions. It further ensures that the property (2.78) can be achieved for any number of discretisation points while letting the inner behaviour adapting towards (2.135). Lastly, the normal component of any function belonging to any developed space being a polynomial whose order is independent of the cell's number of faces, a nodal characterisation ensuring Proposition 2.4.1 to hold and the plain Flux Reconstruction formulation of the presented scheme (2.88) to be valid on polygonal shape can naturally be defined.
<figure><img src="image_33.png" /></figure>
<figure><img src="image_34.png" /></figure>
$H(\operatorname{div})$-conforming elements on
general polytopal meshes

Motivated by developing conservative and stable numerical schemes in the
scope of the Residual Distribution framework, we aim here to establish the
theoretical foundations ensuring the existence and expressibility of correction
functions fulfilling the relations (2.74), (2.78), and (2.82).

In that respect, the pursued strategy is to express the correction functions
directly from their gradients and to discretise them in a finite dimensional
vector space whose properties enforce the desired boundary behaviour. Spaces
that are here suitable typically provide a low order discretisation of the $H(\operatorname{div})$
Sobolev space, ensuring in particular a robust representation of the discretised functions' normal component. When coupled with a given characterisation of the discretised functions, their definitions give rise to a range of
*so-called* $H(\operatorname{div})$-conforming elements that are classically used *e.g.* in the context of mixed and hybrid Finite Elements methods. Also providing a simple
representation framework, they therefore appear as a relevant scope for our
correction functions. However, those spaces are generally defined for specific
geometries and have a rigid variational structure, while the entropy and stable schemes introduced in Section 2.5 are to be set up on any hybrid mesh
containing different cell types, requiring the correction function to exist and
be rigorously expressible independently from the local geometry. We therefore
developed a class of $H(\operatorname{div})$-conforming spaces and elements that are geometrically flexible, obtaining as a by-product an improved variational flexibility
that preserves the classical boundary properties of the discretised functions' normal component.

After briefly describing the $H(\operatorname{div})$-conforming notions in Section 3.1, we
describe in Section 3.2 two classical elements from which stem our investigations. We then introduce the developed geometrically and variationally flexible construction in Section 3.3 providing in particular a rigorous framework
for Flux Reconstruction schemes on arbitrary polytopes.
Chapter contents

<table><tbody><tr><td>3.1</td><td>Conforming discretisation strategies</td><td>103</td></tr><tr><td>3.1.1</td><td>Discrete features and conformity</td><td>104</td></tr><tr><td>3.1.2</td><td>The variational space H(div) and its specificities</td><td>105</td></tr><tr><td>3.1.3</td><td>Application of H(div)-conforming discretisations</td><td>106</td></tr><tr><td>3.1.4</td><td>Classical polytopal and H(div)-conforming elements</td><td>108</td></tr><tr><td>3.1.5</td><td>Scope of the presented method</td><td>109</td></tr><tr><td>3.2</td><td>Description of two conforming elements</td><td>112</td></tr><tr><td>3.2.1</td><td>Simplicial Raviart-Thomas elements</td><td>112</td></tr><tr><td>3.2.2</td><td>Quadrilateral Raviart-Thomas elements</td><td>120</td></tr><tr><td>3.2.3</td><td>The classical Virtual Elements</td><td>122</td></tr><tr><td>3.2.4</td><td>The H(div)-conforming virtual element space</td><td>124</td></tr><tr><td>3.3</td><td>A framework for arbitrary polytopes</td><td>127</td></tr><tr><td>3.3.1</td><td>The general idea in a nutshell</td><td>128</td></tr><tr><td>3.3.2</td><td>Necessary conditions over spaces and elements</td><td>128</td></tr><tr><td>3.3.3</td><td>A class of admissible approximation spaces</td><td>129</td></tr><tr><td>3.3.4</td><td>Definition of H(div)-conforming elements</td><td>135</td></tr><tr><td>3.3.5</td><td>Summary of the construction</td><td>153</td></tr><tr><td>3.3.6</td><td>Two examples in two dimensions</td><td>155</td></tr><tr><td>3.3.7</td><td>Link to other discretisation spaces</td><td>158</td></tr><tr><td>3.3.8</td><td>Numerical results</td><td>161</td></tr><tr><td>3.3.9</td><td>Critical assessment of the proposed approach</td><td>165</td></tr></tbody></table>
## 3.1 Conforming discretisation strategies

Enforcing known-properties of the dynamics or *a-priori* known features of the solution directly in the discretisation setting strengthens the robustness of the approximation and is therefore particularly popular in modelling physical behaviour as *e.g.* optics, electrostatic and fluid dynamics.

Interested here in solving hyperbolic problems with particular applications in fluid dynamics, we first present some classical ways that could use the $H(\text{div})$-conforming elements introduced in this chapter to enforce particular properties. We therefore go slightly beyond our immediate goal of laying a theoretical framework for the correction function (2.88), aiming to provide a quick overview of the classical scope of the presented work. To that end, we consider our initial value problem

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial Y}{\partial t} + \nabla \cdot \mathcal{F}(Y) = 0 \\ Y|_{\partial \Omega} = g, \end{array} \right. \qquad (3.1) $$

over a domain $\Omega \subset \mathbb{R}^d$ and where $Y \in \mathbb{R}^q$, $q \in \mathbb{N}$, is the time-space dependent state vector known only at the initial time. Behaviourally speaking, the flux $\mathcal{F}: \mathcal{D} \to \mathbb{R}^q$ describes the system dynamics over a domain $\mathcal{D}$ representing the subset of $\mathbb{R}^q$ for which the state values $Y$ have a physical meaning. Complementary, the function $g: \mathbb{R}_+ \times \partial\Omega \to \mathbb{R}$ prescribes the boundary conditions, exemplified here as being of Dirichlet type.

As seen in the preliminary Section 2.3.1, a typical approximation of this problem is obtained by deriving a weak formulation, whose expression is later split on small cells $K$ forming a mesh $\tau^h$ that approximates the domain $\Omega$ and defines the computational domain as $\Omega = \cup_{K \in \tau^h} \{K\}$. The weak formulation is then discretised by choosing local finite dimensional trial $U$ and test $V$ spaces within which the solution and the local test functions are respectively represented. Choosing further a unisolvent set of degrees of freedom $\{\sigma\}$ that characterises the discretised solution in the trial space $U$, the set $(K, U, \{\sigma\})$ defines a so called *element*, and a finite element formulation can be derived. Writing in a semi-discrete formulation as

$$ \frac{d}{dt} \left( \int_K Y \, dx \right) + \int_K \nabla \varphi \cdot \mathcal{F}(Y^h) \, dx + \int_{\partial K} g \varphi \, dx = 0, \quad \forall K \in \tau^h, \forall \varphi \in V, $$

any consistent flux approximation would then define a numerical scheme. Hence, in such discretisations the definition of the trial space predetermines the type of features that can be represented at the discrete level, while the chosen degrees of freedom prescribe the actual representation of the discrete solution and can therefore enforce the known features directly in its structure.
### 3.1.1 Discrete features and conformity

Features that are enforced at the discrete level typically relates to the global regularity of the discrete solution, that is, the regularity of

$$Y: \Omega \to \mathbb{R}^q \\ x \mapsto Y^K(x), \quad K \text{ s.t. } x \in K$$

given the local reconstructions $Y^K: K \to \mathbb{R}^q, K \in \tau^h$. The ability of a compactly and locally defined element to transfer its properties to the global level is subsumed by the notion of conformity.

**Definition 3.1.1 Conforming elements**

Let $\Omega$ be a computational domain belonging to $\mathbb{R}^d, d \in \mathbb{N}$ that is tessellated by a mesh $\tau^h = \{K\}$, and define

* $W(\Omega)$ a variational Sobolev vector space whose functions are defined on $\Omega$ and take values in $\mathbb{R}^q, q \in \mathbb{N}$,

* $V$ a finite dimensional subspace contained in $W(\Omega)$,

* $E_K := (K, V, \{\sigma\})$ an element defined on the cell $K$.

Let also $\mathcal{E}_K$ be the space of discrete functions $p: K \to \mathbb{R}^q$ that are contained in $V$ and represented by the values of $\{\sigma\}$, and consider the global approximation $V^h = \{v: \Omega \to \mathbb{R}^q \text{ s.t. } \forall K \in \tau^h, v|_K \in \mathcal{E}_K\}$. Then, if

$$V^h \subset W(\Omega),$$

the element $(K, V, \{\sigma\})$ is said to be $W$-conforming.

In essence, the conforming property only states that when the definition of $V$ confers a certain regularity to the local solution reconstructions $Y_K$, the set of selected degrees of freedom $\{\sigma\}$ ensures that the global reconstruction $Y$ benefits from the same regularity on the full computational domain $\Omega$. When locally seen, the property therefore translates to the local boundary treatments being compatible with the asked regularity on a global scale.

The chosen regularity thus give rise to different classes of conforming elements. In particular, asking continuity over the entire domain for each variable typically yields Continuous Galerkin approaches, and asking for higher regularity may yield to using $H^1$-conforming techniques. More elaborated definitions of regularities can also be asked and are typically related to physical principles, possibly involving interactions between variables. The $H(\operatorname{curl})$ and $H(\operatorname{div})$ conforming elements enter naturally in this scope, as their definition ensure that no jump will occur in the curl or in the divergence of the discretised quantity across the element faces, respectively. Let us now focus on the particular scope of the presented work, defining the $H(\operatorname{div})$ space and exemplifying possible $H(\operatorname{div})$-conforming discretisations and practical use.
### 3.1.2 The variational space $H(\text{div})$ and its specificities

Before going to its application, let us specify the definition of the variational space $H(\text{div})$ and point out a geometrical particularity of $H(\text{div})$–conforming elements that one has to be aware of in view of designing numerical schemes.

**The variational space $H(\text{div})$** In order to meet the expectations of the previous paragraph, one requires $H(\text{div})$ to be a variational space that provides regularity to its function’s normal components and from which finite dimensional subspaces and conforming finite elements can be defined. A natural idea is thus to choose a Sobolev vectorial space and ask for a $L^2$-regularity.

**Definition 3.1.2** *$H(\text{div})$ space*

Let $K$ be a general simplicial reference shape contained in $\mathbb{R}^d$, where $d \in \mathbb{N}$.
It is set

$$H(\text{div}, K) = \{v \in (L^2(K))^d; \nabla \cdot v \in L^2(K)\}.$$

As $H(\text{div})$ is also an Hilbert space, we can endow it with the following norm.

$$\forall v \in H(\text{div}, K), \quad \|v\|_{\nabla, K}^2 = \|v\|_{L^2(K)}^2 + \|\nabla \cdot v\|_{L^2(K)}^2.$$

So defined, the space $H(\text{div}, K)$ imposes regularity of the functions’ normal component within the cell $K$ and on its faces $f \in \partial K$. However, given any finite dimensional subspace of $H(\text{div})$, the continuity of the discrete functions’ normal component across the faces is not guaranteed by the sole definition of the space, but requires a proper definition of degrees of freedom for their representation to be identical on the faces regardless the element the functions have been discretised in. This is particularly delicate when dealing with cells having various number of faces, and will be the topic of the entire Section 3.3.

**Geometrical specificity of $H(\text{div})$–conforming elements** Let us now assume that for any polytopal shape $K \subset \mathbb{R}^d$, $d \in \mathbb{N}$, the triplet $(K, V, \{\sigma\})$ defines a $H(\text{div})$–conforming element, that is, $V$ is a finite dimensional subspace of $H(\text{div}, K)$ and the set $\{\sigma\}$ enforces the $L^2$-regularity of the normal components across the faces. Let us also assume that for the sake of computational handiness, each polytopal element can be identified to a reference one.

Then, the definition of the geometrical transformation from a polytope $K$ to its reference shape is not straightforward. Indeed, usual transformations would prevent the flux of vector fields through the faces to be conserved, thus modifying the local degrees of freedom supported on the boundary and breaking the $H(\text{div}, K)$–conformity. Specific mappings have therefore to be employed, as e.g. the Piola transformation for simplicial and quadrilateral shapes (see Section 2.2.4). However, such transformations are extremely delicate to define in a general setting, especially in the non-convex case where the element orientation itself is questioned.
### 3.1.3 Application of $H(\text{div})$-conforming discretisations

Having defined the variational $H(\text{div})$ space and explained the role of conforming settings, let us now present two classical uses of $H(\text{div})$-conforming elements, one for discretising the local solution, the other the flux itself.

**Mixed finite elements method** Initially though for elliptic problems, the mixed finite elements method [77, 78, 79] has been generalised to applications where divergence free properties are encountered or when extra divergence regularity is required. Hence popular for solving problems involving the Hermolz, magnetohydrodynamics or Stokes equations [80, 81, 82, 83, 84], it is more generally considered for approximating solutions to physical models involving incompressibility conditions or direct divergence constraints.

Monitoring and enforcing the gradient behaviour directly in the discretisation of the variables that are impacted by the divergence conditions, the mixed elements method approximates simultaneously scalar and vector variables. More precisely, the vector variables are defined as subsets of the state variables that are ruled by divergence-based relations, as e.g. the velocity vector. The scalar variables can be either simply defined as those that are not ruled by their divergence in the dynamics, or as extra variables resulting from a system augmentation. In both cases, they define Lagrange multipliers and are considered as dual variables. A $H(\text{div})$-conforming discretisation is then applied to the vector variables, while a standard discretisation defined over a finite dimensional subspace of $L^2(K)$ is used for the scalar variables.

**Example**

For the sake of clarity, let us present an example of the mixed finite elements method applied to the incompressible Euler equations as presented in [85]. There, the derivation of the mixed formulation is based on the velocity-pressure formulation

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{\partial u}{\partial t} + u \cdot \nabla u + \nabla p = 0 \\[1em]
\displaystyle \nabla \cdot u = 0 \\[1em]
\displaystyle u \cdot n_{|\partial\Omega} = 0,
\end{array}
\right. 
\qquad (3.2) $$

where $u \in \mathbb{R}^2$ and $p \in \mathbb{R}$ are respectively the velocity vector and pressure variables, sought with regularity $u \in H(\text{div}, \Omega)$ and $p \in L^2(\Omega)$. Denoting $\otimes$ the dyadic product and recalling that for any vector $v_1$ and $v_2$ in $\mathbb{R}^2$, $\nabla \cdot (v_1 \otimes v_2) = v_2 \cdot \nabla v_1 + \nabla \cdot (v_2) v_1$, the relations (3.2) then also reads

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{\partial u}{\partial t} + \nabla \cdot (u \otimes u) + \nabla p = 0 \\[1em]
\displaystyle \nabla \cdot u = 0 \\[1em]
\displaystyle u \cdot n_{|\partial\Omega} = 0.
\end{array}
\right. 
\qquad (3.3) $$
Meshing then the domain $\Omega$ with a covering of $\tau^h = \{K\}$ of cells $K$, and choosing two finite dimensional test spaces $V \subset H(\text{div}, K)$ and $W \subset L^2(K)$ on each cell, one can then derive the cell-wise discrete weak formulation

$$ \left\{ \begin{array}{ll} \displaystyle \int_K \frac{\partial u}{\partial t} \cdot v - \int_K (u \otimes u) \tilde{\nabla} v - \int_K p \nabla \cdot v + \int_{\partial K} F(u, p) n \cdot v = 0, & \forall v \in V, \quad (3.4a) \\[1.5em] \displaystyle \int_K w \nabla \cdot u = 0, & \forall w \in W, \quad (3.4b) \end{array} \right. $$

where the variable $p$ can be seen as the Lagrange multiplayer associated to the incompressibility constraint, where the flux $F$ is chosen consistently with $F(u, p) = u \otimes u + pI$, and where $\tilde{\nabla}$ represents the broken gradient. One can then straightforwardly discretise by choosing the trial space of $u$ to be a finite dimensional subset of $H(\text{div}, \Omega)$. The relation (3.4b) can then be expressed directly at the discretisation level whenever definition of the degrees of freedom yields a conforming element and allows a direct control of the velocity vector flux at the element boundaries. Strictly speaking, (3.4b) is exactly retrieved whenever a $H(\text{div})$-conforming element is used.

Lastly, let us point out that among other techniques yielding a mixed finite element technique, the flux-vector reformulation [86, 87] seems promising and allows entropy considerations alongside with a conforming flux discretisation.

**Flux discretisation** Another application of the $H(\text{div})$-conforming theory lies in flux discretisation techniques and is thus particularly suitable for the approximation of solutions to hyperbolic systems. Indeed, there, not only the solution should be discretised, but one should also provide a consistent and physically relevant discretisation of the flux across the cell boundaries, as readily visible when applying the Stokes theorem on (3.2) to obtain

$$ \frac{d}{dt} \int_K Y \, dx + \int_{\partial K} F(Y) \cdot n \, dx + \int_{\partial K} \varphi, g \, dx = 0 \quad \forall K \in \{\tau^h\}. $$

In particular, by discretising the flux in a finite dimensional subspace $V$ of $H(\text{div}, K)$ and selecting a set of degrees of freedom $\{\sigma\}$ that makes the triplet $(K, V, \{\sigma\})$ conforming, a global regularity of the flux's normal component can be achieved and the pointwise Rankine-Hugoniot conditions directly enforced. More generally, it is possible to control the shape and characterisation of any discretised function's normal component by adjusting the definition of the degrees of freedom, and thus to enforce specific physical or numerical properties to the discrete flux. Our immediate goal presented in Section 2.5.6 therefore typically enters this scope, as we aim the discrete flux function to ensure the specific properties of the numerical scheme presented earlier in Section 2.3.2.
### 3.1.4 Classical polytopal and $H(\text{div})$-conforming elements

Solving a problem by using the techniques succinctly mentioned in Section 3.1.3 requires to define a $H(\text{div}, K)$-conforming element for each geometry that exists in the spatial discretisation $\tau^h$. In particular, when using polytopal hybrid meshes one should be able to define a local $H(\text{div}, K)$ element that provides a same discretisation nature and quality regardless the element geometry. Very delicate to obtain in practice, especially when allowing the use of non-convex cells, we aim here to highlight the design challenges by overviewing the main existing strategies, observing that classical constructions of $H(\text{div})$-conforming elements are specific to given geometries while geometrically flexible elements are usually non-conforming.

**Geometry flexible elements** Although less common than simplicial and quadrilateral elements, polytopal ones have grown interest and are now employed in a wide range of numerical schemes, from Mimetic Finite Differences [88] to polygonal Finite Elements [89, 90, 91] and Hybrid High order methods [92]. However, geometry flexible definitions usually make use of generalized barycentric coordinates or relies on projectors that echoes the investigated problem. They are therefore delicate to handle in distorted non-convex elements and lack variational flexibility. A more recent approach introduced in [7] under the name of *Virtual Elements Method* rather defines the approximation spaces through solutions to Poisson problems, improving both the variational and geometrical flexibility. Not relying on barycentric coordinates and allowing customisable projectors, the proposed elements are however scalar and not $H(\text{div})$-conforming.

**$H(\text{div})$-conforming elements** Mostly initiated by Raviart and Thomas [93], the definition of $H(\text{div})$-conforming elements usually relies on the particular structure of polynomial spaces and are hence difficult to define beyond simplicial and quadrilateral shapes, where the dimension of polynomial spaces connects with the number of faces. Their variational definitions can however be extended, as shown by the works of Nédélec [94] and Brezzi, Douglas and Marini [95] in the context of Mixed Finite Elements method. More recently, a mixed Petrov-Galerkin scheme using Raviart–Thomas elements has also been investigated [96], though, up to the authors knowledge, the proposed elements are still limited to simplicial and quadrangular shapes.

**Geometry flexible and conforming elements** The difficulty in defining $H(\text{div})$-conforming elements on general polytopes lies both in the dimension of the discretisation space, that should be adaptive to the number of faces, and in the control of the discretised quantities' divergence, as the degrees of free-
dom should be defined so that the global regularity of the discrete divergence
and normal component is ensured.

To overcome the latter issue, a first polygonal $H(\operatorname{div})$-conforming element
has been proposed in [97] using gradient reconstruction and pyramidal sub-
meshes tessellation. However, its construction requires some shape regularity
within the mesh. Not requiring any sub-mesh tessellation or further geomet-
rical constraint, the approach proposed in [98] then allowed to build minimal
degree conforming elements through Whitney forms, though being limited to
convex geometries. An enhanced flexibility on the element shape has been
obtained in [99], but its variational definition remains rigid and requires sta-
bilisation techniques. A first promising $H(\operatorname{div})$-conforming discretisation that
can be defined on any polytopal shape has been recently proposed in [100]
as an extension of the Virtual Elements method. There, the conformity is
enforced directly on the normal components, not leaving any coordinate-wise
freedom to the discretised function on the boundary. Therefore, the boundary
setting may appear quite restrictive, especially in applications for which a
direct component-wise characterisation or tangential information is suitable.

**3.1.5 Scope of the presented method**

In this chapter, we propose to construct a $H(\operatorname{div}, K)$-conforming element that
can be set up on any non-degenerated polytopal shape $K \subset \mathbb{R}^d, d \in \mathbb{N}$, even
non-convex, providing the same discretisation nature and quality regardless
the element geometry.

**Background choices and starting point** Primarily aiming to provide a
theoretical framework for the correction functions (2.77), we focused on dis-
cretisation techniques that would ensure the properties (2.78), (2.82), and
(2.135) to any discretised quantity. Hence, the $H(\operatorname{div})$-conforming approaches
appears natural, being able to control the discrete vector fields across the
cells' faces. Within this scope, the lowest order Raviart – Thomas elements
are of particular interest, as the degrees of freedom that are supported on
the boundary ensure a constant normal component of discretised quantities
on the boundary, making the realisation of the constraint (2.135) straightfor-
ward. However, while the variational setting of the Raviart – Thomas ele-
ments suits our constrained problem, its lack of flexibility towards the geome-
try makes it not directly applicable in the general scope of Residual Distribu-
tion schemes.

Indeed, while providing theoretical foundations to building entropy stable
and conservative schemes, we also aim to preserve the geometrical flexibility
conferred by the Residual Distribution framework. The use of finite dimen-
sional polynomial spaces of discretisation is therefore excluded, their dimen-
sions not being adaptive to the number of the cell's faces. In this respect, and as we do not require the explicit knowledge of the correction functions within the element, the Virtual Elements framework seems relevant. However, its conforming version do not benefit from a flexible variational framework, and the lowest order elements do not guarantee a constant value of the normal components of discretised quantities. We therefore focus on creating $H(\mathrm{div})$-conforming finite elements enhancing a boundary characterisation, built on a space whose construction combines the outer structure of the Raviart–Thomas space with the definition of the scalar Virtual Element space.

**Objectives and achievements** We propose here a construction that inherits from the interface properties of the Raviart–Thomas elements and benefits from the shape flexibility of the Virtual Elements discretisation. In particular, the developed elements can be built on any polytopal shape, even non convex, and are available in any dimension. Equally offering a flexible variational representation of the discretised quantity, they provide more than a theoretical framework on which the correction functions (2.77) can be decomposed by defining a new element class that can be used as such in the construction of further numerical schemes. Furthermore, used in this latter context one would only need the knowledge of the degrees of freedom, making the proposed construction also close in the spirit to the Virtual Elements framework.

**Outline** In order to ease the presentation of the proposed construction and observe more in details the intricate specificities encountered in designing finite dimensional spaces of discretisation and elements, we start in Section 3.2 by describing the two fundamental works from which is based our construction. We then briefly outline our developments in Section 3.3.1, before introducing our new class of discretisation spaces in Sections 3.3.2 and 3.3.3. In a third time, we detail possible definitions of $H(\mathrm{div})$-conforming elements in Sections 3.3.6 and 3.3.5, and compare them to an existing $H(\mathrm{div})$-conforming framework in Section 3.3.7. Finally, we test the presented elements through the behaviour of their corresponding basis functions in the numerical results Section 3.3.8. For interested readers, more details on the construction and projectors test can be found in the extended technical report [63].
### 3.2 Description of two conforming elements

As the proposed discretisation framework combines the definitions of both the Raviart – Thomas and Virtual elements, we first describe their respective constructions and detail their fundamental properties. Starting with the simplicial Raviart – Thomas elements in Section 3.2.1 to emphasise the challenges involved in defining discretisation spaces, we present its extension to the quadrilateral elements in Section 3.2.2. We then make a step towards flexible geometries in Section 3.2.3 with the Virtual Element approach.

#### 3.2.1 Simplicial Raviart – Thomas elements

One of the pioneering works investigating $H(\operatorname{div})$-conforming elements has been conducted by Raviart and Thomas [93] and led to the definition of simplicial and quadrilateral elements. Defining the finite dimensional space of discretisation as a vectorial polynomial subspace of $H(\operatorname{div})$ by exploiting directly the structure of scalar polynomial spaces, the conformity is enforced by a specific polynomial projection of the derivatives. We describe here the construction of the $d$-dimensional simplicial elements, built from the scalar space $\mathbb{P}_k(\mathbb{R}^d)$ of $k^{th}$-order polynomials. Presenting them from their origins to their nowadays' expression, we aim to highlight the fundamental notions to consider and difficulties to face when designing tailored discretisation spaces. We then describe the chosen set of degrees of freedom, concluding the definition of the Raviart – Thomas elements.

##### 3.2.1.1 Historical framework of Nédélec

Initially introduced by Raviart and Thomas in 1977, the Raviart – Thomas space and elements have been popularised by Nédélec's 1980 work [94] where the origins of the construction have been shown and a practical formulation derived. After constructing the Raviart – Thomas space along those lines for observing its essence, we outline its modern formulation and present recent proofs of its crucial properties by following [101]. For a deeper view of technicalities and detailed examples, we refer the interested reader to [102].

**Discretisation structure** Let us consider a simplex $K \subset \mathbb{R}^d$ where to build a $H(\operatorname{div}, K)$-conforming discretisation. The main idea of the Raviart–Thomas framework is to build a finite dimensional subspace $V$ contained in $H(\operatorname{div})$ in which the functions belong to $L^2(K)$ and see any linear combination of their $k^{th}$ order derivatives vanish for a given $k \in \mathbb{N}$. In practice, the discretisation space is defined as the kernel of a symmetrised linear differential operator of order $k$ applied to a vectorial polynomial space, the order driving the overall approximation quality offered by the discretisation space.
**Symmetrised differential operator** Aiming to specify the differential operator used in the Nédélec's definition [94] of the Raviart–Thomas space, let us introduce a vector valued function $u: \mathbb{R}^d \to \mathbb{R}^d \in L^2(K)$ whose components write $u = (u_1, \dots, u_d)^T$. There, one can obtain the disappearance of any linear combination of each component's $k^{th}$ order derivative by asking the $k^{th}$ order directional derivative of $u$ to vanish along any set of $(k+1)$ directional vectors $\chi = (\chi_1, \dots, \chi_{k+1}) \in \bigtimes_{l=1}^{k+1} \mathbb{R}^d$. In view of specifying the associated operator, let us denote the projection of the vectorial function $u$ on any directional vector $\chi_1 = (\chi_1^1, \dots, \chi_1^d)^T \in \mathbb{R}^d$ by $u \cdot \chi_1 = u_1\chi_1^1 + \dots + u_d\chi_1^d \in \mathbb{R}$ and define its $k^{th}$ order directional derivative along the given directional vectors through the $(k+1)$-differential form

$$ D^k u(\chi_1, \dots, \chi_{k+1}) : K \to \mathbb{R} \\ x \mapsto \begin{cases} u(x) \cdot \chi_1 & \text{if } k=0 \\ \bigcirc_{j=2}^{k+1} \left( \nabla(u(x) \cdot \chi_1) \cdot \chi_j \right) & \text{else,} \end{cases} \qquad (3.5) $$

where $\bigcirc_{j=2}^{k+1}$ denotes the composition operator

$$ \bigcirc_{j=2}^{k+1} : L^2(K) \to \mathbb{R} \\ u \mapsto \nabla(\nabla(\dots(\nabla(u \cdot \chi_1) \cdot \chi_2)\dots)) \cdot \chi_k \cdot \chi_{k+1}. $$

The differential operator $D^k$ is then symmetrised in $\mathbb{R}^d$ for considering the behaviour of each partial derivative along all the given directions. More specifically, writing $\zeta_{k+1}$ the set of all permutations of $(1, \dots, k+1)$, the symmetrised differential operator as initially considered by Nédélec writes

$$ \varepsilon^k: \mathbb{R}^d \times \dots \times \mathbb{R}^d \to \mathbb{R} \\ (\chi_1, \dots, \chi_{k+1}) \mapsto \frac{1}{(k+1)!} \sum_{\zeta \in \zeta_{k+1}} D^k u(\chi_{\zeta(1)}, \dots, \chi_{\zeta(k+1)}). \qquad (3.6) $$

Going further in the analysis by remembering that the function $u$ lies in $L^2(K)$, one can reduce the above expression by using the commutativity of the partial differential operators. Introducing then the vector of the $d$ Fréchet derivatives

$$ d^k u(\chi_2, \dots, \chi_{k+1}) = \times_{l=1}^{d} \left( \bigcirc_{j=2}^{k+1} \left( \frac{\partial}{\partial x} u_l \right) \cdot \chi_j \right), \qquad (3.7) $$

one can reduce the expression (3.6) to the following definition [101].

**Definition 3.2.1** *Symmetrised differential operator*

$$ \varepsilon^k: \mathbb{R}^d \times \dots \times \mathbb{R}^d \to \mathbb{R} \\ (\chi_1, \dots, \chi_{k+1}) \mapsto \frac{1}{k+1} \sum_{j=1}^{k+1} \chi_j \cdot d^k u(\chi_1, \dots, \chi_{j-1}, \chi_{j+1}, \dots, \chi_{k+1}) $$
### 3.2.1.2 Definition of the simplicial Raviart – Thomas space

Having specified the symmetrised differential operator $\varepsilon^k$ in Definition 3.2.1, we can now define an approximation space by selecting an application domain and expressing the corresponding kernel.

**Abstract definition of the space** The application domain coinciding with the space where the discretised functions lie, we consider the finite dimensional polynomial space $(\mathbb{P}_k(K))^d$. There, to further enforce $H(\text{div})$ – conformity in the discretisation, we restrict the space $(\mathbb{P}_k(K))^d$ to a subspace $V \subset H(\text{div}, K)$ where the $k^{\text{th}}$ order derivative of the functions vanish, ensuring that $\nabla \cdot u \in L^2(K)$. With respect to the above developments, we then get the definition of the simplicial Raviart – Thomas space.

**Definition 3.2.2** *Simplicial Raviart – Thomas space*

$$RT(K) = \{u \in (\mathbb{P}_k(K))^d, \varepsilon^k u = 0\}$$

This early definition immediately shows the implications $V \subset (\mathbb{P}_k(K))^d$ and $\mathbb{P}_{k-1}(K) \subset V$, but is not convenient in practice. We therefore explicit the kernel of $\varepsilon^k$ and present the classical definition of the Raviart – Thomas space, relying on the kernel’s parametrisation.

**Space characterisation** Applied to the polynomial space $(\mathbb{P}_k(K))^d$, the differential operator $\varepsilon^k$ has a kernel that can be expressed by the polynomial’s derivatives, yielding the following parametrisation.

**Proposition 3.2.1** *Kernel of the symmetrised differential operator*

$$\text{Ker}(\varepsilon^k) = (\mathbb{P}_{k-1}(K))^d \oplus \{u \in (\mathbb{P}_{[k]}(K))^d \text{ s.t. } \sum_{i=1}^{d} x_i u_i = 0\}$$

In order to prove this kernel parametrisation, we look for a characterisation of the functions $u \in (\mathbb{P}_k(K))^d$ satisfying $\varepsilon^k u(\chi_1, \cdots, \chi_{k+1}) = 0$ for any set of directions $(\chi_1, \cdots, \chi_{k+1}) \in \mathbb{R}^{d \times (k+1)}$. To this aim, we start by noticing that the space $(\mathbb{P}_k(K))^d$ can be written as the direct sum $(\mathbb{P}_{k-1}(K))^d \oplus (\mathbb{P}_{[k]}(K))^d$, where $\mathbb{P}_{[k]}(K)$ is the space of polynomials having for degree exactly $k$. Therefore, the kernel of the operator $\varepsilon^k$ acting on functions belonging to $(\mathbb{P}_k(K))^d$ can be seen as the union of the kernels of the operator $\varepsilon^k$ acting either on functions of $(\mathbb{P}_{k-1}(K))^d$ or on functions of $(\mathbb{P}_{[k]}(K))^d$. Thus, for proving Proposition 3.2.1 it is enough to show that

$$(\mathbb{P}_{k-1}(K))^d \subset \text{Ker}(\varepsilon^k) \qquad (3.8a)$$

$$\begin{cases} (\mathbb{P}_{k-1}(K))^d \subset \text{Ker}(\varepsilon^k) \\ \forall u \in (\mathbb{P}_{[k]}(K))^d, \quad \sum_{i=1}^{d} x_i u_i = 0 \quad \Leftrightarrow \quad \varepsilon^k u = 0. \end{cases} \qquad (3.8b)$$
**Lemma 3.2.1** **First kernel's subspace determination**

$$ \forall u \in (\mathbb{P}_{k-1}(K))^d, \quad \varepsilon^k(u) = 0 $$

**Proof.**

Let us first point out that finding the functions $u \in (\mathbb{P}_{k-1}(K))^d$ verifying $\varepsilon^k u(\chi_1, \dots, \chi_{k+1}) = 0$ for any set of directions $(\chi_1, \dots, \chi_{k+1}) \in \mathbb{R}^{d \times (k+1)}$ reduces to find the set of $u \in (\mathbb{P}_{k-1}(K))^d$ so that $\varepsilon^k u(b^{p(1)}, \dots, b^{p(k+1)}) = 0$ for a given basis $\mathcal{B} = \{b^i\}_{i=1}^d$ of $\mathbb{R}^d$ and all mappings $p: \llbracket 1, k+1 \rrbracket \to \times_{i=1}^d \llbracket 1, d \rrbracket$ generating the $k+1$-multiset of $\llbracket 1, d \rrbracket$. By (3.7) and choosing the canonical basis for $\mathcal{B}$, expressing the kernel of $\varepsilon^k$ thus comes down to find every polynomial vector $u \in (\mathbb{P}_k(K))^d$ that verify

$$ \sum_{j=1}^{k+1} \frac{\partial^k u_{p_m(j)}}{\partial x_1^{\alpha_1^m} \cdots \partial x_d^{\alpha_d^m}} = 0 \quad \forall p_m \in \left\{ p : \llbracket 1, k+1 \rrbracket \to \times_{i=1}^d \left( \llbracket 1, j-1 \rrbracket \cup \llbracket j+1, d \rrbracket \right) \right\}, $$

where $\alpha^m = (\alpha_1^m, \dots, \alpha_d^m)$ is a multi-index such that $|\alpha^m| = k$ and $\alpha_i^m$ represents the number of occurrences of $i$ in $\{p_m(l), l \in \llbracket 1, j-1 \rrbracket \cup \llbracket j+1, d \rrbracket\}$. As this specification holds for any integers $k \in \mathbb{N}$, $d \ge 2$ and for any set $\{p_m\}_{m=1}^{(k+1)!}$, any polynomial $u \in (\mathbb{P}_{k-1}(K))^d$ belongs to the kernel of $\varepsilon^k$.

The inclusion (3.8a) follows immediately from Lemma 3.2.1. Let us now investigate the relation (3.8b) by reproducing the work of Nédélec [94].

**Theorem 3.2.1** **Second kernel's subspace determination**

$$ \{u \in (\mathbb{P}_{[k]}(K))^d, \sum_{i=1}^{d} x_i u_i \equiv 0\} \Leftrightarrow \{u \in (\mathbb{P}_{[k]}(K))^d \text{ s.t. } \varepsilon^k u \equiv 0\} $$

**Proof.**

▷ Let us first show that for any $u \in (\mathbb{P}_{[k]}(K))^d$ and $x \in \mathbb{R}^d$, $\varepsilon^k u \equiv 0 \Rightarrow x \cdot u = 0$. By assumption, we have for any $\chi = (\chi_1, \dots, \chi_{k+1}) \in \mathbb{R}^{d \times (k+1)}$ and for any $x \in \mathbb{R}^d$ that $\varepsilon^k u(x)(\chi) = 0$. Thus, taking in particular $\chi = \times_{i=1}^{k+1} x$, it comes $\varepsilon^k u(x, \dots, x) = 0$, which by definition of the differential operator reads

$$ \frac{1}{k+1} \sum_{j=1}^{k+1} x \cdot d^k u(\underbrace{x, \dots, x}_{k \text{ times}}) = x \cdot d^k u(\underbrace{x, \dots, x}_{k \text{ times}}) = 0. $$

In addition, since $u$ is homogeneous we obtain by the Euler's identity that $d^k u(\underbrace{x, \dots, x}_{k \text{ times}}) = k! u$. Plugging it in the expression of $\varepsilon^k$, we then get

$$ \varepsilon^k u(\underbrace{x, \dots, x}_{k+1 \text{ times}}) = 0 \quad \Rightarrow \quad d^k(\underbrace{x, \dots, x}_{k \text{ times}}) \cdot x = 0 $$

for any $x \in \mathbb{R}^d$, finally yielding $k! x \cdot u = 0$. Thus, for any $u \in (\mathbb{P}_{[k]}(K))^d$ and any $x \in \mathbb{R}^d$, $\varepsilon^k u \equiv 0 \Rightarrow \sum_{i=1}^d x_i u_i = 0$.
► Conversely, let us show that for any $u \in (\mathbb{P}_{[k]}(K))^d$ s.t. $x \cdot u \equiv 0, \varepsilon^k u = 0$. By assumption, we have that for any $x \in \mathbb{R}^d$, $x \cdot u = 0$. Hence, deriving $\varepsilon^k u$ in the particular direction $\chi = (x, \dots, x)$ yields

$$ \varepsilon^k u(\underbrace{x, \dots, x}_{k+1 \text{ times}}) = \frac{1}{k+1} \sum_{j=1}^{k+1} x \cdot d^k u(\underbrace{x, \dots, x}_{k \text{ times}}) = x \cdot d^k u(x, \dots, x), $$

which as $u$ is a polynomial of order exactly $k$ reduces to $\varepsilon^k u(\underbrace{x, \dots, x}_{k+1 \text{ times}}) = 0$. Expressing then $x$ on the canonical basis $B = \{b_i\}_{i=1}^d$ of $\mathbb{R}^d$, we can write $x \cdot u$ as $x \cdot u = \sum_{i=1}^d x_i b^i u_i$. Being is an homogeneous polynomial for which the term $\varepsilon^k u$ is identically vanishing, we can extend the above equality to any direction $\chi$, concluding the proof. ■

### Modern definition of the simplicial Raviart – Thomas space
Using the kernel parametrisation given in Proposition 3.2.1, one can now derive the classical expression of the Raviart–Thomas space.

**Proposition 3.2.1** *Simplicial Raviart–Thomas space*

The kernel characterisation Proposition 3.2.1 can also be expressed as

$$ RT_k = (\mathbb{P}_k(K))^d \oplus x \mathbb{P}_{[k]}(K), $$

defining the Raviart–Thomas space in its classical formulation.

**Proof.**

As described in Proposition 3.2.1, the Raviart–Thomas space of discretisation reads

$$ RT(K) = (\mathbb{P}_{k-1}(K))^d \oplus \{u \in (\mathbb{P}_{[k]}(K))^d, x \cdot u \equiv 0\}. $$

Thus, as the latter subspace only contains polynomials of degree $k$ whose gradient components are homogeneous polynomials of degree $k-1$, we obtain by Euler's identity that

$$ u \in (\mathbb{P}_{[k]}(K))^d \quad \text{and} \quad x \cdot u \equiv 0 \quad \Rightarrow \quad u = (k+1) x \nabla u. $$

Thus, as $\nabla p \in \mathbb{P}_{[k-1]}(K)$, one can identify the latter subspace with $x \mathbb{P}_{[k-1]}$ and derive immediately the classical expression of the Raviart–Thomas space of discretisation by shifting the order index by one. ■

The Raviart–Thomas space so defined, any element $p \in RT_k(K)$ writes under the form

$$ p = \begin{pmatrix} p_1 + x_1 q \\ p_2 + x_2 q \\ \vdots \\ p_d + x_d q \end{pmatrix} := \sum_{i=1}^{d} (p_i + x_i q) = \sum_{i=1}^{d} p_i + \sum_{i=1}^{d} x_i q \quad (3.9) $$

for some $p_i \in \mathbb{P}_k(K)$, $i \in [1, d]$ and $q \in \mathbb{P}_{[k]}(K)$.
**Raviart – Thomas space properties** For the space $RT_k(K)$ to be suitable for designing $H(\text{div}, K)$–conforming elements, it should allow the definition of a discretisation subspace on each face $f$ of the simplex’s boundary $\partial K$ where degrees of freedom enforcing conformity can be specifically designed. We thus present a dimension split and conforming properties guaranteeing it.

**Proposition 3.2.2** *Simplicial Raviart – Thomas space’s dimension*

$$ \dim RT_k(K) = \dim(\mathbb{P}_{k-1}(K)^d) + (d+1) \dim \mathbb{P}_k(f) $$

*Proof.* Combining the known results about the dimension of polynomial spaces with the classical expression of $RT_k(K)$ yield directly

$$ \dim RT_k(K) = \dim(\mathbb{P}_k(K)^d) + \dim \mathbb{P}_{[k]}(K) = \frac{1}{k!d!} (d(k+d)! + d(k+d-1)!) $$

Reorganising then the terms, it comes

$$ \begin{aligned} \dim RT_k(K) &= \frac{d}{dk!(d-1)!}k(k+d-1) + \frac{1}{k!(d-1)!}(1+d)(k+d-1)! \\ &= d \binom{k-1+d}{k-1} + (d+1) \binom{k+d-1}{k}, \end{aligned} $$

concluding the proof by definition of the $\mathbb{P}_k(K)$ and $\mathbb{P}_k(f)$ spaces.

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

In light of this dimension split, the boundary subspace is simply set as

$$ \mathcal{R}_k(\partial K) = \{p \in L^2(K), p|_{f_i} \in \mathbb{P}_k(f_i) \text{ for every face } f_i \in \partial K\}. $$

Let us now derive the divergence properties yielding $RT_k(K) \subset H(\text{div}, K)$.

**Theorem 3.2.2** *Divergence properties*

For any dimension $d \in \mathbb{N}$, $d \ge 2$ and any simplex $K \subset \mathbb{R}^d$, it holds

$$ \forall q \in RT_k(K), \quad \begin{cases} \text{div } q \in \mathbb{P}_k(K) \\ q \cdot n|_{\partial K} \in \mathcal{R}_k(\partial K) \end{cases} $$

Moreover, the divergence is onto from $RT_k(K)$ onto $\mathbb{P}_k(K)$.

*Proof.*

* Any $q \in RT_k(K)$ can be written $q = q_0 + xp_k$ with some $q_0 \in (\mathbb{P}_k(K))^d$ and $p_k \in \mathbb{P}_{[k]}(K)$. Thus, by Euler’s identity we have that $\text{div } q = \text{div } q_0 + \text{div}(x p_k)$ is a polynomial of degree $k$, and $\text{div } q \in \mathbb{P}_k(K)$.

* Let $f \in \partial K$ be a face of $K$ and $n$ its normal. By definition of the normal vector, we know that there exists a constant $c \in \mathbb{R}$ so that $\forall x \in f$, $x \cdot n \equiv c$, and thus $q \cdot n = q_0 \cdot n + p_k(x \cdot n)$. Recalling that $q_0 \in \mathbb{P}_k(K)$, $p_k \in \mathbb{P}_{[k]}(K)$ and that the definition of $\mathcal{R}_k(\partial K)$ does not require continuity across the faces, we obtain that $q \cdot n|_{\partial K} \in \mathcal{R}_k(\partial K)$.

* A proof of the surjectivity of the divergence operator from $RT_k(K)$ onto $\mathbb{P}_k(K)$ can be found in [103] and will not be repeated.

<figure><img src="image_36.png" /></figure>
### 3.2.1.3 Definition of simplicial Raviart – Thomas elements

The simplicial Raviart – Thomas space being a subset of $H(\mathrm{div}, K)$, a conforming element can be completely defined simply by designing a unisolvent set $\{\sigma_i\}_i$ of degrees of freedom, linear forms on $K$ that map $RT_k(K)$ to $\mathbb{R}$ and uniquely determines any discretised function (see Section 2.2.2).

**Degrees of freedom** To be unisolvent, the set of degrees of freedom must have the same dimension as $RT_k(K)$ and be linearly independent. Moreover, in our scope its definition should ensure the $H(\mathrm{div}, K)$ – conformity of the discretisation setting. Hence, one can take advantage of the dimension split obtained in Proposition 3.2.2 and notice that any face $f \in \partial K$ is contained in $\mathbb{R}^{d-1}$. Thus, $(1+d) \dim \mathbb{P}_k(f) = \dim \mathcal{R}_k(\partial K)$, and we can split the degrees of freedom into a normal and internal classification characterising functions living in $\mathbb{P}_{k-1}(K)^d$ and $\mathcal{R}_k(\partial K)$, respectively. There, the discretised functions are represented disjointly by moment-based degrees of freedom enforcing $H(\mathrm{div})$–conformity on the simplex boundaries and preserving it in the interior.

**Theorem 3.2.3** *Degrees of Freedom in the simplicial case*

For any $q \in RT_k(K)$, we set as degrees of freedom:

$$q \mapsto \int_f q \cdot n p_k \, d\gamma(x), \quad \forall f \in \partial K, \forall p_k \in \mathbb{P}_k(f) \qquad (3.10a)$$

$$q \mapsto \int_K q \cdot p_{k-1} \, dx, \quad \forall p_{k-1} \in (\mathbb{P}_{k-1}(K))^d \qquad (3.10b)$$

where $\gamma$ represents the paths skimming the faces is unisolvent for $RT_k(K)$ and $(K, RT_k(K), \{\sigma\})$ completely defines an element on $K$.

The unisolvence is immediate as the degrees of freedom are linear forms, $\dim(\{\sigma\}) = \dim RT_k(K)$ and $\mathrm{Ker}(\{\sigma\}) = \{0\}$. The kernel determination can be found e.g. in [103], not be repeated here for the sake of concision. ■

With such a clear separation of the degrees of freedom and the nature of the normal set, the $H(\mathrm{div})$–conformity of the elements immediately follows.

**Theorem 3.2.4** *Proposition*

The element $(K, RT_k(K), \{\sigma\})$ is $H(\mathrm{div})$–conforming.

By definition of the degrees of freedom and $\mathcal{R}_k(\partial K)$, the discretised function is entirely determined on a given face $f \in \partial K$ by the subset of (3.10a) containing moments supported on $f$ alone. Thus, any discretised quantity's normal component is determined identically on $f$ from any element $K$ sharing the face $f$. Since $q \cdot n \in \mathbb{P}_k(f)$, the $H(\mathrm{div})$–conformity follows. ■

Endowed with the set of degrees of freedom (3.10), the $RT_k(K)$ space is the smallest polynomial space within which the divergence maps $RT_k(K)$ to $\mathbb{P}_k(K)$.
**A word on the associated basis functions** Even though the set of basis functions associated to the simplicial Raviart – Thomas elements have not been explicitly shown here for the sake of concision [102], let us point out that they are the dual forms of $\{\sigma\}$ over the element $(K, RT_k(K), \{\sigma\})$. In particular, the normal component of basis functions associated to the internal degrees of freedom will vanish on the boundary, not containing any discretisation information there. Not vanishing inside the simplex, the normal component of the normal basis function will generate the space $\mathbb{P}_k(f)$ on each face.

In practice, a set of basis functions can be built by considering any basis of $\mathbb{H}_k(K)$ and tune it via the transfer matrix (2.58) built from their values through the set of degrees of freedom defining the Raviart–Thomas elements. It is indeed enough to transform the representation of $\mathbb{H}_k(K)$ by to fit the degrees of freedom representation (3.10a) and (3.10b) given a choice of basis $\{p_k\}$ generating $\mathcal{R}_k(\partial K)$ and $\{p_{k-1}\}$ generating $(\mathbb{P}_{k-1}(K))^d$, respectively. More convenient, a direct construction is equally possible when using a reference element through the use of Whitney forms (see e.g. [102] for more details).

**Graphical representation** Let us conclude the presentation of the simplicial Raviart – Thomas element by depicting its classical representation in Figure 3.1. Summarising the construction, the degrees of freedom' type and numbers are shown for a given discretisation order, the normal moments being represented by an arrow and the number of internal ones simply encircled.

<figure><img src="image_37.png" /><figcaption>**Figure 3.1.** Simplicial $RT_2(K)$ in two dimensions.</figcaption></figure>

### 3.2.1.4 Strengths and drawbacks

The simplicial Raviart–Thomas framework make use of a finite dimensional vectorial polynomial space whose dimension can be split into expressions that match the dimension of polynomial subspaces in $d$ and $d-1$ dimensions. The elements are thus built through the definition of two sets of degrees of freedom, characterising independently the boundary and inner parts of any discretised object, the $H(\mathrm{div})$-conformity being enforced by definition of the normal degrees of freedom. Defining a robust setting, one can either directly work with the set of degrees of freedom when defining numerical schemes, or design associated basis functions.

However, the design of the discretisation space is extremely rigid, dedicated to simplicial shapes only and preventing any direct adaptation to other geometries. Furthermore, when focusing on the discretisation quality, it can be seen that by construction the lowest order does not have internal moment, which can be troublesome in some applications. Retrieving the basis functions from a reference element may also be delicate, as a conforming mapping such as the Piola transform is required to preserve the element properties.
### 3.2.2 Quadrilateral Raviart – Thomas elements

The key property of the previously presented simplicial Raviart – Thomas elements is the adequacy between their dimension and the number of faces, allowing a dimension split yielding a possible construction of normal degrees of freedom on each face. If the space can be set up on any shape, such an element can only be built on simplices. A change of the discretisation space is therefore needed, done directly at the level of the scalar polynomial spaces.

**Raviart – Thomas spaces for quadrilateral geometries**

Considering the $d$-dimensional case and quadrilateral shapes, having $n = 2d$ faces. The Raviart – Thomas space is set up so that a dimension split similar to Proposition 3.2.2 can be established. A simple change allowing it is to consider the polynomial space $\mathbb{Q}_k$ instead of $\mathbb{P}_k$ when defining the scalar polynomial spaces, yielding to define $RT_k(K) \subset H(\text{div}, K)$ as

**Definition 3.2.3** *Quadrilateral Raviart – Thomas space*

Assuming $x \in \mathbb{R}^d$ and defining $\{\zeta_i\}_i$ is the set of the $d$ cyclic permutations of $\{k+1, \underbrace{k, \dots, k}_{d-1 \text{ times}}\}$,

$$RT_k(K) = \{u \in (\mathbb{Q}_k(K))^d, \varepsilon^k u = 0\} \\ = (\mathbb{Q}_k(K))^d + x \mathbb{Q}_{[k]}(K), \qquad (3.11a)$$

$$= \bigtimes_{i=1}^{d} \mathbb{P}_{\zeta_i(k+1, \dots, k, k)} \qquad (3.11b)$$

Indeed, it can be shown [94] that the kernel of $\varepsilon^k$ is parametrised by (3.11a), yielding the classical Raviart – Thomas space definition. However, as this definition is no longer in direct sum, one can reformulate it to (3.11b) for practical uses. So defined, the dimension of the space reads as follows.

**Proposition 3.2.2** *Dimension of $RT_k(K)$ on quadrilateral shapes*

$$\dim RT_k(K) = 2d \dim \mathbb{Q}_k(f) + \dim \bigtimes_{i=1}^{d} \mathbb{P}_{\zeta_i(\{k-1, \dots, k, k\})}$$

We can then define as in the simplicial case the boundary space as

$$\mathcal{T}_k(\partial K) = \{p \in L^2(K), p|_{f_i} \in \mathbb{Q}_k(f_i) \text{ for every face } f_i \in \partial K\}$$

and derive the following divergence properties

**Theorem 3.2.1** *Divergence properties*

For any $q$ belonging to $RT_k(K)$, it holds:

$$\begin{cases} \text{div } q \in \mathbb{Q}_k \\ q \cdot n|_{\partial K} \in \mathcal{T}(\partial K) \end{cases}$$

Having then $RT_k(K) \subset H(\text{div}, K)$, conforming elements can be designed.
**Degrees of freedom** To build the Raviart – Thomas elements on a quadrilateral of reference, we start by defining degrees of freedom ensuring the $H(\mathrm{div}, K)$ – conformity. In particular, by the above dimension split as $2d \dim \mathbb{Q}_k(f) = \dim \mathcal{T}_k(\partial K)$, we can write

$$ \dim RT_k(K) = \dim \bigtimes_{i=1}^{d} \mathbb{P}_{\zeta_i(\{k-1, \dots, k, k\})} + \dim \mathcal{T}_k(\partial K). $$

Therefore similarly as in the simplicial case, we can set up degrees of freedom respectively constructed from $\mathcal{T}_k(\partial K)$ and $\times_{i=1}^d \mathbb{P}_{\zeta_i(\{k-1, k, \dots, k\})}$, ensuring that the generated set is linearly independent. The shape of the degrees of freedom built from $\mathcal{T}_k(\partial K)$ are then selected so as to enforce $H(\mathrm{div})$ – conformity, yielding a similar result as in the simplicial case.

**Proposition 3.2.3** Degrees of freedom for quadrilateral elements

Representing by $\gamma$ the paths skimming the faces and considering any function $q \in RT_k(K)$, the set of degrees of freedom

$$ q \mapsto \int_{\partial K} q \cdot n p_k \, d\gamma(x), \quad \forall p_k \in \mathcal{T}_k(\partial K) \qquad (3.12a) $$

$$ q \mapsto \int_K q \cdot p_k \, dx, \quad \forall p_k \in \bigtimes_{i=1}^{d} \mathbb{P}_{\zeta_i(\{k-1, \underbrace{k, \dots, k}_{d-1 \text{ times}}\})} \qquad (3.12b) $$

define $(K, RT_k(K), \{\sigma\})$ as a finite element that is $H(\mathrm{div})$-conforming.

As in the simplicial case, the applications (3.12a) are called normal degrees of freedom, and (3.12b) internal degrees of freedom. Furthermore, any discretised function $q \in RT_k(K)$ is vectorial and can be decomposed on normal and internal basis functions, satisfying respectively $b \cdot n|_{\partial K} \in \mathcal{T}_k(\partial K)$ and $b \cdot n|_{\partial K} = 0$. The representation of the Raviart – Thomas element on quadrilaterals is given in Figure 3.2.

<figure><img src="image_38.png" /><figcaption>Figure 3.2. Representation of $RT_k(K)$ on a square, here for $k = 2$</figcaption></figure>

**Limitations of the Raviart – Thomas elements** As in the simplicial case, both the element geometry and discretisation design are extremely rigid. In particular, the quadrilateral Raviart – Thomas elements are specific to quadrilateral geometry without showing any possibility for geometrical extensions without having to change the definition of the discretisation space itself. Let us therefore present a more recent approach, which, if not directly $H(\mathrm{div})$–conforming, benefits from a great geometrical flexibility.
### 3.2.3 The classical Virtual Elements

Ab-initio developed in the context of two-dimensional elliptic problems, the Virtual Elements method has been introduced by Beirão da Veiga, Brezzi, Cangiani, Manzini, Russo and Marini in the recent work [7]. Further analysed and extended to higher dimensions and in [104, 105], the approach is in essence a generalization of the Finite Element method that borrows the point of view of modern Mimetic Finite Differences [106, 88]. Let us detail here in the two-dimensional case only for the sake of clarity.

**Virtual element spaces** Although using classically scalar discretisation spaces of finite dimension to approximate solutions to weak formulations of elliptic problems, the original virtual elements method does not rely on polynomial discretisation subspace of $Q_k(\mathbb{R}^d)$, but rather on subsets of solutions to Poisson problems whose second member are polynomials. So defined, the Virtual Element spaces of discretisation are naturally suitable to any geometry, the element shape being implicit in the definition of the Poisson problems.

More specifically, let us consider any polytope shape $K \subset \mathbb{R}^2$, possibly non-convex, on which to define a virtual element of order $k \in \mathbb{N} \setminus \{0\}$. A geometrically flexible finite-dimensional discretisation space is then defined as follows.

**Definition 3.2.4 Virtual element space**

For any $K \subset \mathbb{R}^2$, it is set

$$V_k(K) = \{v \in H^1(K), \Delta v \in \mathbb{P}_{k-2}(K), v|_{\partial K} \in \mathcal{V}_k(\partial K)\},$$

where

$$\mathcal{V}_k(\partial K) = \{v \in C^0(\partial K), \forall f \in \partial K, v \in \mathbb{P}_k(f)\}$$

Let us first observe that the considered boundary space $\mathcal{V}_k(\partial K)$ is non-empty, as containing all the continuous functions defined on $\partial K$ that are piecewise polynomial of degree $k$ on each face. In addition, as functions are uniquely determined by their values taken at the vertices and at $k-1$ distinct points within each face, it naturally comes $\dim \mathcal{V}_k(\partial K) = nk \neq 0$. Consequently, the Poisson problems are well-defined for any polytope having a reasonable aspect ratio, and the discretisation space $V_k(K)$ has a valid definition.

**Definition 3.2.5 Dimension of the virtual element space**

For any polytope $K \subset \mathbb{R}^2$ and any $k \in \mathbb{N} \setminus \{0\}$,

$$\dim V_k(K) = nk + \frac{k(k-1)}{2} = n + n(k-1) + \dim \mathbb{P}_{k-2}(K).$$

The space structure therefore allows, as in the Raviart – Thomas framework, to define degrees of freedom that are specifically designed for characterising either the inner or the boundary part of discretised quantities.
Moreover, the characterisation of the discretised functions on the boundary is local to each face, conferring a global continuity of any function discretised on the full computational domain $\Omega$. Polynomial discretisation spaces can also be embedded in the virtual element space, providing natural projectors and estimates [104]. In other words, it holds the following properties [7].

**Property 3.2.1** *Properties of the Virtual Elements discretisation*

For any $k \in \mathbb{N} \setminus \{0\}$, it holds

$$P_k \subset V_k.$$

Furthermore, for any two polytopes $K$ and $K^-$ belonging to $\mathbb{R}^2$ and such that $K \cap K^- = f \subset \mathbb{R}$, any two characterised functions $v_1 \in V_k(K)$ and $v_2 \in V_k(K^-)$ satisfy

$$v_{1|f} = v_{2|f}.$$

The discretisation space so defined, let us now introduce the virtual elements.

**Degrees of freedom** The virtual elements are classically defined through the set of degrees of freedom (3.13), providing here a mixed characterisation.

**Proposition 3.2.4** *Degrees of freedom defining the virtual elements*

Let $K \subset \mathbb{R}^2$ be a polytope and $k \in \mathbb{N} \setminus \{0\}$. Then, the set of degrees of freedom

$$q \mapsto q(v) \quad \forall v \text{ vertex of } K \qquad (3.13a)$$

$$q \mapsto q(x_i) \quad \forall \{x_i\}_{i=1}^{k-1} \text{ distinct points on the face } f, \forall f \in \partial K \qquad (3.13b)$$

$$q \mapsto \int_K q p_k \quad \forall p_k \in P_{k-2}(K) \qquad (3.13c)$$

is unisolvent in $V_k(K)$ and defines the Virtual element of order $k$.

More than separating the boundary characterisation (3.13a)-(3.13b) from the inner characterisation (3.13c), this choice also allows to work from a weak formulation of an elliptic problem without having to know the expression of the basis or even an explicit expressions for the discretised functions.

**Method strengths and drawbacks** Very promising, the recent Virtual Elements method is geometrically flexible and allows the particular use of general polytopal meshes. Its remarkable flexibility also extends to the variational level, as one can e.g. design arbitrary high order discretisations that are globally continuous. More involved settings can also be developed, as by example the serpentineity [107], divergence free [108], and $H(\operatorname{div})$-conforming ones. However, as any Finite Element-type method, this variational flexibility is somewhat restricted by the nature of the discretisation space $V_k(K)$, not allowing a lot of liberty in the type of degrees of freedom that one can use for defining an element by guaranteeing the set's unisolvence.
### 3.2.4 The $H(\operatorname{div})$-conforming virtual element space

More in line with our objective, the virtual elements are also available in a $H(\operatorname{div})$-conforming version, that has been originally developed for solving diffusion problems in mixed spaces [109, 110]. There, the discretisation space is specifically designed to enforce the regularity properties within the element, the definition of the degrees of freedom enforcing the cross-cell continuity of the normal component as in the Raviart – Thomas framework.

**Definition 3.2.6** $H(\operatorname{div})$-conforming virtual element space

For any $K \subset \mathbb{R}^2$, it is set

$$V_{2,k}^{face}(K) = \{v \in H(\operatorname{div}, K) \cap H(\operatorname{curl}, K) \quad \text{s.t.} \quad \forall f \in \partial K, v \cdot n|_f \in \mathbb{P}_k(f), \\ \nabla (\nabla \cdot v) \in \nabla \mathbb{P}_{k-1}(K), \quad \text{and} \quad \nabla \times v \in \mathbb{P}_{k-1}(K)\}.$$

Note that in this definition, and contrarily to the classical virtual element space, the discretised quantities are not necessarily continuous at the polytope's vertices. Furthermore, the space is generated by the set of solutions to problems

$$\left\{ \begin{array}{ll} \nabla \cdot v = p_1 & \text{on } K \\ \nabla \times v = p_2 & \text{on } K \\ v \cdot n = p_b, & \text{on } \partial K, \end{array} \right. \qquad (3.14)$$

defined for any $p_1, p_2 \in \mathbb{P}_k(K)$ and any $p_b$ face-wise polynomial of order $k$. The problem admitting a unique a solution for any polytope $K$ having a reasonable aspect ratio only if $\int_K p_1 \,dx = \int_{\partial K} p_b \,d\gamma(x)$, the space dimension reduces to the following expression.

**Property 3.2.2** *Dimension of the $H(\operatorname{div})$-conforming virtual space*

For any polytope $K \subset \mathbb{R}^2$ and any $k \in \mathbb{N} \setminus \{0\}$,

$$\dim V_{2,k}^{face}(K) = n \dim \mathbb{P}_k(f) + 2 \dim \mathbb{P}_{k-1}(K) - 1$$

The space $\mathbb{P}_{k-1}(K) \setminus \mathbb{R}$ being isomorphic to $\nabla \mathbb{P}_{k-2}$ through the gradient operator, a choice of degrees of freedom is thus naturally given by the Proposition 3.2.5.

**Proposition 3.2.5** *Degrees of freedom for conforming elements*

Let $K \subset \mathbb{R}^2$ be a polytope and $k \in \mathbb{N} \setminus \{0\}$. Then, defining the sets $G_k = \nabla \mathbb{P}_k$ and $G_k^{\perp}$ its orthogonal in $\mathbb{P}_k(K) \times \mathbb{P}_k(K)$, the set of degrees of freedom

$$q \mapsto \int_{f} q \cdot np_{k} \qquad \forall p_{k} \in \mathbb{P}_{k}(f) \qquad (3.15)$$

$$q \mapsto \int_{K} \operatorname{div}(q) p_{k-1} \qquad \forall p_{k-1} \in \mathcal{G}_{k-2}(K) \qquad (3.16)$$

$$q \mapsto \int_{K} \operatorname{div}(q) p_{k} \qquad \forall p_{k} \in \mathcal{G}_{k-2}^{\perp}(K) \qquad (3.17)$$

is unisolvent in $V_{2,k}^{face}(K)$ and defines a $H(\operatorname{div})$-conforming virtual element.
A graphical representation of the obtained conforming elements is pro-
vided in Figure 3.3. Lastly, let us mention that the well-definition of the space
and the unisolvence of the elements defined in Proposition 3.2.5, derived in
[109], allows in particular to define a natural projection operator onto the
space $\mathbb{P}_k(K) \times \mathbb{P}_k(K)$ and provide a rigorous finite element framework.

<figure><img src="image_39.png" /><figcaption>**Figure 3.3.** Graphical representation of a virtual conforming element for the orders $k = 1, k = 2$, and $k = 3$, respectively.</figcaption></figure>

**Method strengths and drawbacks** Comparably to the classical frame-
work detailed in Section 3.2.3, the $H(\operatorname{div})$-conforming version of the virtual
elements seems promising. However, even though the geometrical flexibil-
ity of the classical setting completely transfers to the conforming case, the
variational setting appears more restrictive and less customisable to further
specific properties or applications. Especially, the divergence properties of the
discretised functions at the boundary are enforced by controlling the regu-
larity and shape of their normal component directly in the definition of the
discretisation space. The uniqueness of the problem (3.14) thus prevents
the possibility of expressing a coordinate-wise behaviour, feature required for
grasping the subtleties of the discrete quantities’ behaviour on the boundary.
On a less critical note, one can also observe that the considered spaces are
built in link with $\mathbb{P}_k(K)$ polynomial spaces, providing less boundary freedom
than the classical quadrilateral Raviart – Thomas elements.

As we aim to provide the theoretical foundations guaranteeing the existence
and representativeness of the correction function (2.77) in a discrete setting
ensuring the conditions (2.78)–(2.82)–(2.135), we typically value this coor-
dinate wise freedom and developed a more variationally flexible framework by
combining the classical Virtual Element and Raviart – Thomas approaches.
## 3.3 A framework for arbitrary polytopes

We would like to extend the previously described Raviart – Thomas elements by drawing a general framework that allows the use of $H(\text{div})$-conforming elements any non-convex polytopes. In other words, for a given polytope $K \subset \mathbb{R}^d$ that is reasonably squeezed, we would like to define a finite dimensional discretisation space of order $k$, $\mathbb{H}_k(K)$, that satisfies the following property.

**Property 3.3.1 Discretisation space suited for $H(\text{div})$-conformity**

For any polytopal shape $K \subset \mathbb{R}^d$ in any dimension $d \in \mathbb{N}, \ge 2$, the space $\mathbb{H}_k(K)$ is a vectorial finite dimensional subset of $H(\text{div}, K)$ where any function $p \in \mathbb{H}_k(K)$ satisfies

$$ \left\{ \begin{array}{l} \text{div } p \in L^2(K) \\ p \cdot n|_{\partial K} \in \mathcal{H}_k(\partial K) \end{array} \right. $$

for some space

$$ \mathcal{H}_k(\partial K) = \{u|_{\partial K} \in L^2(\partial K), u|_f \in \mathcal{P}_{b,k}(f), \forall f \in \partial K\} $$

where $P_{b,k}(f)$ represents any polynomial subspace of $\mathbb{Q}(\mathbb{R}^{d-1})$ providing a discretisation of order $k$ on the face $f$.

Once one could design a discretisation space that fulfils the Property 3.3.1, it is enough for constructing a $H(\text{div})$-conforming element to define a set of degrees of freedoms whose expressions guarantee the continuity of the normal component of any discretised quantity across faces. In this setting, they can simply be defined face-wise accordingly to the expression of $P_{b,k}(f)$, the dimension of $\mathcal{H}_k(\partial K)$ being given by $\dim \mathcal{H}_k(\partial K) = n \dim(P_{b,k}(f))$ and dependent on the number of faces $n$.

In the following, we start by outlining the general idea in Section 3.3.1 deriving in Section 3.3.2 admissibility conditions for the discretisation space to be used in a $H(\text{div})$-conforming framework. We then develop a class of admissible discretisation spaces in Section 3.3.3, deriving simultaneously their associated boundary space $P_{b,k}(f)$. Possible sets of degrees of freedom yielding $H(\text{div})$-conforming elements are then proposed in Section 3.3.4, and examples of practical constructions are proposed in Section 3.3.6, on which numerical tests are performed in Section 3.3.8.
### 3.3.1 The general idea in a nutshell

As we aim to provide a general setting, we develop a flexible class of discretisation spaces whose dimension is adaptive to the number of faces rather than considering spaces whose nature change with the geometry. Furthermore, we describe the quality of the developed discretisations by a meaningful *order*, whose definition is linked to a preselected refinement pattern in the definitions of the space and degrees of freedoms. Hence, as the choice of this refinement pattern reduces to defining a series of discretisation spaces, we use the flexibility of the proposed elements to allow the definition of different frameworks fulfilling specific properties, as *e.g.* whether the lowest order space contains internal degrees of freedom or not.

To practically obtain those properties, one can however not directly extend the Raviart–Thomas elements to polytopes, the dimension of $RT_k(K)$ only depending on the desired order $k$ and not being adaptive to the number of the element faces. The Virtual elements framework is not directly suitable either, as if it allows the use of any general polytopal shapes, the variational definition of the elements appear to be rigid and restrictive, especially in their $H(\operatorname{div})$-conforming version.

We therefore develop geometrically and variationally flexible elements by combining the two approaches. More precisely, we construct the discretisation spaces by applying the vectorial structure of Raviart–Thomas to the geometrically flexible Virtual elements’ scalar space. The elements are then built through the definition of internal and normal degrees of freedoms, the latter enforcing the $H(\operatorname{div})$-conformity by imposing the normal component to be cross-face continuous. Lastly, series of spaces are defined by linking together the indices generating the second members of the Poisson problems defining the scalar Virtual element spaces.

### 3.3.2 Necessary conditions over spaces and elements

To begin with, let us derive step by step the admissibility conditions to build discretisation spaces that can endow $H(\operatorname{div})$-conforming elements. Recalling first that we aim to construct an unifying discretisation framework, we require the discretisation space $\mathbb{H}_k(K)$ to fulfil Admissibility conditions 3.3.1.

#### Admissibility conditions 3.3.1 *Requirements on the discretisation space*

The space $\mathbb{H}_k(K)$ is a finite dimensional vectorial subspace of $H(\operatorname{div}, K)$ whose dimension adapts to both the number of the polygon's faces and the discretisation order.

In addition, to be able to endow $H(\operatorname{div}, K)$-conforming elements $E_k(K)$ by defining a set of degrees of freedom, we have to keep in mind that the confor-
mity is only enforced by the definition of the normal degrees of freedom. We
therefore further ask Admissibility conditions 3.3.2.

Admissibility conditions 3.3.2 *Requirements on the elements*

1. For any space $\mathbb{H}_k(K)$, there exists a unisolvent set of degrees of freedom $\{\sigma\}$ that can be split into internal and normal subsets so that both the number of internal degrees of freedom and of the normal degrees of freedom per face do not depend on the shape of $K$.

2. The number of internal and normal degrees of freedom both increase strictly monotonically with the discretisation order.

Lastly, to ensure the existence of a split into internal and normal subsets of
degrees of freedom that match the normal and internal classification at the
level of the dual basis functions, the feasibility of Admissibility conditions
3.3.3 in $\mathbb{H}_k(K)$ is also required.

Admissibility conditions 3.3.3 *Requirement on the basis functions*

For any polytope $K$, the normal component of internal basis functions van-
ish on every face of the element, and the normal component of normal basis
function belong to $\mathcal{H}_k(\partial K)$.

Optionally, as it may be useful for some applications to work with an exclu-
sively normal characterisation, we aim to provide the possibility of designing
a discretisation framework within which the lowest order element is analo-
gous to $RT_0(K)$ for any number of faces, not having any internal moment.

Admissibility conditions 3.3.4 *Optional requirement on the basis functions*

The lowest order element has no internal degrees of freedom.

Given those requirements, let us first design a class of admissible discreti-
sation spaces before deriving possible definitions of degrees of freedom.

**3.3.3 A class of admissible approximation spaces**

Let us now progressively derive the discretisation spaces fulfilling the condi-
tions 3.3.1–3.3.4, before deriving their properties and extracting the subset of
them where $H(\mathrm{div})$-conforming elements can be defined.

**3.3.3.1 Construction of the $\mathbb{H}_k(K)$ discretisation spaces**

In order to design a subspace of $H(\mathrm{div}, K)$ that satisfies Admissibility condi-
tions 3.3.1, we are led to define a space $\mathbb{H}_k(K)$ with the same architecture as
the classical Raviart–Thomas space. We therefore look for spaces in the form
$\mathbb{H}_k(K) = (A_k)^d + x B_k$ for two given functional sets $A_k$ and $B_k$.
There, we start by observing that the use of polynomial spaces is excluded by the Admissibility conditions 3.3.3, being required for any number of faces. Therefore, we consider the spaces $A_k$ and $B_k$ based on solutions to Poisson's problems as in the context of the Virtual elements method [7]. There, a way to allow the existence in $\mathbb{H}_k(K)$ of smooth internal basis functions is to use the set of solutions to the boundary problems $\{u|_{\partial K} = 0, \Delta u = p_k\}$ for any $p_k$ belonging to $Q_m(K)$, $m \in \mathbb{N} \cup \{-1\}$. In addition, as the $H(\mathrm{div}, K)$-conformity will be enforced by normal quantities that are tested only on the boundaries, we also consider for given values of $l \in \mathbb{N} \cup \{-1\}$ the sets of Poisson's problems

$$ \left\{ \begin{array}{ll} \Delta u = 0 & \forall f \in \partial K, \forall p_k \in Q_l(f), \\ u|_{\partial K} = p_k \mathbb{1}_f & \end{array} \right. \qquad (3.18) $$

Thus, seeing the boundary $\partial K$ face-wise, we define the set

$$ \mathcal{H}_k(\partial K) = \{u|_{\partial K} \in L^2(\partial K), u|_f \in Q_k(f), \forall f \in \partial K\} \qquad (3.19) $$

and build the space $\mathbb{H}_k(K)$ as follows, for integers $l_1, l_2, m_1$ and $m_2$ whose values depend on $k$ and that will be determined in Section 3.3.3.3.

**Definition 3.3.1** *Definition of $\mathbb{H}_k(K)$ spaces*

$$ \mathbb{H}_k(K) = \{u \in H^1(K), u|_{\partial K} \in \mathcal{H}_{l_1}(\partial K), \Delta u \in Q_{m_1}(K)\}^d + x \{u \in H^1(K), u|_{\partial K} \in \mathcal{H}_{l_2}(\partial K), \Delta u \in Q_{[m_2]}(K)\} $$

### 3.3.3.2 Properties of the $\mathbb{H}_k(K)$ discretisation spaces

Having constructed the shape of the discretisation space so that the conditions 3.3.1 are fulfilled, let us now justify that $H(\mathrm{div})$-conforming elements can indeed be built by investigating its properties.

**Existence of any discretisation space** First of all, let us point out that as the boundary functions are in $L^2(\partial K)$, we get existence and uniqueness of the solution to any Poisson problems (3.18) for any polygon $K$ that is reasonably squeezed (see e.g. Theorem 3.5 in [111]). The scalar spaces encountered in Definition 3.3.1 are therefore never empty for any $(l_1, l_2, m_1, m_2) \in (\mathbb{N} \cup \{-1\})^4$, and any space $\mathbb{H}_k(K)$ is properly defined.

**Admissibility** Any space $\mathbb{H}_k(K)$ is constructed from four independent blocks whose definitions are driven by the independent coefficients $l_1, l_2$ and $m_1, m_2$. There, the couple $(m_1, m_2)$ drives the discretisation quality exclusively within the cell, while $(l_1, l_2)$ takes care only of the boundary. The internal and normal subspace classification is thus natural, and the Admissibility conditions 3.3.3 is fulfilled by construction. Furthermore, the Property 3.3.1 holds, emphasising that the $H(\mathrm{div}, K)$-conformity is ensured by the definition of $\mathcal{H}_k(\partial K)$, while the inner smoothness is provided through the Laplacian.
**Property 3.3.1** Divergence properties

For any function $q$ belonging to any space $\mathbb{H}_k(K)$, it holds

$$ \left\{ \begin{array}{l} q \cdot n_{|\partial K} \in \mathcal{H}_{\max\{l_1, l_2\}}(\partial K) \\ \nabla \cdot q \in L^2(K). \end{array} \right. $$

**Proof.**

• We start by deriving the first statement. By construction, any $q \in \mathbb{H}_k(K)$ can be decomposed into $q = q_0 + x q_1$ for some $q_0 \in (A_k)^d$ and $q_1 \in B_k$. Therefore, on the boundary of $K$ one has $q \cdot n_{|\partial K} = q_0 \cdot n_{|\partial K} + (x q_1 \cdot n)_{|\partial K}$. As the functions $q_1$ is scalar, this quantity can also read

$$ q \cdot n_{|\partial K} = q_0 \cdot n_{|\partial K} + q_1 (x \cdot n)_{|\partial K} $$

by linearity and commutativity of the dot product. Thus, Since for every face $f$ of $K$ the term $x \cdot n_{|f}$ is constant, it reduces to $q \cdot n_{|f} = q_0 \cdot n_{|f} + c_f q_{1|f}$ on each face $f$ for a constant $c_f \in \mathbb{R}$ depending only on the face layout and position with respect to the axes and origin. Therefore, since $q_{0|f} \in (\mathbb{Q}_{l_1}(f))^d$ and $q_{1|f} \in \mathbb{Q}_{l_2}(f)$, it comes $q \cdot n_{|f} \in \mathbb{Q}_{\max\{l_1, l_2\}}(f)$. And since it is valid for any face $f \in \partial K$, we finally get that $q \cdot n_{|\partial K} \in \mathcal{H}_{\max\{l_1, l_2\}}(\partial K)$

• Let us now derive the divergence property within the cell. Any $u \in \mathbb{H}_k(K)$ can be written under the form $u = \tilde{q} + x q$ for some functions $q \in H^1(K)$ and $\tilde{q} = (\tilde{q}_1, \dots, \tilde{q}_d)^T \in (H^1(K))^d$ such that

$$ \left\{ \begin{array}{l} \Delta q \in \mathbb{Q}_{[m_2]}(K) \\ q_{|\partial K} \in \mathbb{Q}_{l_2}(\partial K) \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} \Delta \tilde{q}_i \in \mathbb{Q}_{m_1}(K) \\ \tilde{q}_{i|\partial K} \in \mathbb{Q}_{l_1}(\partial K), \end{array} \right. \quad \forall i \in [1, d]. \quad (3.20) $$

We have

$$ \begin{aligned} \mathrm{div}(u) &= \sum_{i=1}^{d} \partial_{x_i} (x_i q) + \sum_{i=1}^{d} \partial_{x_i} \tilde{q}_i = \sum_{i=1}^{d} (q + x_i \partial_{x_i} q) + \sum_{i=1}^{d} \partial_{x_i} \tilde{q}_i \\ &= \underbrace{dq}_{\in L^2(K)} + \sum_{i=1}^{d} \underbrace{(x_i \partial_{x_i} q)}_{\in L^2(K)} + \sum_{i=1}^{d} \underbrace{\partial_{x_i} \tilde{q}_i}_{\in L^2(K)}. \end{aligned} $$

Since by (3.20) we have $\nabla \cdot q \in L^2(K)$, it comes that for any $i \in [1, d]$; $x_i \partial_{x_i} q \in L_{\mathrm{loc}}^2(K)$. As $K$ is compact and bounded, we have $L_{\mathrm{loc}}^2(K) = L^2(K)$ and $\mathrm{div} \, q \in L^2(K)$. As a by-product, note that we can derive $\nabla \cdot (x \nabla q) = \nabla \cdot q + x \Delta q$, where $\Delta q \in \mathbb{Q}_{\max[m_1, m_2+1]}$ and $x \Delta q \in C^\infty(K)$. ■

As for the Raviart–Thomas space, the above properties yield the following inclusion, allowing to design $H(\mathrm{div})$-conforming elements. The practical design of conforming element is then eased by a dimension split.

**Corollary 3.3.1** Subspace of $H(\mathrm{div})$

For any couples $(l_1, l_2)$ and $(m_1, m_2)$, $\mathbb{H}_k(K) \subset H(\mathrm{div}, K)$.
**Dimension** As in the Raviart – Thomas setting, the proposed construction benefits from a dimension split that will allow the classification of degrees of freedom into the internal and normal subcategories.

**Property 3.3.2 Dimension of the space $\mathbb{H}_k(K)$**

Whenever $l_1 \le 0$, the space $\mathbb{H}_k(K)$ has the following dimension.

$$ \dim \mathbb{H}_k(K) = n(d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1}) + (d(m_1 + 1)^d + (m_2 + 1)^d - m_2^d) $$

**Proof.**

As for any $l_1 \le 0$, the two natural subspaces are in direct sum, recalling the block construction of $\mathbb{H}_k(K)$ allows the dimension of the space $\mathbb{H}_k(K)$ to be easily derived. We can simply add the dimension of the two main subspaces $(A_k)^d$ and $x B_k$ to retrieve the dimension of $\mathbb{H}_k(K)$.

▷ First, we compute the dimension of $A_k$. In the way presented in [7], one can get it by using the superposition theorem. Indeed, for any second member belonging to $\mathbb{Q}_{m_1}$ and any boundary function $p_k \mathbb{1}_f \in L^2(K)$, there exists a unique solution to the Poisson's problems defining $A_k$, see e.g. [111]. Thus, reading out the structure of the set $A_k$ implies the following relation.

$$ \dim A_k = \dim \mathcal{H}_{l_1}(\partial K) + \dim \mathbb{Q}_{m_1}(K) = n(l_1 + 1)^{d-1} + (m_1 + 1)^d $$

Therefore, as $(A_k)^d$ is a simple Cartesian product of $d$ copies of $A_k$, we have immediately $\dim A_k = d(\dim A_k) = d(n(l_1+1)^{d-1}+(m_1+1)^d)$. We then retrieve in the exact same way the dimension of $B_k$, reading

$$ \dim B_k = \dim \mathcal{H}_{l_2}(\partial K) + \dim \mathbb{Q}_{[m_2]}(K) = n(l_2 + 1)^{d-1} + (m_2 + 1)^d - m_2^d. $$

▷ Lastly, we recall that the space $x B_k$ simply corresponds to an identical $d$-duplication of the space $B_k$ where each coordinate has been multiplied by the corresponding spatial variable. Thus, there is no liberty adjunction during its construction, and the dimension of $x B_k$ equals the one of $B_k$. By combining this, we finally get

$$ \begin{aligned} \dim \mathbb{H}_k(K) &= d \dim A_k + \dim B_k \\ &= d(n(l_1 + 1)^{d-1} + (m_1 + 1)^d) \\ &\quad + n(l_2 + 1)^{d-1} + (m_2 + 1)^d - m_2^d. \end{aligned} $$

One can observe in Property 3.3.2 that the term $(d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1}$ relates to the dimension of polynomial spaces in $d-1$ dimensions while $(d(m_1 + 1)^d + (m_2+1)^d - m_2^d)$ is linked to the dimension of polynomial spaces in $d$ dimensions. The structure of the space is therefore *a-priori* suitable for designing unisolvent sets of internal and normal degrees of freedom. Let us now make a step further and restrict range of $l_1, l_2, m_1$ and $m_2$ so that the generated discretisation spaces can endow $H(\mathrm{div}, K)$-conforming elements.
### 3.3.3.3 Admissible discretisation spaces

In order to define $H(\mathrm{div}, K)$-conforming elements from the definition of degrees of freedom in the spirit of Raviart-Thomas, we need to set $(d(l_1+1)^{d-1} + (l_2+1)^{d-1})$ normal degrees of freedom per face and $d(m_1+1)^d + (m_2+1)^d - m_2^d$ internal degrees of freedom. While this requirement has no impact on the admissibility of the coefficients $(m_1, m_2)$, it does reduce the range of admissible values that one can use for the coefficients $(l_1, l_2)$.

**Structure of admissible spaces** Any space $\mathbb{H}_k(K)$ is constructed from four independent blocks providing two distinct discretisations; on the boundary and within the element. Thus, when testing a function of $\mathbb{H}_k(K)$ through normal degrees of freedom, one can only retrieve the polynomial obtained from the two boundary conditions defining the sets $A_k$ and $B_k$. Especially, on each face, this polynomial is of the form $p = p_{k, A} + p_{k, B}$, where the function $p_{k, A} \in (A_k)_{|f}^d$ reads

$$p_{k, A} = \times_{j=1}^{d} \left( \sum_{|\alpha_i| \le l_1} a_{ij} x^{\alpha_i} \right) \qquad (3.21)$$

for a given set of multi-index $\{\alpha_i\}_i$ and coefficients $\{a_{i,j}\}_{i,j}$ that depend on the coordinates $x_j$. The function $p_{k, B} \in x B_k|_f$ reads however

$$p_{k, B} = \times_{j=1}^{d} \left( x_j \sum_{|\beta_i| \le l_2} b_i x^{\beta_i} \right) \qquad (3.22)$$

for a given set of multi-indices $\{\beta_i\}_i$ and coefficients $\{b_i\}_i$ independent of the coordinates $x_j$. Therefore, denoting by $\{\xi_j\}_{j \in \llbracket 1, d \rrbracket}$ the coordinates permutation that allows to shift the lowest orders terms of $x B_k|_f$ to $(A_k)^d|_f$, $p \in \mathbb{H}_k(K)|_f$ can be written as follows.

$$ \begin{aligned} \text{If } l_2 \ge l_1, \qquad & p = \times_{j=1}^{d} \left( \sum_{|\alpha_i| \le l_1} a_{ij} x^{\alpha_i} \right) + \times_{j=1}^{d} \left( x_j \sum_{|\beta_i| \le l_2} b_i x^{\beta_i} \right) \\ & = \times_{j=1}^{d} \left( \sum_{\substack{|\alpha_i| \le l_1 \\ |\alpha_i| \ne 0}} (a_{ij} + b_{\xi_j(i)}) x^{\alpha_i} + x_j \sum_{l_1 \le |\beta_i| \le l_2} b_i x^{\beta_i} \right) + \times_{j=1}^{d} a_{0j} x_j^0. \end{aligned} \qquad (3.23) $$

$$ \text{If } l_1 \ge l_2 + 1, \qquad p = \times_{j=1}^{d} \left( \sum_{\substack{|\alpha_i| \le l_1 \\ |\alpha_i| \ne 0}} (a_{ij} + b_{\xi_j(i)}) x^{\alpha_i} + a_{0j} x_j^0 \right). \qquad (3.24) $$

The structure of those relations implies that the terms $a_{ij}$ and $b_{\xi_j(i)}$ are combined into a single coefficient and cannot be specified individually from further normal degrees of freedom. Indeed, the remaining freedom can only be seen inside the polytope, as a consequence of the boundary conditions on the Poisson's solutions in either $A_k$ or $x B_k$. To prevent any over-determination
by the normal degrees of freedom in $\mathbb{H}_k(K)|_f$, we therefore have to make sure that the dimension of the boundary part (3.23)-(3.24) of any function living in $\mathbb{H}_k(K)$ is larger than the number of wished normal degrees of freedom per face. By reading out the structure of (3.23)-(3.24) it comes

$$ \dim \mathbb{H}_k(K)|_f = \begin{cases} (l_2 + 1)^{d-1} & \text{if } l_1 = -1, \\ d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1} - l_1^{d-1} & \text{if } l_2 \ge l_1, \\ d(l_1 + 1)^{d-1} & \text{otherwise.} \end{cases} \quad (3.25) $$

We thus restrict the admissible couples $(l_1, l_2)$ to those verifying the Admissibility conditions 3.3.5, preventing any over-determination.

### Admissibility conditions 3.3.5

#### Spaces allowing the use of conforming elements

Denoting $\dim \mathcal{N}$ the wished number of normal moments per face, and considering $\dim \mathbb{H}_k(K)|_f$ the number of coefficients that can be tuned from the face $f$, there should be

$$ \dim \mathcal{N} \le \dim \mathbb{H}_k(K)|_f. $$

In the case $l_2 \ge l_1$, it reduces to

$$ d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1} \le d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1} - l_1^{d-1} \quad (\Leftrightarrow \quad l_1^{d-1} \le 0), $$

while otherwise it comes

$$ d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1} \le d(l_1 + 1)^{d-1} \quad (\Leftrightarrow \quad l_2 = -1). $$

The Admissibility conditions 3.3.5 themselves are enough to ensure that the given $\mathbb{H}_k(K)$ space can endow a $H(\mathrm{div})$-conforming element. Indeed, by definition of $\mathbb{H}_k(K)$, the two subspaces $A_k$ and $x B_k$ do not intersect within the polytope, making any pair of coefficients $(m_1, m_2)$ suitable for providing the inner characterisation.

### Definition of series of spaces

While fulfilling the above conditions, one can set a specific discretisation framework within which the spaces share a predefined structure. By example, defining the four coefficients $l_1, l_2, m_1$ and $m_2$ through affine relations of the type $l = ak + b$ for some index $k \in \mathbb{N}$, the range of discretisation qualities achievable within the framework is predetermined by a refinement sequence in each block, and the order of each space can be simply defined as the index $k$ generating each of the four coefficients. A typical example is obtained by letting $m_1 = m_2 = k - 1, l_1 = 0$ and $l_2 = k$, leading to a series of discretisation spaces of order $k$ that always provide both internal and normal characterisations. Another natural definition would be to restrict the boundary characterisation to $l_1 = -1, l_2 = k$, yielding any normal component to be discretised as in the Raviart–Thomas setting. Fundamental, their complete descriptions are given in Sections 3.3.6.1 and 3.3.6.2, respectively.
### 3.3.4 Definition of $H(\mathrm{div})$-conforming elements

Under the Admissibility conditions 3.3.5, any space $\mathbb{H}_k(K)$ allows the construction of $H(\mathrm{div})$-conforming elements through the definition of normal and internal degrees of freedom, respectively enforcing and preserving conformity. We propose here a possible construction of such sets.

#### 3.3.4.1 Selection of Normal degrees of freedom

The role of the normal degrees of freedom is to define the characterisation of vectorial polynomials on the boundaries so that the $H(\mathrm{div})$-conformity of the element is enforced. We propose here a collection of suitable degrees of freedom any given coefficients $l_1$ and $l_2$, and derive admissibility conditions on extracted subsets that guarantee their unisolvence in $\mathbb{H}_k(K)$.

**Suitable definitions** Focusing on moment-based degrees of freedom whose support are restricted to a single face, we consider in particular definitions that are based on the projection of the tested quantity's normal component on polynomials living in $\mathbb{H}_k(K)|_f$. Namely, we consider the following expressions, applying to any discretised function $q \in \mathbb{H}_k(K)$.

1. The face integral of coordinate-wise components that are tested against polynomials

$$q \mapsto \int_f q_i n_i p \, d\gamma(x), \quad \forall p \in \mathbb{Q}_{\max\{l_1, l_2\}}(f), \qquad (3.26a)$$

2. The face integral of a function in $\mathbb{H}_k(K)$ projected onto the face normal, and tested against polynomials

$$q \mapsto \int_f q \cdot n_p \, d\gamma(x), \quad \forall p \in \mathbb{Q}_{\max\{l_1, l_2\}}(f), \qquad (3.26b)$$

$$q \mapsto \int_f q \cdot n_p \, d\gamma(x), \quad \forall p \in \{x_i p_i, p_i \in \mathbb{Q}_{\zeta_i([l_2], l_2, \dots, l_2)}(f)\}_{i \in [1, d-1]} \qquad (3.26c)$$

3. The pointwise values of the discretised quantity tested against the face's normal

$$q \mapsto q(x_{im}) \cdot n_i, \quad \text{for sampling points } \{x_{im}\}_m \text{ on the face } f_i. \qquad (3.26d)$$

Defining the normal degrees of freedom then reduces to choosing $d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1}$ of them among the possibilities (3.26) so that their set is unisolvent for $\mathbb{H}_k(K)|_f$. All the extracted degrees of freedom applied to any discretised function $q \in \mathbb{H}_k(K)$ will then provide the conformal representation associated to the generated element definition.
**Admissible sets of normal degrees of freedom** Extracting as many degrees of freedom as the dimension of $\mathbb{H}_k(K)|_f$ on each face among the above presented degrees of freedom yields automatically the unisolvence of the selected set in $\mathbb{H}_k(K)|_{\partial K}$ whenever no underdetermination is encountered. Selecting an unisolvent set of degrees of freedom among the above expressions thus reduces to avoiding the selection of linearly dependent projectors while paying attention to determining both global and coordinate-wise behaviours of any vector polynomial $q \in \mathbb{H}_k(K)|_f$.

In practice, it thus reduces to choosing the polynomials $p$ on which the discretised functions $q \in \mathbb{H}_k(K)$ will be tested coordinate-wise. The other polynomials $p$ are test functions for the global normal component $q \cdot n$. The unisolvence of the set is therefore ensured by the following admissibility conditions.

**Admissibility conditions 3.3.6** *Selection of degrees of freedom*

1. The projection polynomials $p$, and all the polynomials $\sigma: q \mapsto \sigma(q(x_{im}))$ that define the point values must be linearly independent.

2. When using a coordinate-wise degree of freedom of the type (3.26a), polygonal shapes $K$ containing a face parallel to any axis are not allowed. One of the terms $\{n_i\}_{i \in [1, d]}$ would indeed vanish for some $i \in [1, n]$, thus not describing any function of $\mathbb{H}_k(K)|_f$.

Note that here, the second point of the Admissibility conditions 3.3.6 may seem unreasonable as it may prevent the use of some geometrical shapes in specific orientations. However, it is always possible to modify the incriminated moments either easily by a local change in the definition of the projection vector, or more globally by redefining the moment so that it makes the element robust to rotations while staying $H(\mathrm{div})$-conforming. In the former case, as the incriminated moment is always of the coordinate-wise type, one can simply selecting another projection vector $v \neq n$.

Indeed, given that the dimension of the full set (3.26b) is always strictly smaller than the one of $\mathbb{H}_k(K)|_f$, there will necessarily exist in the selected set of degrees of freedom a global moment that involves a projection against the normal vector $n$, selected either from (3.26a), (3.26b) or even (3.26d). Each moment being supported on a single face, changing the projection direction in coordinate-wise moment has therefore no impact on the $H(\mathrm{div})$-conformity. We refer the interested reader to [63] for a detailed explanation.
**Impact of the normal degrees of freedom's selection** Due to the dual relationship between the degrees of freedom and the basis function, the selection of a particular set of degrees of freedom fulfilling the Admissibility conditions 3.3.6 rules the shape of the dual basis functions. One would then expect as many normal basis function as normal degrees of freedom, and respectively for the internal case. However, it is not necessary the case in the proposed setting, since a reclassification of some normal basis functions and associated degrees of freedom as internal ones is possibly encountered. Indeed, as the face-wise normal component of any function $q$ in $H_k(K)|_f$ is only of degree $\max\{l_1, l_2\}$, the term $q \cdot n|_f$ requires only $(\max\{l_1, l_2\}+1)^{d-1}$ basis functions to be decomposed on. Therefore, up to $d (l_1+1)^{d-1}$ basis functions may see their global normal component vanishing on every face. Their coordinate-wise components will however not vanish, as they take care of the coordinate-wise behaviours that cannot be determined solely through the expression of $q \cdot n|_f$. Typically, the more global degrees of freedom are designed, the more the representation of $p \cdot n$ is completed globally. As a consequence, more basis functions have a vanishing normal component as they are forced to take care only of coordinate-wise behaviours, forcing them to be reclassified into internal basis function. A reverse scenario may equally be considered.

**Example**

Let us consider the two-dimensional case where $l_1 = l_2 = 0$. There, any discretised function $p \in H_k(K)|_f$ reads

$$q = \begin{pmatrix} A \\ B \end{pmatrix} + C \begin{pmatrix} x \\ y \end{pmatrix}$$

for some constants $A, B$ and $C$ belonging to $\mathbb{R}$. The characterisation of $q$ on a face $f$ can be then done by selecting two component-wise moments involving $A n_1$ or $B n_2$ tested against the constant polynomial $p = 1$, and one global moment that tests $q \cdot n = C (n_1 + n_2) + A n_1 + B n_2$ against the polynomial $p = x$. The representation of the function's normal component being completed component-wise, this choice would yield one basis function that vanishes on the boundary, and two basis functions having a constant normal component. Reversely, one could also choose to select two global and one coordinate-wise moments, yielding the normal component of two basis functions to vanish on the boundary, the third one being constant. Note however that all basis functions whose normal component vanishes on the face $f$ have a first order polynomial behaviour component-wise.

Hence, one can select the degrees of freedom given the properties wished for the associated basis functions. Furthermore, should one like to avoid any reclassification, allowing to draw a parallel with the Raviart-Thomas framework from the lowest order space on, one can meet the Admissibility conditions 3.3.4 by choosing wisely the value of $l_2$, as discussed in Section 3.3.6.2.
**A first example of the normal degrees of freedom' selection** In order to support the above statements and in view of later use, let us detail a possible selection of normal degrees of freedom for any face $f \in \partial K$ in the case $l_2 \ge l_1$, where every function in $\mathbb{H}_k(K)_{|f}$ takes the form (3.23).

We here choose to extract a subset of moment-based degrees of freedom from (3.26) so that the functions discretised in $\mathbb{H}_k(K)_{|f}$ are determined as much as possible from their normal component, the remaining freedom being characterised by few coordinate-wise moments. We thus select for each face

$$ \sigma: q \mapsto \int_f q_i n_i x_i^{l_1+1} d\gamma, \quad \forall i \in [1, d], \qquad (3.27a) $$

$$ \sigma: q \mapsto \int_f q \cdot n p_k d\gamma, \quad \forall p_k \in \mathbb{Q}_{l_2}(f) \setminus \mathbb{Q}_{l_1}(f), \qquad (3.27b) $$

$$ \sigma: q \mapsto \int_f q \cdot n x_j x_j^{l_2} \tilde{x} d\gamma, \quad \forall j \in [1, d], \forall \tilde{x} \in \mathbb{Q}_{l_2}(\partial_j K), \qquad (3.27c) $$

where $\tilde{x} \in \mathbb{Q}_{l_2}(\mathbb{R}^{d-2})$ is not involving the variable $x_j$ so that the moment (3.27c) has for integrands the second terms of the right hand side of (3.23) when $|\beta_i| = l_2$. Note that the set (3.27) is of dimension $d(l_1+1)^{d-1} + (l_2+1)^{d-1} - (l_1+1)^{d-1} + (d-1)(l_2+1)^{d-2}$, while we require $d(l_1+1)^{d-1} + (l_2+1)^{d-1}$ moments. This configuration can therefore only be used when $l_1$ and $l_2$ satisfy

$$ (l_1 + 1)^{d-1} \le (d-1)(l_2 + 1)^{d-2} \qquad (3.28) $$

criterion that is a simple reduction of the Admissibility conditions 3.3.5. In two dimensions, the above relation reduces to an equality, and all the degrees of freedom presented in (3.27) are considered. In higher dimensions, a further selection from the set (3.27) is required. There, we consider the sets (3.27a)-(3.27b) fully and select any $(l_1 + 1)^{d-1}$ moments from (3.27c).

**Definition 3.3.2 Moment based normal degrees of freedom set**

Any choice of $(l_1 + 1)^{d-1}$ moments among (3.27c) is denoted as the “configuration $Ia$”.

So defined, the configuration $Ia$ yield the non-degenerating normal basis function to match the ones of the Raviart – Thomas setting, though up to a different scaling. To also obtain a similar scaling, one should rather scale the above degrees of freedom with respect to each face's area and orientation, or consider in place of the moments (3.27c) the pointwise values

$$ q \mapsto q(x_{im}) \cdot n, \forall i \in [1, d], m \in [1, (1+l_1)^{d-1}], \qquad (3.29) $$

$x_{im}$ being any sampling point on the face $f_i$, allowing to control directly both the scaling and offset of the basis function's normal components.

**Definition 3.3.3 Mixed-type normal degrees of freedom set**

Any selection of $(l_1+1)^{d-1}$ degrees of freedom among the sets (3.27a),(3.27b) and (3.29) is labelled as the “configuration $Ib$”.
**A second example of the normal degrees of freedom' selection** In order to corroborate the observation made on the impact of the degrees of freedom's choice, let us provide a second example where the focus is made on a coordinate-wise characterisation of the normal component, the remaining freedom being characterised by few global moments. Still considering $l_2 \ge l_1$, a selection along those lines arises naturally from the equation (3.23), providing for each face $f \in \partial K$ the following degrees of freedom.

$$ \sigma: q \mapsto \int_f q_i n_i d\gamma, \quad \forall i \in [1, d] \qquad (3.30a) $$

$$ \sigma: q \mapsto \int_f q \cdot v x^{\alpha_l} d\gamma, \quad \forall |\alpha_l| = l_1 + 1, \text{ given } v \in \mathbb{R}^d, v \ne n \qquad (3.30b) $$

$$ \sigma: q \mapsto \int_f q \cdot n x_i x^{\alpha_l} d\gamma, \quad \forall i \in [1, d-1], \forall l_1 < |\alpha_l| \le l_2, \qquad (3.30c) $$

There, the moment (3.30b) has been slightly modified from (3.30c) for the Admissibility conditions 3.3.6.2 to be satisfied, without impacting the $H(\text{div})$-conformity. So defined, the total number of degrees of freedom available in the layout (3.30) is $\dim((3.30)) = d(l_1 + 1)^{d-1} + (d-1)(l_2 + 1)^{d-1} - (d-1)l_1^{d-1}$, larger than the required $d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1}$. We therefore need a further extraction, strengthening in the process the Admissibility conditions 3.3.5 to

$$ (d-1)l_1^{d-1} \le (d-2)(l_2+1)^{(d-1)}. \qquad (3.31) $$

As $l_1 \le 0$, this inequality always holds. In particular, in two dimensions the bound (3.31) reduces to an equality, and the set (3.30) has to be considered fully. In higher dimensions, not achieving the bound provides a choice in the projection space, as any extraction of $(l_2 + 1)^{d-1}$ moments is admissible. A natural choice is to follow the relation (3.31) and select the sets (3.30a) and (3.30b) entirely, while extracting any $(l_2 + 1)^{d-1}$ elements from (3.30c).

**Definition 3.3.4 Moment based normal degrees of freedom set**

Any extraction of $(l_2 + 1)^{d-1}$ moments among (3.30) is referred to as the "configuration $IIa$".

As in the previous example, one can consider a definition using point values instead of the moments (3.30a), using in place the expressions

$$ \sigma: q \mapsto q_i(x_m) n_i. \forall i \in [1, d], x_m \text{ middle or barycentric point in } f. \qquad (3.32) $$

More than directly scaling the dual basis functions and their offset, the above expressions provide a direct control of their shape. One can thus e.g. weight or offset the above expression to enforce positivity of the basis functions.

**Definition 3.3.5 Mixed-type normal degrees of freedom set**

Any extraction of $(l_2+1)^{d-1}$ degrees of freedom among the sets (3.32), (3.30b) and (3.30c) is labelled as the "configuration $IIb$".
### 3.3.4.2 Selection of internal degrees of freedom

Let us now derive a characterisation of any function living in $\mathbb{H}_k(K)$ within the element through the definition of admissible internal degrees of freedom. Their design is here flexible, but delicate to determine. Indeed, we have to make sure that the corresponding internal basis functions fulfil the Admissibility conditions 3.3.3 for any number of faces $n$. We therefore stick to the idea of Raviart–Thomas and define moment based degrees of freedom that read for any $q \in \mathbb{H}_k(K)$

$$ \sigma(q) \mapsto \int_K q \cdot p_k \, dx, \quad \text{for all } p_k \in \mathcal{P}(K) \qquad (3.33) $$

for some function space $\mathcal{P}(K)$ of dimension $((m_1 + 1)^d + (m_2 + 1)^d)$. So defined, the internal degrees of freedom are moments supported strictly within the element. Thus, given the definition of $\mathbb{H}_k(K)$ strictly within the polytope $K$, being in particular independent of the number of face, the dual basis functions will necessarily vanish at the boundaries.

In practice, considered as a test space, the space $\mathcal{P}(K)$ may simply gather polynomial functions used in the definition of the Poisson's problems generating $\mathbb{H}_k(K)$. The discretised quantities would then be determined through their polynomial projections. Another possible choice is to test the functions living in $\mathbb{H}_k(K)$ against a set of Poisson's solutions to the problems

$$ \left\{ \begin{array}{l} \Delta p_k \in \mathcal{P}(K) \\ p_{k|_{\partial K}} = 0 \end{array} \right. $$

Using one or the other possibility for defining the projector space, the unisolvency of the internal degrees of freedom set in $\mathbb{H}_k(K)_{|\tilde{K}}$ is ensured by the following admissibility conditions (see *Proof of Proposition 3.3.1*, third part).

**Admissibility conditions 3.3.7 Selection of internal degrees of freedom**

1. The polynomials $\{p_k\}_k$ generating $\mathcal{P}(K)$ are linearly independent.

2. No polynomial belonging to the set $\{p_k\}_k$ generating $\mathcal{P}(K)$ has a degree larger than $\max\{m_1, m_2 + 1\}$.

Having stated admissibility conditions for designing admissible sets of normal and internal degrees of freedom, let us define $H(\mathrm{div})$-conforming elements over the discretisation space $\mathbb{H}_k(K)$. To conclude the construction, unisolvency proofs will be given, both for the four presented examples and for general sets of degrees of freedom satisfying the derived admissibility conditions.
### 3.3.4.3 $H(\mathrm{div})$-conforming elements

Let us now combine the two previous subsections with the definition of the space $\mathbb{H}_k(K)$ and define a class of $H(\mathrm{div}, K)$-conforming polytopal elements.

**Proposition 3.3.1** *$H(\mathrm{div})$-conforming elements*

* Let $K$ be any polytope and $\mathbb{H}_k(K)$ be a finite dimensional space given by Definition 3.3.1 that verifies the Admissibility conditions 3.3.5.

* Let also $\{\sigma_N\}$ be any selection of $d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1}$ degrees of freedom from the set (3.26) fulfilling first item of the Admissibility conditions 3.3.6, and $\{\sigma_I\}$ the set of internal degrees of freedom built through the expression (3.33) for any projection set $\mathcal{P}(K)$ fulfilling the admissibility conditions 3.3.7.

Then, the set $\{\sigma_N\} \cup \{\sigma_I\}$ is unisolvent for the discretisation space $\mathbb{H}_k(K)$ and defines a $H(\mathrm{div}, K)$-conforming element.

The well definition of the elements is an immediate corollary of the following proposition, while the $H(\mathrm{div})$-conformity is immediate by definition of the face-wise support of the normal degrees of freedom and the polynomial nature of the space $\mathbb{H}_k(K)$ on the boundary.

**Proposition 3.3.2** *Unisolvence of the derived elements*

Let $q \in \mathbb{H}_k(K)$, and denote $\sigma_N(q)$ the n-tuple of normal degrees of freedom extracted from the set (3.26) for each face $f \in \partial K$. Then,

$$ \left\{ \begin{array}{l} \sigma_N(q) = 0 \\ \displaystyle\int_K q \cdot p_k \,dx = 0 \quad \text{for all } p_k \in \mathcal{P}, \end{array} \right. $$

implies $q = 0$.

**Example 3.3.1** *Unisolvence of the four $H(\mathrm{div})$-conforming example elements*

Combined with any selection of internal degrees of freedom fulfilling the Admissibility conditions 3.3.7, the sets of degrees of freedom defining the configurations $I$ and $II$ in their variants $a$ and $b$, given respectively in Definitions 3.3.2, 3.3.3 and Definitions 3.3.4,3.3.5, are unisolvent for $\mathbb{H}_k(K)$. ♦

Let us now prove Proposition 3.3.2 by first showing the unisolvence of the four sets of degrees of freedom proposed in Example 3.3.1. The arguments are then generalised and the proof extended to any admissible configuration in Proof of Proposition 3.3.2.
### Proof of the element unisolvence stated in the Example 3.3.1

All along this proof, we refer to the functions $p_k$ by the term “kernel”, while using the term “integrand” to represent the term $q \cdot p_k$. Immediate transfer of this designation apply to the normal moment based degrees of freedom.

Let us first sketch the proof. To begin with, let us point out that the key lies in the Admissibility conditions 3.3.6 ensuring the linear independence of the set of pointwise values and moment's integrands. The linearity of the integral operators transfers then this independence to the moments themselves, characterising any function of $\mathbb{H}_k(K)$ independently on the boundary and within the cell. We proceed in three steps.

1. First, we show that the internal characterisation of the function does not impact the normal one, allowing the determination to be done distinctively within the element and on the boundary.

2. Then, we show that selecting the appropriate number of degrees of freedom in any of the sets $Ia$ or $Ib$ ensures a unique characterisation on the boundary. We use the fact that the kernels are scalar polynomials while the functions of $\mathbb{H}_k(K)$ are vector polynomials.

3. Lastly, we consider the interior of the element where the characterisation is done through projections over linearly independent sets. Those projections of functions in $\mathbb{H}_k(K)$ are indeed neither identically null nor identically identical (i.e. they differ at least on a subset of non-zero measure).

▷ **Step 1.** Let us first recall that the space $\mathbb{H}_k(K)$ is built from blocks of independent functions. In particular, the boundary behaviour of functions living in $\mathbb{H}_k(K)$ is independent of their behaviour within the inner cell. Therefore, by the structure of $\mathbb{H}_k(K)$ and making use of the superposition theorem, any function $q \in \mathbb{H}_k(K)$ reads

$$q = h \mathbb{1}_{\partial K} + g \mathbb{1}_{\mathcal{K}} \qquad (3.34)$$

for two functions $h$ and $g$ belonging to $\mathbb{H}_k(K)$. As a consequence, characterising a function $q \in \mathbb{H}_k(K)$ comes down to characterising the independent functions $f$ and $g$ on the distinct supports $\partial K$ and $\mathcal{K}$, respectively. Note also that necessarily, $h|_f \in \mathbb{X}_{i=1}^d \mathbb{Q}_{\max\{l_1, l_2+1\}}(\mathbb{R}^{d-1})$ for any face $f \in \partial K$. We show that any admissible extraction in the sense of the Admissibility conditions 3.3.6 from either of the two sets of degrees of freedom ($Ia$, internal), ($Ib$, internal) fully characterises the functions $h$ and $g$, independently. In all the following, the notation $Ia$ or $Ib$ refers to the corresponding set of normal degrees of freedom while internal refers to the set (3.33) and is identical to any of the two configurations under consideration.
We first show that any above defined set of degrees of freedom preserve the independence of the boundary and inner characterisations. To this aim, we combine the relation (3.34) with the all possible definitions of the degrees of freedom. It comes that all global normal moments lead to an expression of the form

$$ \sigma(q) = \int_f q \cdot n p_k = \int_f (h \mathbb{1}_f + g \mathbb{1}_{\hat{K}}) \cdot n p_k = \int_f h \cdot n p_k $$

for some polynomial function $p_k$ living on $\partial K$. On the other side, as $x_m \in f$, the global degrees of freedom that are built from pointwise values read

$$ \sigma(q) = q(x_m) \cdot n = h(x_m) \cdot n \mathbb{1}_f(x_m) + g(x_m) \cdot n \mathbb{1}_{\hat{K}}(x_m) = h(x_m) \cdot n. $$

Similar relations for coordinate-wise degrees of freedom can be derived, i.e.

$$ \sigma(q) = \int_f q_i n_i p_k = \int_f (h_i \mathbb{1}_f + g_i \mathbb{1}_{\hat{K}}) n_i p_k = \int_f h_i n_i p_k $$

$$ \text{and} \quad \sigma(q) = q_i(x_m)n_i = h_i(x_m)n_i \mathbb{1}_f(x_m) + g_i(x_m)n_i \mathbb{1}_{\hat{K}}(x_m) = h_i(x_m)n_i, $$

where here the terms $q_i$ simply represent the $i^{th}$ component of the function $q$. Therefore, in any of the configurations $Ia$ and $Ib$ no contribution of the function $g$ representing the inner part of the cell is involved in the normal degrees of freedom. The mirror case is obtained with the internal moments, leading via (3.34) to

$$ \sigma(q) = \int_K q \cdot p_k = \int_K (h \mathbb{1}_{\partial K} + g \mathbb{1}_{\hat{K}}) \cdot p_k = \int_K g \cdot p_k, $$

where $p_k$ stands for any Poisson's function living in $\mathbb{H}_k(K)$ or any polynomial function defining the second member of a Poisson's problem involved in the definition of $\mathbb{H}_k(K)$. There, the function $h$ representing the boundary part of the function $q$ is not involved, that for any definition of the space $P_k$ generating the internal moments. Thus, by linearity we can decompose the degrees of freedom $\{q \mapsto \sigma_i(q)\}_i$ in the following matrix.

$$ \text{Normal Dofs values} \quad \left\{ \begin{array}{c} \sigma_1 \\ \vdots \\ \sigma_{N_N} \\ \sigma_{N_N+1} \\ \vdots \\ \sigma_{N_I} \end{array} \right\} = \quad \left( \begin{array}{c} \text{Normal moments} \\ \text{applied to } h \\ 0 \\ \text{Internal moments} \\ \text{applied to } g \end{array} \right) \quad \left( \begin{array}{c} h \\ g \end{array} \right) \\ \text{Internal Dofs values} $$

Clearly, there is no interconnection between the function's characterisation on the boundaries and the one performed within the element. Thus, showing Proposition 3.3.2 reduces to prove independently that either $(Ia) = 0$, $(Ib) = 0$, $(IIa) = 0$ or $(IIb) = 0$ implies $h|_{\partial K} = 0$, and that $\int_K g \cdot p_k x = 0$ for all $p_k \in P_k$ implies $g|_{\hat{K}} = 0$.
▷ **Step 2.** Let us now consider the boundary characterisation. There, by definition of the spaces $\mathcal{H}_{l_1}$ and $\mathcal{H}_{l_2}$, the function $q|_{\partial K}$ is discontinuous at the polytope's vertices and can be decomposed into $n$ vectorial polynomial functions $\{h_j\}_j$ with distinct supports, each of them matching one particular face of the polytope. Thus, we can write

$$q|_{\partial K} = \sum_{j=1}^{n} r_j \mathbb{1}_f$$

with $r_j \in \times_{i=1}^d \mathbb{Q}_{\max\{l_1, l_2+1\}}(f)$ and its support $f$ any face belonging to $\partial K$. With a similar argument than in the previous point, the characterisation of $q|_{\partial K}$ can therefore be done face-wise, and the determination matrix becomes block-diagonal. We discuss here the characterisation on one particular face $f$ by showing the invertibility of the corresponding matrix block, the arguments naturally transposing to the other ones. In this perspective, we show that for any $r_j \in \times_{i=1}^d \mathbb{Q}_{\max\{l_1, l_2+1\}}(f)$, it holds $\{(Ia)|_f = 0 \text{ or } (Ib)|_f = 0 \text{ or } (IIa)|_f = 0 \text{ or } (IIB)|_f = 0\} \Rightarrow r_j = 0$, where $(\cdot)|_f$ represents the subset of the degrees of freedom $(\cdot)$ whose support (or evaluation point for point-values) matches (or lies on) the face $f \in \partial K$.

To begin with, let us recall that on the face $f$ the function $r_j$ is a multi-valued polynomial of the form

$$r_j|_f = \begin{pmatrix} a_{0,1} \\ \vdots \\ a_{0,d} \end{pmatrix} + \sum_{i=\dim(\mathcal{H}_0)}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} \begin{pmatrix} b_{\xi_1(i)} + a_{i,1} \\ \vdots \\ b_{\xi_d(i)} + a_{i,d} \end{pmatrix} m_{\alpha_i(x)} + \sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} \begin{pmatrix} x_1 b_i \\ \vdots \\ x_d b_i \end{pmatrix} m_{\alpha_i(x)},$$

where $m_{\alpha_i}$ represents a monomial of $\mathbb{Q}_{\max\{l_1, l_2\}}$ of multi-index degree $\alpha_i$ such that the set $\{m_{\alpha_i}\}_{i=\dim(\mathcal{H}_0)}^{\dim(\mathcal{H}_m)}$ forms a base of $\mathcal{H}_m \setminus \mathcal{H}_l$. Note that the coefficients $\{a_{ij}\}_{i,j}$ are defined coordinate-wise while the coefficients $\{b_i\}_i$ are identical for all the components. The function $r_j$ is therefore determined by

$$\dim(\{\{a_{i,m}\}_{i \in [0, \dim(\mathcal{H}_{l_1} \setminus \mathcal{H}_0)]}, \{b_i\}_{i \in [\dim(\mathcal{H}_{l_1} \setminus \mathcal{H}_0), \dim(\mathcal{H}_{l_2})]}\})$$

coefficients. Going further by exploiting the fact that in all configurations, the function $r_j$ is determined only through its normal components, we can specify the above equation. Denoting $n = (n_1, \dots, n_d)$ the normal vector to the $j$th face $f$, it comes

$$r_j \cdot n|_f = \sum_{m=1}^{d} a_{0,m} n_m + \sum_{i=\dim(\mathcal{H}_0)}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} \sum_{m=1}^{d} (b_{\xi_m(i)} + a_{i,m}) n_m m_{\alpha_i}(x) + \sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} c_j b_i m_{\alpha_i}(x),$$

where $c_j = x \cdot n$ is a constant term on the face $f$. Reordering the terms, we
end up with the formulation

$$
r_j \cdot n_{|f} = \sum_{m=1}^{d} \left( \left( a_{0,m} + \sum_{i=\dim(\mathcal{H}_0)}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} a_{i,m} m_{\alpha_i}(x) \right) n_m \right) \\
+ \left( \sum_{m=1}^{d} n_m \right)^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} \sum_{i=1}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} b_{\xi_m(i)} m_{\alpha_i}(x) + c_j \left( \sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} b_i m_{\alpha_i}(x) \right). \tag{3.35}
$$

The structure of the retrieved form makes clearly emerge the coefficients
that should be used depending on the coordinate-wise behaviour of the
polygon. In addition, as all the coefficients determining $r_j$ appear in this
expression, using degrees of freedom defined only from the normal compo-
nents of tested functions is admissible. Thus, the two configurations fitting
this framework, we only have to make sure that the set of extracted degrees
of freedom are uniquely characterising each of the involved coefficients. To
this aim, we explicit all the possible degrees of freedom when applied to
the function $r_j$. For the sake of clarity, we denote by $\{\sigma_{M_{i,l}}\}_{il}$ the moments
designed coordinate-wise, being of the form

$$
\sigma_{M_{i,l}} : q \mapsto q_i(x_{jl}) n_{jx_i} \quad \text{or} \quad \sigma_{M_{i,l}} : q \mapsto \int_f q_i n_{jx_i} p_l
$$

for some scalar polynomial $p_l$, and by $\{\sigma_{T_l}\}_l$ the ones acting globally, reading

$$
\sigma_{T_l} : q \mapsto q(x_l) \cdot n \quad \text{or} \quad \sigma_{T_l} : q \mapsto \int_f q \cdot n \, p_l
$$

for some scalar polynomial $p_l$. Further, for convenience we denote by $\{\sigma_{V_l}\}_l$
the global degrees of freedom that comes into play to determining the
coordinate-wise coefficients, whose expressions are done in the same way
as $\{\sigma_{T_l}\}_l$. We now express those degrees of freedom depending on the co-
efficients $\{b_{i,m}\}_{im}$ and $\{a_i\}_i$. Using the linearity of the degrees of freedom,
plugging the expression (3.35) in place of $q$ and setting the permutation op-
erator directly on the multi - indices $\alpha_i$ instead of the coefficients $a_i$, we can
rewrite the moments as follows.

$$
\sigma_{M_{m,l}} : (\{a_{i,m}\}, \{b_i\}) \mapsto b_{0,m} \int_f n_m p_l + \sum_{i=\mathcal{H}_0}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} a_{i,m} \int_f m_{\alpha_i} n_m p_l \\
+ \sum_{i=\mathcal{H}_0}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} b_i \int_f n_m m_{\xi_m(\alpha_i)} p_l \\
+ \sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} b_i \int_f (x_m n_m m_{\alpha_i}) p_l
$$
$$
\begin{align*}
\sigma_{T_l} : (\{a_{i,m}\}, \{b_i\}) &\mapsto \sum_{m=1}^{d} a_{0,m} \int_f n_m p_l + \sum_{m=1}^{d} \sum_{i=\dim(\mathcal{H}_0)} a_{i,m} \int_f m_{\alpha_i} n_m p_l \\
&\qquad + \sum_{i=\mathcal{H}_0}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} b_i \sum_{m=1}^{d} \left( \int_f n_m m_{\xi_m(\alpha_i)} p_l \right) \\
&\qquad + \sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} b_i \int_f c_{m\alpha_i} p_l
\end{align*}
$$

Thus, defining the component-wise parts of the global moments $\sigma_{T_l}$ by
$\sigma_{T_{m,l}}(q) = \int_f n_m q p_l$ such that $\sigma_{T_l} = \sum_{m=1}^d \sigma_{T_{m,l}}$, one can express any degrees
of freedom of the two considered sets as

$$
\sigma_{M_{m,l}} : q \mapsto a_{0,m} \sigma_{M_{m,l}}(1) + \sum_{i=\mathcal{H}_0}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} a_{i,m} \sigma_{M_{m,l}}(m_{\alpha_i}) + \sum_{i=\mathcal{H}_0}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} b_i \sigma_{M_{m,l}}(m_{\xi_m(\alpha_i)}) \\
\qquad + \sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} b_i \sigma_{M_{m,l}}(x_m m_{\alpha_i})
$$

$$
\text{and} \quad \sigma_{T_l} : q \mapsto \sum_{i=\mathcal{H}_0}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} b_i \sum_{m=1}^{d} \sigma_{T_{m,l}}(m_{\xi_m(\alpha_i)}) + \sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} a_i c_j \sigma_{T_l}(m_{\alpha_i}) \\
\qquad + \sum_{i=\mathcal{H}_0}^{\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})} \sum_{m=1}^{d} a_{i,m} \sigma_{T_{m,l}}(m_{\alpha_i}) + \sum_{m=1}^{d} a_{0,m} \sigma_{T_{m,l}}(1).
$$

Note that in view of deriving the determination matrix, the last term can
also be decomposed as follows.

$$
\sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} b_i c_j \sigma_{T_l}(m_{\alpha_i}) = \sum_{i=\dim(\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2})}^{\dim(\mathcal{H}_{l_2})} b_i \sum_{m=1}^{d} \sigma_{T_{l,m}}(x_m m_{\alpha_i}).
$$

Similar relations for $\sigma_V$ can be derived from the expression of $\sigma_T$. Thus, we
can rewrite the degrees of freedom as a dot product and derive the charac-
terisation matrix $\Sigma$

$$
(\sigma_{M_{1,1}}, \dots, \sigma_{M_{d,l}}, \sigma_{V_1}, \dots, \sigma_{V_l}, \sigma_{T_1}, \dots, \sigma_{T_l})^T = \Sigma (\{a_{i,m}\}_{im}, \{b_i\}_i)^T
$$

which shape is given in the figure 3.4. We now investigate its structure.
First of all, as the number of extracted degrees of freedom from the two sets
$Ia$ and $Ib$ matches the number of coefficients determining $r_j$, the matrix $\Sigma$
is a square matrix.
<figure><img src="image_40.png" /><figcaption>C: Group of coordinate-wise moments associated to a same quantification process (same testing monomial in moments or evaluation at a same same point-value).
(*) Tuning globally coefficients that have repercussion on the coordinate-wise coefficients $a_i$ (misc moments).
As many rows as moments $\sigma_{v_{i,m}}$ with $m \in [\dim(\mathcal{H}_1 \setminus \mathcal{H}_0) + 1, \dim(\mathcal{H}_1)]$</figcaption></figure>

<caption>**Figure 3.4.** Layout of the determination matrix associated to the proposed shapes of normal degrees of freedom. Each highlighted block corresponds to a given type of degrees of freedom, characterising either the global behaviour of their normal component or the coordinate-wise one.</caption>
Let us focus on the top two-by-two left blocks, surrounded in blue. They correspond to the coefficients that should be determined coordinate-wise. Thus, by construction, there are $\dim(\{a_{im}\}_{im}) = d(l_1 + 1)^{d-1}$ columns. And by the definition of the configurations *II* and *I*, we have $\dim(\{\sigma_{M_{i,j}}\}_{ij}) = d + d((l_1 + 1)^{d-1} - 1) = d(l_1 + 1)^{d-1}$. Therefore, this submatrix is a square matrix. Furthermore, each subblock corresponds to one member of the decomposition of $q$ tested through coordinate-wise degrees of freedom whose kernels are built on the same monomial. Therefore, as the degrees of freedom $\{\sigma_{M_{i,j}}\}$ consider one normal component only, the coefficients $\{a_{im}\}_{im}$ for $m \neq j$ are not involved, and the subblocks are diagonal. Thus, those submatrices are invertible and in particular their columns and rows are linearly independent.

On the other side of the matrix, the last bottom block surrounded in deep red matches the Raviart – Thomas moments tuning members of $\mathbb{H}_k(K)|_{\partial K}$ living exclusively in $xB_{k|f}$. It is then a submatrix of the classical Raviart – Thomas’ one, and is thus invertible. In particular, its rows and columns are linearly independent. The extended bottom right submatrix highlighted in dashed red corresponds to the previously described high-order submatrix of the Raviart – Thomas’s setting, enriched by the moments $\{\sigma_V\}$ tuning the behaviour of members of $\mathbb{H}_k(K)|_f$ falling in the intersection $\mathcal{H}_{l_1} \cap \mathcal{H}_{l_2}(f)$. This extended submatrix is thus equivalent to the full Raviart – Thomas setting, and therefore invertible. Indeed, even if the moments $\{\sigma_V\}$ have to be slightly modified from the Raviart – Thomas setting in the configuration *II*, this modification leaves the projection order unchanged, and the integrand still belongs to $\mathbb{H}_k(K)|_f$.

There is thus only left to show that there is no linear dependence between rows of different row blocks. And there, as the degrees of freedom are linear forms, it is enough to show that the integrand of moments (or polynomials constructing the pointwise values) that involve the same monomial are linearly independent. Indeed, being linear forms whose kernels are polynomials, the degrees of freedom can combine each other only if their integrand ($q$ tested against the kernel) involve – up to constants – the same monomials. We then have to show that in both configurations, the rows involving terms whose projection onto the kernel can be expressed from a same monomial are linearly independent.

In the configuration *I*, this property comes automatically. Indeed, the only interaction between degrees of freedom having integrands sharing the same monomial order (and then possibly being based on the same monomial) is possible between (3.27a) and (3.27b) when $|p_k| = l_1 + 1$. Indeed, by definition of $\mathbb{H}_k(K)$, the polynomial $p_k \cdot n$ in (3.27a) is only of order $l_1$.
However, no combination of (3.27b) can form the moments (3.27a). Indeed, for any real coefficients $c_i$ it holds

$$ \sum c_i q_i n_i x_i^{l_1+1} \neq q \cdot n p_k $$

for any monomial $p_k$ such that $|p_k| = l_1 + 1$. Note that in the left hand side, all the $c_i$ should be non-null to reconstruct the term $p \cdot n$. However, doing so no factorisation by a single monomial such that

$$ \sum_{i=1}^{d} c_i q_i n_i x_i^{l_1+1} = \left( \sum_{i=1}^{d} c_i q_i n_i \right) p_k $$

is possible. Thus, the designed moments are linearly independent, and no row combination can occur for any tested polynomial belonging to $\mathbb{H}_k(K)|_{\partial K}$.

In the configuration $II$, the integrands involving the same monomials are found in the definitions (3.32) and (3.30c) when $l_1 = 0$ (resp. (3.30a) and (3.30c)). As the problem arises on constants, linear combinations would be possible in a classical Raviart – Thomas setting as the terms $p_i n_i$ and $p_i n_i$ can be combined to form $p \cdot n$. However, the use of the vector $v$ not collinear to $n$ in place of $n$ for the lowest order moment in the equation (3.30c) makes the setting fulfil the Assumption 3.3.6. Therefore, those degrees of freedom are not linearly dependent. All in all, for both configurations all the rows are linearly independent. As by construction we have as many relations as unknowns, the matrix is invertible. Thus, we get a null kernel, meaningly,

$$ (Ia) = 0 \text{ or } (Ib) = 0 \text{ or } (IIa = 0) \text{ or } (IIb = 0) \text{ implies } h|_{\partial K} = 0. $$

The proposed boundary characterisations are therefore comprehensive.

### ▷ Step 3.
Let us now consider the internal characterisation of functions living in $\mathbb{H}_k(K)$, and from the first step, it is enough to study the characterisation of $g$ within the inner cell. By definition of $\mathbb{H}_k(K)$, any function $g \in \mathbb{H}_k(K)|_{\partial K}$ can be decomposed over a set of Poisson's solutions as follows.

$$ g = \sum_{i=1}^{d} \sum_{j=1}^{\dim A_k} a_{i,j} e_j u_i + \sum_{i=1}^{\dim B_k} b_j x \tilde{u}_i $$

Here, the vector $e_j$ stands for $e_j = (0, \dots, 1, 0, \dots, 0)^T$ where the 1 is in the $j$th position. The functions $u_i$ and $\tilde{u}_i$ represents the Poisson's solutions of the problems

$$ \left\{ \begin{array}{l} \Delta u_i = p_i, p_i \in \mathbb{Q}_{m_1}(K) \\ u_i|_{\partial K} = 0 \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} \Delta \tilde{u}_i = \tilde{p}_i, \tilde{p}_i \in \mathbb{Q}_{[m_2]}(K) \\ \tilde{u}_i|_{\partial K} = 0 \end{array} \right. \qquad (3.36) $$

where $\{p_i\}_i$ and $\{\tilde{p}_i\}_i$ form respectively a basis of $\mathbb{Q}_{m_1}(K)$ and $\mathbb{Q}_{[m_2]}(K)$. In
any presented definition of the degrees of freedom, the internal characteri-
sation is done through moment-based degrees of freedom of the form

$$
\sigma_{I_k} : q \mapsto \int_K q \cdot p_k \, dx
$$

where the kernels $p_k \in \mathcal{P}_k$ consist of linearly independent polynomials be-
longing to $\mathbb{Q}_{\max\{m_1, m_2+1\}}(K)$, or of the solution of their corresponding prob-
lems of the form (3.36). Therefore, we can derive a characterisation matrix
in the same spirit as in the case of the normal characterisation.

<figure><img src="image_41.png" /><figcaption>**Figure 3.5.** Layout of the determination matrix associated to the proposed shapes of internal degrees of freedom.</figcaption></figure>

Let us consider the case where $\mathcal{P}_k$ forms a polynomial projection space.
There, none of the $p_k \in \mathcal{P}_k$ is the zero function. In the same time, the func-
tions $\{\{u_i\}_i, \{x\tilde{u}_i\}_i\}$ are linearly independent, and being solutions to some
Poisson's problem with non-zero second member, they are by construction
not identically vanishing on $K$. Indeed, even when $m_2 < m_1$ where sec-
ond members of the problems (3.36) lives both in $\mathbb{Q}_{m_1}$ and $\mathbb{Q}_{[m_2]}$, it holds
$\Delta(x\tilde{u}_i) = 2\nabla \cdot \tilde{u}_i + \Delta(\tilde{u}_i)$. Thus, it is impossible to combine linearly the func-
tion $x\tilde{u}_i$ with functions of the set $\{u_i\}_i$. Furthermore, the degrees of the
polynomials belonging to the space $\mathcal{P}_k$ are lower or equal than the highest
degree of the second members of the Poisson's problem defining the space
$\mathbb{H}_k(K)$. Thus, every projection of function of $\mathbb{H}_k(K)$ onto the space $\mathcal{P}_k$ is not
null. And as the internal moments are linear forms, any linear combina-
tion of those moments at fixed $p_k$ could have its integrand factorised by the
kernel $p_k$ for any $p_k \in \mathcal{P}_k$, transferring the linear independency of the set
$\{\{u_i\}_i, \{x\tilde{u}_i\}_i\}$ to the terms $\{\int u_i \cdot p_k\}_i$ for any fixed $p_k \in \mathcal{P}_k$.

Lastly, as the space $\mathcal{P}_k$ contains only linearly independent functions
the previous argument can be repeated for each row of the matrix shown
in Figure 3.5. And as by construction the number of internal degrees of
freedom matches the dimension of the space $\mathbb{H}_k(K)|_{\check{K}}$, the linear indepen-
dence of functions of $\mathcal{P}_k$ combined with the linear independence of the tested
functions transfer automatically to the moments tested against a basis of
$\mathbb{H}_k(K)|_{\check{K}}$. Thus, the internal submatrix is invertible. The same reasoning
can be applied when $\mathcal{P}_k$ is built from the Poisson's solutions themselves, as
the projections of functions would decompose the functions directly.
▷ **Wrap-up.** Merging the above points together, we get that $Ia = 0$ or $Ib = 0$ or $IIa = 0$ or $IIb = 0$ implies $h|_{\partial K} = 0$ and $\int_K g \cdot p_k \,dx = 0$ for all $p_k \in \mathcal{P}_k$ implies $g|_{\tilde{K}} = 0$. Thus, from the first step we get that for $q \in \mathbb{H}_k(K)$ $Ia = 0$ or $Ib = 0$ or $IIa = 0$ or $IIb = 0$ and $\int_K q \cdot p_k \,dx = 0$ for all $p_k \in \mathcal{P}_k$ implies $q = 0$, concluding the proof.

Having obtained the unisolvency of the given for example, let us generalise the arguments to any configuration satisfying the Admissibility conditions 3.3.6 and Admissibility conditions 3.3.7.

**Proof of Proposition 3.3.2 and Proposition 3.3.1**

The proof of Proposition 3.3.2 is a straightforward generalisation of the one presented in the context of the four examples $Ia, Ib, IIa$ and $IIb$ in Proof of Proposition 3.3.1. The only change lies in the extraction of the degrees of freedom, which impacts the matrix only on the top left two by two blocks describing the coordinate-wise behaviours. As the extraction fulfils the Admissibility conditions 3.3.6, the rows involving terms whose projection onto the kernel can be expressed from a same monomial are linearly independent. Thus, the same arguments as above can be applied.

In particular, by this admissibility criterion there cannot be more than $d+1$ polynomials reducing to the same moments' kernel or to an equivalent point-value quantifier. Thus, as we have $d+1$ coordinate-wise moments to tune per decomposed monomial, there is no over-determination at a fixed polynomial degree. The constraint on the extraction of degrees of freedom ensures the non over-determination overall. Furthermore, the linear independence of the sub-matrix's columns is ensured as those polynomials cannot be linearly dependent. Thus, by linearity of the degrees of freedom, the independence of the kernels transfers to the moments and there is no row dependency. The submatrix block corresponding to any specific order is therefore invertible, and the same conclusion as in the Proof of Proposition 3.3.1 follows.

The proposition 3.3.1 therefore holds by the number of selected degrees of freedom, matching the dimension of the space $\mathbb{H}_k(K)$. Furthermore, one can observe that the admissible normal degrees of freedom are face-wise supported, and that the space $\mathbb{H}_k(K)|_f$ is of polynomial nature. The constructed elements are therefore $H(\mathrm{div})$-conforming, concluding the proof and the construction of the elements.

Considering the definition of $\mathbb{H}_k(K)$ as a discretisation space, any admissible definition of degrees of freedom thus yields a $H(\mathrm{div}, K)$-conforming element.
### 3.3.5 Summary of the construction

Before going to practical case studies and comparison with existing $H(\mathrm{div})$-conforming settings, let us summarise the construction of the $\mathbb{H}_k(K)$ space and recapitulate the example definitions of normal degrees of freedom.

**Discretisation space** The proposed class of discretisation spaces reads

$$ \mathbb{H}_k(K) = \{u \in H^1(K), u|_{\partial K} \in \mathcal{H}_{l_1}(\partial K), \Delta u \in \mathbb{Q}_{m_1}(K)\}^d \oplus x \{u \in H^1(K), u|_{\partial K} \in \mathcal{H}_{l_2}(\partial K), \Delta u \in \mathbb{Q}_{[m_2]}(K)\}, $$

with the convention that $\mathbb{Q}_{-1} = \{0\}$ and where the integers $l_1, l_2, m_1$ and $m_2$ verify
$$ m_1, m_2, l_2 \ge -1 \quad \text{and} \quad -1 \le l_1 \le 0. $$

So defined, it holds

$$ \dim \mathbb{H}_k(K) = n(d(l_1 + 1)^{d-1} + (l_2 + 1)^{d-1}) + ((m_1 + 1)^d + (m_2 + 1)^d - m_2^2) $$

$$ \forall q \in \mathbb{H}_k(K), q \cdot n|_{\partial K} \in \mathcal{H}_{\max\{l_1, l_2\}}(\partial K) \quad \text{and} \quad \mathbb{H}_k(K) \subset H(\mathrm{div}, K). $$

Thus, conforming elements can be defined through normal degrees of freedom enforcing the $H(\mathrm{div})$-conformity and internal ones preserving it, provided that the polytope $K$ satisfies the two conditions

* The polytope $K$ has a reasonable aspect ratio, so that the Poisson problems defining the underlying Virtual elements spaces are well-posed.

* No face is parallel to any axis, ensuring the unisolvence of the presented sets of degrees of freedom when using component-wise expressions.

The above condition on the orientation of the face with respect to the axis raises stability issues when dealing with element whose faces are almost parallel to the axis and using component-wise degrees of freedom. However, this issue is easily avoidable by selecting at least a global degree of freedom involving the term $p \cdot n$ in the moment's integrand and changing the testing vector in the coordinate-wise degrees of freedom to any vector $v \neq n$.

**Internal degrees of freedom** The internal characterisation of functions $q \in \mathbb{H}_k(K)$ relied entirely on moment-based degrees of freedom taking the form

$$ \sigma: q \mapsto \int_K q \cdot p_k \, dx, \quad \forall p_k \in \mathcal{P}, $$

where $\mathcal{P}$ is a $(m_1+1)^d + (m_2+1)^d - m_2^d$ dimensional space that can be defined either as a polynomial space or as a set of solutions to Poisson problems having a polynomial second member. Either case, its definition should be such that

* The polynomials $\{p_k\}_k$ generating $\mathcal{P}(K)$ are linearly independent.

* Any polynomial $p \in \{p_k\}_k$ satisfies $\deg(p) \le \max\{m_1, m_2+1\}$.
**Normal degrees of freedom** The boundary characterisation of functions $q \in \mathbb{H}_k(K)$ relies on global and coordinate-wise degrees of freedom that can be either moment-based or point values. If typical choices are directly reported in Table 3.1, the general selection rule reads as follows.

In order to ensure the unisolvence of the normal degrees of freedom in $\mathbb{H}_k(K)|_{\partial K}$ and the $H(\text{div})$-conformity of the elements, their expressions are considered face-wise. Exactly $d(l_1 + 1)^{d-1} + (l_2 + 1)^d$ of them are thus selected among the set (3.26) for each face $f \in \partial K$.

<table><caption>**Table 3.1.** Summary of the used degrees of freedom for the four considered configurations. There, the integrals sum up along the boundary path $\gamma$.</caption><thead><tr><th>Configuration</th><th>Ia</th><th>Ib</th><th>IIa</th><th>IIb</th></tr></thead><tbody><tr><td>Available when</td><td colspan="2">(d - 2)(l<sub>2</sub> + 1)<sup>d-1</sup> ≥ l<sub>1</sub><sup>d-1</sup>(d - 1)</td><td colspan="2">(d - 1)(l<sub>2</sub> + 1)<sup>d-2</sup> ≥ (l<sub>1</sub> + 1)<sup>d-1</sup></td></tr><tr><td>Free choice</td><td colspan="2">Select (l<sub>2</sub>+1)<sup>d-1</sup> moments per face from the bold ones</td><td colspan="2">Select (l<sub>1</sub>+1)<sup>d-1</sup> moments per face from the bold ones</td></tr><tr><td>Low order representation</td><td>∫<sub>f</sub> q<sub>i</sub> n<sub>i</sub><br>∀i ∈ ⌊1, d⌋<br>∀f ∈ ∂K</td><td>q<sub>j</sub>(x<sub>im</sub>) n<sub>jx<sub>i</sub></sub><br>∀i ∈ ⌊1, d⌋<br>∀i ∈ ⌊1, n⌋<br>x<sub>m</sub> midpoint</td><td>Inherited from the highest order representation</td><td>Inherited from the highest order representation</td></tr><tr><td>Representation of elements in A<sub>k</sub> ∩ B<sub>k|∂K</sub></td><td>∫<sub>f</sub> q<sub>i</sub> n<sub>i</sub> x<sup>β<sub>j</sub></sup><br>∫<sub>f</sub> q ⋅ v x<sup>α<sub>j</sub></sup><br>∀i ∈ ⌊1, d⌋<br>∀f ∈ ∂K<br>∀1 ≤ |β<sub>j</sub>| ≤ l<sub>1</sub>,<br>|α<sub>j</sub>| = l<sub>1</sub> + 1,<br>v ≠ n<sub>i</sub></td><td>∫<sub>f</sub> q<sub>i</sub> n<sub>i</sub> x<sup>β<sub>j</sub></sup><br>∫<sub>f</sub> q ⋅ v x<sup>α<sub>j</sub></sup><br>∀i ∈ ⌊1, d⌋<br>∀f ∈ ∂K<br>∀1 ≤ |β<sub>j</sub>| ≤ l<sub>1</sub>,<br>|α<sub>j</sub>| = l<sub>1</sub> + 1,<br>v ≠ n<sub>i</sub></td><td>∫<sub>∂K</sub> q<sub>i</sub> n<sub>i</sub> p<sub>k</sub><br>∫<sub>f</sub> q<sub>i</sub> x<sub>i</sub><sup>l<sub>1</sub>+1</sup> n<sub>i</sub><br>∀i ∈ ⌊1, d⌋<br>∀f ∈ ∂K,<br>p<sub>k</sub> ∈ H<sub>l<sub>1</sub></sub>(∂K) \ H<sub>0</sub>(∂K)</td><td>∫<sub>∂K</sub> q<sub>i</sub> p<sub>k</sub> n<sub>i</sub><br>∫<sub>f</sub> q<sub>i</sub> x<sub>i</sub><sup>l<sub>1</sub>+1</sup> n<sub>i</sub><br>∀i ∈ ⌊1, d⌋<br>∀f ∈ ∂K,<br>p<sub>k</sub> ∈ H<sub>l<sub>1</sub></sub>(∂K) \ H<sub>0</sub></td></tr><tr><td>Higher orders representation</td><td>∫<sub>f</sub> q ⋅ n x<sub>j</sub> x<sup>α<sub>j</sub></sup><br>∀f ∈ ∂K<br>∀l<sub>1</sub> < |α<sub>j</sub>| ≤ l<sub>2</sub>,<br>∀j ∈ ⌊1, d - 1⌋</td><td>∫<sub>f</sub> q ⋅ n x<sub>j</sub> x<sup>α<sub>j</sub></sup><br>∀f ∈ ∂K<br>l<sub>1</sub> < |α<sub>j</sub>| ≤ l<sub>2</sub>,<br>∀j ∈ ⌊1, d - 1⌋</td><td>∫<sub>∂K</sub> q ⋅ n p<sub>k</sub><br>∫<sub>f</sub> q ⋅ n x<sub>j</sub> x<sub>j</sub><sup>l<sub>2</sub></sup> x̃<br>∀p<sub>k</sub> ∈ H<sub>l<sub>2</sub></sub>(∂K) \ H<sub>l<sub>1</sub></sub>(∂K)<br>∀f ∈ ∂K<br>∀x̃ ∈ Q<sub>l<sub>2</sub></sub>(∂<sub>j</sub>K)</td><td>∫<sub>∂K</sub> q ⋅ n p<sub>k</sub><br>q(x<sub>jm</sub>) ⋅ n<br>∀p<sub>k</sub> ∈ H<sub>l<sub>2</sub></sub>(∂K) \ H<sub>l<sub>1</sub></sub>(∂K)<br>∀j ∈ ⌊1, n⌋,<br>{x<sub>jm</sub>}<sub>m</sub> sampling points</td></tr></tbody></table>

In two dimensions, the proposed setting reported above is fixed and all the mentioned degrees of freedom have to be considered. In higher dimensions however, the selection of degrees of freedom among the bold ones is a matter of taste, possibly directed by properties of the discretised quantities that are a priori known and that one wants represented down to the discrete level. Lastly, note also that in practice one could project on any other polynomial basis rather than using projections over the canonical one, substantially improving the conditioning of the characterisation matrix.
### 3.3.6 Two examples in two dimensions

Let us now detail examples of a discretisation framework that is contained in the previously presented setting. We start with a general framework for which a parallel with the Raviart–Thomas elements can be drawn from the order $k = 1$ on, and present in a second time an example of a reduced setting where a parallel with the Raviart–Thomas is achieved at any order. We then end the two-dimensional section on a generalising note.

#### 3.3.6.1 An example of a general setting

We consider a series of discretisation spaces by indexing the coefficients as $l_1 = 0, l_2 = k, m_1 = k - 1$ and $m_2 = k - 1$ for any $k \in \mathbb{N}$, seen here as the space order. The space $\mathbb{H}_k(K)$ is then defined as

$$ \mathbb{H}_k(K) = \{u \in H^1(K), u|_{\partial K} \in \mathcal{H}_0(\partial K), \Delta u \in \mathbb{Q}_{k-1}(K)\}^2 \\ \oplus x \{u \in H^1(K), u|_{\partial K} \in \mathcal{H}_k(\partial K), \Delta u \in \mathbb{Q}_{[k-1]}(K)\}. \quad (3.37) $$

By a straightforward application of the previous section, it comes

$$ \dim \mathbb{H}_k(K) = n(k+3) + 2k(k+1) - 1_{k>0}. \quad (3.38) $$

**Two-dimensional $H(\mathrm{div})$-conforming element** The examples of selected normal degrees of freedom defining the elements $E_k = (K, \mathbb{H}_k(K), \{\sigma\})$ presented in Section 3.3.4 reduce here to the expressions given in Table 3.2, here detailed only in their configuration “a” for the sake of concision.

<table><caption>**Table 3.2.** Example definitions of normal degrees of freedom in the 2D case, where $\gamma$ represent any polynomial parametrisation of the face $f$ of first order.</caption><thead><tr><th>Representation</th><th>of low order</th><th>of (A<sub>k</sub> ∩ B<sub>k</sub>)<sub>|∂K</sub></th><th>of higher orders</th></tr></thead><tbody><tr><td>Moments Ia</td><td>∫<sub>f</sub> q·n &nbsp;&nbsp;&nbsp;&nbsp; ∀f ∈ ∂K</td><td>∫<sub>f</sub> q<sub>i</sub> n<sub>i</sub> x<sub>i</sub> &nbsp;&nbsp;&nbsp;&nbsp; ∀i ∈ [1, 2],<br>&nbsp;&nbsp;&nbsp;&nbsp; ∀f ∈ ∂K</td><td>∫<sub>f</sub> q·n p<sub>k</sub> &nbsp;&nbsp;&nbsp;&nbsp; ∀f ∈ ∂K,<br>&nbsp;&nbsp;&nbsp;&nbsp; ∀p<sub>k</sub> ∈ Q<sub>k</sub>(f)<br>\Q<sub>0</sub>(f)</td></tr><tr><td>Moments IIa</td><td>∫<sub>f</sub> q<sub>i</sub> n<sub>i</sub> &nbsp;&nbsp;&nbsp;&nbsp; ∀f ∈ ∂K,<br>&nbsp;&nbsp;&nbsp;&nbsp; ∀i ∈ [1, 2]</td><td>∫<sub>f</sub> q·v p<sub>k</sub> &nbsp;&nbsp;&nbsp;&nbsp; ∀f ∈ ∂K,<br>&nbsp;&nbsp;&nbsp;&nbsp; ∀p<sub>k</sub> ∈ Q<sub>[1]</sub>(f)<br>given v≠n</td><td>∫<sub>f</sub> q·nγ(x) p<sub>k</sub> &nbsp;&nbsp;&nbsp;&nbsp; ∀f ∈ ∂K,<br>&nbsp;&nbsp;&nbsp;&nbsp; ∀p<sub>k</sub> ∈ Q<sub>k</sub>(f)<br>\Q<sub>0</sub>(f)</td></tr></tbody></table>

The internal degrees of freedom are simply set as

$$ \sigma(q) \mapsto \int_K q \cdot p_k \, dx, \quad \text{for all } p_k \in \mathcal{P}, \quad (3.39) $$

where $\mathcal{P}$ is chosen as the vectorial polynomial space

$$ \mathcal{P} = \mathbb{P}_{k, k-1} \times \mathbb{P}_{k-1, k} \setminus (\mathbb{P}_{[k], [k-1]} \times \mathbb{P}_{[k-1], [k]}) \cup \{(x, y)^T \mapsto \begin{pmatrix} x^k y^{k-1} \\ x^{k-1} y^k \end{pmatrix}\}. $$

Though the internal projection space is less refined than the one set on the edges, this is not bothersome as the impact of the divergence within the cell is less dramatic for the quality of the discretisation. Note also that in practice, the projections (3.39) can be performed on any basis of $\mathcal{P}$.
The designed element is thence $H(\mathrm{div})$-conforming by definition, and con- fers a similar discretisation framework for any cell shape, its expression being adaptive towards the cell geometry without changing the nature of the varia- tional discretisation, as depicted in Figure 3.6.

<figure><img src="image_42.png" /><figcaption>**Figure 3.6.** Adaptivity of the exemplified element towards the discretisation order and the cell shape. Note that same results are obtained for non-convex elements.</figcaption></figure>

### 3.3.6.2 An example of a reduced setting

As quickly addressed in Section 3.3.4.1 and as it will be shown in the numer- ical results, a classical construction of the space $\mathbb{H}_k(K)$ implies the degenera- tion of some normal basis functions into internal ones. This is a consequence of the coordinate-wise freedom provided on the boundary from the definition of the set $A_k$. Therefore, to allow a parallel with the Raviart – Thomas ele- ments from the lowest order on and to fulfil the optional Admissibility condi- tions 3.3.4, one can consider replacing the boundary conditions $u|_{\partial K} \in \mathcal{H}_{l_1}(\partial K)$ in $A_k$ to obtain the following reduced space.

**Definition 3.3.6 Reduced discretisation space**

$$ \mathbb{H}_k(K) = \{u \in H^1(K), u|_{\partial K} \equiv 1, \Delta u \in \mathbb{Q}_{m_1}(K)\}^d \oplus x \{u \in H^1(K), u|_{\partial K} \in \mathcal{H}_{l_2}(\partial K), \Delta u \in \mathbb{Q}_{[m_2]}(K)\} $$

There, the coordinate-wise freedom on the boundary is reduced and the nor- mal degrees of freedom can be set as in the classical Raviart – Thomas setting. Furthermore, contrarily to the general case, any definition of $l_2$, $m_1$ and $m_2$ leads to an $H(\mathrm{div})$-conformity ready space.

**An example of a reduced $H(\mathrm{div})$-conforming element** To emphasise the parallel with the Raviart – Thomas setting on the boundary, let us de- rive a two-dimensional example in the context of the reduced discretisation framework, similar in the spirit with the two examples presented in Sec- tion 3.3.6.1 for a general framework. Choosing in particular the coefficients generating the Virtual elements spaces accordingly to the sequence $l_2 = 1$, $m_1 = m_2 = k - 1$, where $k \in \mathbb{N}$ is seen as the space order, it comes

$$ \mathbb{H}_k(K) = \{u \in H^1(K), u|_{\partial K} \equiv 1, \Delta u \in \mathbb{Q}_{k-1}(K)\}^2 \oplus x \{u \in H^1(K), u|_{\partial K} \in \mathcal{H}_k(\partial K), \Delta u \in \mathbb{Q}_{[k-1]}(K)\} $$

whose dimension naturally reads $\dim \mathbb{H}_k(K) = n(k+1) + 2k(k-1) - 1_{k>0}$.
Therefore, exactly $k + 1$ normal functions per edge can be designed, fitting the framework of Raviart – Thomas. As this matches the dimension of $\mathbb{Q}_k(f)$, all the freedom is required to entirely determine the global normal component. Thus, as a straightforward reduction of the general case, the $H(\mathrm{div}, K)$-conforming element presented in Section 3.3.6 simplifies to the one defined through the normal degrees of freedom given in Table 3.3 and the internal degrees of freedom (3.39), where the projections can be considered against any basis of $\mathcal{P}$, not necessarily the canonical one.

<table><caption>**Table 3.3.** degrees of freedom of the element $Ia$ defined within the reduced setting.</caption><thead><tr><th>Core normal Moments</th><th>Misc moment</th><th>Core internal moment</th><th>Misc internal moment</th></tr></thead><tbody><tr><td>$\int_f p_k q \cdot n,$<br> $\forall f \in \partial K$<br> $\forall p_k \in \mathbb{Q}_k(f) \setminus \mathbb{Q}_0(f)$</td><td>$\int_f q \cdot n$<br> $\forall f \in \partial K$</td><td>$\int_K q \cdot \begin{pmatrix} x_1^l x_2^m \\ 0 \end{pmatrix}$ and $\int_K q \cdot \begin{pmatrix} 0 \\ x_1^m x_2^l \end{pmatrix}$<br> $\forall l \in [0, k], m \in [0, k-1]$<br> s.t. $(l, m) \neq (k, k-1)$</td><td>$\int_K q \cdot \begin{pmatrix} x_1^k x_2^{k-1} \\ x_1^{k-1} x_2^k \end{pmatrix}$</td></tr></tbody></table>

### 3.3.6.3 Many more elements

As examples, we only detailed here two declinations of two possible configuration of degrees of freedom. Needless to say, there exists many more possibilities, and their choice impact the properties of the elements, from its robustness to the properties of the associated basis functions.

In particular, we saw that it is possible to focus on a coordinate-wise boundary characterisation of the quantities living in a general space $\mathbb{H}_k(K)$ rather than focusing exclusively on their normal component in itself. Though not explicitly stated here, it is equally possible to shape the basis functions of $\mathbb{H}_k(K)$ without impacting the regularity, conformity, and order of the discretisation, simply by choosing a specific basis for the projectors $p$ involved in the normal degrees of freedom. On the same vein, one can also consider shifting the basis functions by modulating the offset directly from the definition of the degrees of freedom to enforce their positivity. More generally, plenty of discretisation focuses and degrees of freedom choices are possible, providing a great flexibility in the element design. However, the choices of the degrees of freedom and projector shapes have a very strong impact on the conditioning number of the linear systems characterising any discretised quantity. The examples presented here are the one offering the best compromise between the quality of the discretisation and the obtained conditioning. They are also the ones that are the closest to the classical Raviart–Thomas framework.

For interested readers, further details and extensive examples for both the general and the reduced space are available in [63], along with further investigations on various configurations and projector spaces.