# Stability estimates for solving Stokes problem with nonconforming finite elements.

Erell Jamelot

► To cite this version:

Erell Jamelot. Stability estimates for solving Stokes problem with nonconforming finite elements.. 2024. cea-03833616v4

HAL Id: cea-03833616

https://cea.hal.science/cea-03833616v4

Preprint submitted on 27 Mar 2024

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

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

<figure><img src="image_1.png" /></figure>
# Stability estimates for solving Stokes problem with nonconforming finite elements

Erell Jamelot

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

## To cite this version:

Erell Jamelot. Stability estimates for solving Stokes problem with nonconforming finite elements. Comptes Rendus. Mathématique, 2025, 363, pp.115-137. 10.5802/crmath.707. cea-03833616v5

HAL Id: cea-03833616

https://cea.hal.science/cea-03833616v5

Submitted on 27 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
STABILITY ESTIMATES FOR SOLVING STOKES PROBLEM
WITH NONCONFORMING FINITE ELEMENTS

ERELL JAMELOT

Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique
des Fluides, 91191 Gif-sur-Yvette, France

ABSTRACT. We propose to analyse the discretization of the Stokes problem with nonconforming finite elements in light of the T-coercivity. First we exhibit a family of operators to prove T-coercivity and we show that the stability constant is equal to the classical one up to a constant which depends on the Babuška-Aziz constant. Then we explicit the stability constants with respect to the shape regularity parameter for order 1 in 2 or 3 dimension, and order 2 in 2 dimension. In this last case, we improve the result of the original Fortin-Soulie paper. Second, we illustrate the importance of using a divergence-free velocity reconstruction on some numerical experiments.

**Keywords.** Stokes problem, T-coercitivity, Fortin-Soulie finite elements, Fortin operator

**2020 Mathematics Subject Classification.** 65N30, 35J57, 76D07

**Funding.** CEA SIMU/SITHY project

# 1. INTRODUCTION

The Stokes problem describes the steady state of incompressible Newtonian flows. It follows from the Navier–Stokes equations [1]. With regard to numerical analysis, the study of Stokes problem helps to build an appropriate approximation of the Navier–Stokes equations. We consider here a discretization with nonconforming finite elements [2, 3]. We propose to state the discrete inf-sup condition in light of the T-coercivity (cf. [4] for Helmholtz-like problems, see [5], [6] and [7] for the neutron diffusion equation), which allows to estimate the discrete error constant. In Section 2.1, we recall the T-coercivity theory [4], which is known to be an equivalent reformulation of the Banach–Nečas–Babuška Theorem and we apply it to the continuous Stokes problem. We give details on the triangulation in Section 2.4, and we apply the T-coercivity to the discretization of Stokes problem with nonconforming mixed finite elements in Section 2.5. For the Stokes problem, in the discrete case, this amounts to finding a Fortin operator. In Section 3, we precise the proof of the well-posedness in the case of order 1 and 2 nonconforming mixed finite elements. In Section 4.2, we illustrate the importance of using a divergence-free velocity on some numerical experiments.

E-mail address: erell.jamelot@cea.fr.

Date: March 27, 2024.
## 2. EXACT AND DISCRETE T-COERCIVITY FOR STOKES PROBLEM

2.1. **T-coercivity and application to Stokes problem.** We recall here the T-coercivity theory as written in [4]. Consider first the variational problem, where $V$ and $W$ are two Hilbert spaces and $f \in V'$:

$$ (2.1) \qquad \text{Find } u \in V \text{ such that } \forall v \in W, a(u, v) = \langle f, v \rangle_V. $$

Classically, we know that Problem (2.1) is well-posed if $a(\cdot, \cdot)$ satisfies the stability and the solvability conditions of the so-called Banach–Nečas–Babuška (BNB) Theorem (see e.g. [8, Thm. 25.9]). For some models, one can also prove the well-posedness using the T-coercivity theory (cf. [4] for Helmholtz-like problems, see [5], [6] and [7] for the neutron diffusion equation).

**Definition 1.** Let $V$ and $W$ be two Hilbert spaces and $a(\cdot, \cdot)$ be a continuous and bilinear form over $V \times W$. It is T-coercive if

$$ (2.2) \qquad \exists T \in \mathcal{L}(V, W), \text{ bijective, } \exists \alpha_T > 0, \forall v \in V, |a(v, Tv)| \ge \alpha_T \|v\|_V^2. $$

It is proved in [4, 9] that the T-coercivity condition is equivalent to the stability and solvability conditions of the BNB Theorem. Whereas the BNB theorem relies on an abstract inf-sup condition, T-coercivity uses explicit inf-sup operators, both at the continuous and discrete levels. Notice that if the pair $(T, \alpha_T)$ satisfies (2.2), then for all $\lambda > 0$, the pair $(T_\lambda, \alpha_{T_\lambda}) := (\lambda T, \lambda \alpha_T)$ also satisfies (2.2). Thus, there exist an infinity of pairs $(T, \alpha_T)$ and $(T_\lambda, \alpha_{T_\lambda})_{\lambda>0}$ satisfying (2.2).

**Theorem 1.** (well-posedness) Let $a(\cdot, \cdot)$ be a continuous and bilinear form. Suppose that the form $a(\cdot, \cdot)$ is T-coercive. Then Problem (2.1) is well-posed.

Let $\Omega$ be a connected bounded domain of $\mathbb{R}^d$, $d = 2, 3$, with a polygonal ($d=2$) or Lipschitz polyhedral ($d=3$) boundary $\partial\Omega$. We consider Stokes problem:

$$ (2.3) \qquad \text{Find } (u, p) \text{ such that } \left\{ \begin{array}{l} -\nu \Delta u + \text{grad } p = \mathbf{f}, \\ \text{div } u = 0. \end{array} \right. $$

with Dirichlet boundary conditions for the velocity $u$ and a normalization condition for the pressure $p$:

$$ u = 0 \text{ on } \partial\Omega, \quad \int_{\Omega} p = 0. $$

The vector field $u$ represents the velocity of the fluid and the scalar field $p$ represents its pressure divided by the fluid density which is supposed to be constant. The first equation of (2.3) corresponds to the momentum balance equation and the second one corresponds to the conservation of the mass. The constant parameter $\nu > 0$ is the kinematic viscosity of the fluid. The vector field $\mathbf{f} \in \mathbf{H}^{-1}(\Omega)$ represents a body forces divided by the fluid density.

Before stating the variational formulation of Problem (2.3), we provide some definition and reminders. Let us set $\mathbf{L}^2(\Omega) = (L^2(\Omega))^d$, $\mathbf{H}_0^1(\Omega) = (H_0^1(\Omega))^d$, $\mathbf{H}^{-1}(\Omega) = (H^{-1}(\Omega))^d$ its dual space and $L_{zmv}^2(\Omega) = \{q \in L^2(\Omega) | \int_{\Omega} q = 0\}$. We recall that $\mathbf{H}(\text{div}; \Omega) = \{v \in \mathbf{L}^2(\Omega) | \text{div } v \in L^2(\Omega)\}$. Let $h_\Omega$ be the diameter of $\Omega$. We recall Poincaré-Steklov inequality:

$$ (2.4) \qquad \exists C_{PS} > 0 |\forall v \in H_0^1(\Omega), \quad \|v\|_{L^2(\Omega)} \le C_{PS} h_\Omega \|\text{grad } v\|_{\mathbf{L}^2(\Omega)}. $$

Thanks to this result, in $H_0^1(\Omega)$, the semi-norm is equivalent to the natural norm, so that the scalar product reads $(v, w)_{H_0^1(\Omega)} = (\text{grad } v, \text{grad } w)_{\mathbf{L}^2(\Omega)}$ and the norm
the neutron diffusion equation), which allows to estimate the discrete error constant. In Section 2, we recall the T-coercivity theory [17], which is known to be an equivalent reformulation of the Banach-Nečas-Babuška Theorem and we apply it to the continuous Stokes problem. We give details on the triangulation, and we apply the T-coercivity to the discretization of Stokes problem with nonconforming mixed finite elements. For the Stokes problem, in the discrete case, this amounts to finding a Fortin operator. In Section 3, we precise the proof of the well-posedness in the case of order 1 and 2 nonconforming mixed finite elements. In Section 4, we illustrate the importance of using a divergence-free velocity on some numerical experiments.

## 2. Exact and discrete T-coercivity for Stokes problem

### 2.1. T-coercivity and application to Stokes problem

We recall here the T-coercivity theory as written in [17]. Consider first the variational problem, where $V$ and $W$ are two Hilbert spaces and $f \in V'$:

$$ \text{Find } u \in V \text{ such that } \forall v \in W, a(u, v) = \langle f, v \rangle_V. \quad (1) $$

Classically, we know that Problem (1) is well-posed if $a(\cdot, \cdot)$ satisfies the stability and the solvability conditions of the so-called Banach-Nečas-Babuška (BNB) Theorem (see e.g. [27, Theorem 25.9]). For some models, one can also prove the well-posedness using the T-coercivity theory (cf. [17] for Helmholtz-like problems, see [20, 33, 36] for the neutron diffusion equation).

**Definition 1.** Let $V$ and $W$ be two Hilbert spaces and $a(\cdot, \cdot)$ be a continuous bilinear form over $V \times W$. It is T-coercive if

$$ \exists T \in \mathcal{L}(V, W), \text{ bijective, } \exists \alpha_T > 0, \forall v \in V, |a(v, T v)| \ge \alpha_T \|v\|_V^2. \quad (2) $$

It is proved in [16, 17] that the T-coercivity condition is equivalent to the stability and solvability conditions of the BNB Theorem. Whereas the BNB Theorem relies on an abstract inf-sup condition, T-coercivity uses explicit inf-sup operators, both at the continuous and discrete levels. Notice that if the pair $(T, \alpha_T)$ satisfies (2), then for all $\lambda > 0$, the pair $(T_\lambda, \alpha_{T_\lambda}) := (\lambda T, \lambda \alpha_T)$ also satisfies (2). Thus, there exists an infinity of pairs $(T, \alpha_T)$ and $(T_\lambda, \alpha_{T_\lambda})_{\lambda>0}$ satisfying (2).

**Theorem 2 (well-posedness).** Let $a(\cdot, \cdot)$ be a continuous and bilinear form. Suppose that the form $a(\cdot, \cdot)$ is T-coercive. Then Problem (1) is well-posed.

Let $\Omega$ be a connected bounded domain of $\mathbb{R}^d$, $d = 2, 3$, with a polygonal ($d = 2$) or Lipschitz polyhedral ($d = 3$) boundary $\partial\Omega$. We consider the Stokes problem:

$$ \text{Find } (u, p) \text{ such that } -v \Delta u + \text{grad } p = f \quad \text{and} \quad \text{div } u = 0, \quad (3) $$

with Dirichlet boundary conditions for the velocity $u$ and a normalization condition for the pressure $p$:

$$ u = 0 \text{ on } \partial\Omega, \quad \int_{\Omega} p = 0. $$

The vector field $u$ represents the velocity of the fluid and the scalar field $p$ represents its pressure divided by the fluid density which is supposed to be constant. The first equation of (3) corresponds to the momentum balance equation and the second one corresponds to the conservation of the mass. The constant parameter $v > 0$ is the kinematic viscosity of the fluid. The vector field $f \in H^{-1}(\Omega)$ represents a body forces divided by the fluid density.

Before stating the variational formulation of Problem (3), we provide some definitions and reminders. Let us set $L^2(\Omega) = (L^2(\Omega))^d$, $H_0^1(\Omega) = (H_0^1(\Omega))^d$, $H^{-1}(\Omega) = (H^{-1}(\Omega))^d$ its dual space and
is $\|v\|_{H_0^1(\Omega)} = \|\text{grad } v\|_{L^2(\Omega)}$. Let $\mathbf{v}, \mathbf{w} \in \mathbf{H}_0^1(\Omega)$. We denote by $(v_i)_{i=1}^d$ (resp. $(w_i)_{i=1}^d$) the components of $\mathbf{v}$ (resp. $\mathbf{w}$), and we set $\text{Grad } \mathbf{v} = (\partial_j v_i)_{i,j=1}^d \in L^2(\Omega)$, where $L^2(\Omega) = [L^2(\Omega)]^{d \times d}$. We have:

$$ (\text{Grad } \mathbf{v}, \text{Grad } \mathbf{w})_{L^2(\Omega)} = (\mathbf{v}, \mathbf{w})_{\mathbf{H}_0^1(\Omega)} = \sum_{i=1}^{d} (v_i, w_i)_{\mathbf{H}_0^1(\Omega)} $$

and:

$$ \|\mathbf{v}\|_{\mathbf{H}_0^1(\Omega)} = \left( \sum_{j=1}^{d} \|v_j\|_{\mathbf{H}_0^1(\Omega)}^2 \right)^{\frac{1}{2}} = \|\text{Grad } \mathbf{v}\|_{L^2(\Omega)}. $$

Let us set $\mathbf{V} = \{\mathbf{v} \in \mathbf{H}_0^1(\Omega) | \text{div } \mathbf{v} = 0\}$. The vector space $\mathbf{V}$ is a closed subset of $\mathbf{H}_0^1(\Omega)$. We denote by $\mathbf{V}^{\perp}$ the orthogonal of $\mathbf{V}$ in $\mathbf{H}_0^1(\Omega)$. We recall that [1, cor. I.2.4]:

**Proposition 1.** The operator $\text{div} : \mathbf{H}_0^1(\Omega) \to L^2(\Omega)$ is an isomorphism of $\mathbf{V}^{\perp}$ onto $L_{zmv}^2(\Omega)$. We call $C_{\text{div}}$ the constant such that:

$$ (2.5) \qquad \forall p \in L_{zmv}^2(\Omega), \exists! \mathbf{v} \in \mathbf{V}^{\perp} \mid \text{div } \mathbf{v} = p \text{ and } \|\mathbf{v}\|_{\mathbf{H}_0^1(\Omega)} \le C_{\text{div}} \|p\|_{L^2(\Omega)} $$

The constant $C_{\text{div}}$ depends only on the domain $\Omega$.

Recall that we have: $\|\mathbf{v}\|_{\mathbf{H}_0^1(\Omega)}^2 = \|\text{curl } \mathbf{v}\|_{L^2(\Omega)}^2 + \|\text{div } \mathbf{v}\|_{L^2(\Omega)}^2 \ge \|p\|_{L^2(\Omega)}^2$. Hence $C_{\text{div}} \ge 1$. Actually, the constant $C_{\text{div}}$ is such that $C_{\text{div}} = 1/\beta(\Omega)$ where $\beta(\Omega)$, known as the Babuška-Aziz constant, is the inf-sup condition<fnref n="1" />:

$$ (2.6) \qquad \beta(\Omega) = \inf_{q \in L_{zmv}^2(\Omega) \setminus \{0\}} \sup_{\mathbf{v} \in \mathbf{H}_0^1(\Omega) \setminus \{0\}} \frac{(q, \text{div } \mathbf{v})_{L^2(\Omega)}}{\|q\|_{L^2(\Omega)} \|\mathbf{v}\|_{\mathbf{H}_0^1(\Omega)}}. $$

Generally, the value of $\beta(\Omega)$ is not known explicitly. In [10], Bernardi et al established results on the discrete approximation of $\beta(\Omega)$ using conforming finite elements. Recently, Gallistl proposed in [11] a numerical scheme with adaptive meshes for computing approximations to $\beta(\Omega)$. In the case of $d = 2$, Costabel and Dauge [12] established the following bound:

**Theorem 2.** Let $\Omega \subset \mathbb{R}^2$ be a domain contained in a ball of radius $R$, star-shaped with respect to a concentric ball of radius $\rho$. Then

$$ (2.7) \qquad \beta(\Omega) \ge \frac{\rho}{\sqrt{2} R} \left( 1 + \sqrt{1 - \frac{\rho^2}{R^2}} \right)^{-\frac{1}{2}} \ge \frac{\rho}{2 R}. $$

Let us detail the bound for some remarkable domains. If $\Omega$ is a ball, $\beta(\Omega) \ge \frac{1}{2}$ and if $\Omega$ is a square, $\beta(\Omega) \ge \frac{1}{2\sqrt{2}}$. Suppose now that $\Omega$ is stretched in some direction by a factor $k$, then $\beta(\Omega) \ge \frac{1}{2k}$. Finally, if $\Omega$ is L-shaped (resp. cross-shaped) such that $L = kl$, where $L$ is the largest length and $l$ is the smallest length of an edge, then $\beta(\Omega) \ge \frac{1}{2\sqrt{2}k}$ (resp. $\beta(\Omega) \ge \frac{1}{4k}$).

The variational formulation of Problem (2.3) reads:

Find $(\mathbf{u}, p) \in \mathbf{H}_0^1(\Omega) \times L_{zmv}^2(\Omega)$ such that

$$ (2.8) \qquad \left\{ \begin{array}{ll} \nu(\mathbf{u}, \mathbf{v})_{\mathbf{H}_0^1(\Omega)} - (p, \text{div } \mathbf{v})_{L^2(\Omega)} & = \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}_0^1(\Omega)} \quad \forall \mathbf{v} \in \mathbf{H}_0^1(\Omega); \\ (q, \text{div } \mathbf{u})_{L^2(\Omega)} & = 0 \quad \forall q \in L_{zmv}^2(\Omega). \end{array} \right. $$

<footnote n="1">¹or Ladyzhenskaya-Babuška-Brezzi condition</footnote>
$L^2_{zmv}(\Omega) = \{q \in L^2(\Omega) | \int_{\Omega} q = 0\}$. We recall that $\mathbf{H}(\text{div}; \Omega) = \{\mathbf{v} \in L^2(\Omega) | \text{div}\mathbf{v} \in L^2(\Omega)\}$. Let $h_{\Omega}$ be the diameter of $\Omega$. We recall the Poincaré-Steklov inequality [26, Lemma 3.24]:

$$ \exists C_{PS} > 0 \text{ such that } \forall v \in H_0^1(\Omega), \quad \|v\|_{L^2(\Omega)} \le C_{PS} h_{\Omega} \|\text{grad } v\|_{L^2(\Omega)}. \quad (4) $$

Thanks to this result, in $H_0^1(\Omega)$, the semi-norm is equivalent to the natural norm, so that the scalar product reads $(v, w)_{H_0^1(\Omega)} = (\text{grad } v, \text{grad } w)_{L^2(\Omega)}$ and the norm is $\|v\|_{H_0^1(\Omega)} = \|\text{grad } v\|_{L^2(\Omega)}$. Let $\mathbf{v}, \mathbf{w} \in H_0^1(\Omega)$. We denote by $(v_i)_{i \neq 1}^d$ (resp. $(w_i)_{i \neq 1}^d$) the components of $\mathbf{v}$ (resp. $\mathbf{w}$), and we set $\text{Grad } \mathbf{v} = (\partial_j v_i)_{i,j=1}^d \in L^2(\Omega)$, where $L^2(\Omega) = (L^2(\Omega))^{d \times d}$. We have:

$$ (\text{Grad } \mathbf{v}, \text{Grad } \mathbf{w})_{L^2(\Omega)} = (\mathbf{v}, \mathbf{w})_{H_0^1(\Omega)} = \sum_{i=1}^{d} (v_i, w_i)_{H_0^1(\Omega)} $$

and:

$$ \|\mathbf{v}\|_{H_0^1(\Omega)} = \left( \sum_{j=1}^{d} \|v_j\|_{H_0^1(\Omega)}^2 \right)^{\frac{1}{2}} = \|\text{Grad } \mathbf{v}\|_{L^2(\Omega)}. $$

Let us set $\mathbf{V} = \{\mathbf{v} \in H_0^1(\Omega) | \text{div}\mathbf{v} = 0\}$. The vector space $\mathbf{V}$ is a closed subset of $H_0^1(\Omega)$. We denote by $\mathbf{V}^{\perp}$ the orthogonal component of $\mathbf{V}$ in $H_0^1(\Omega)$. We recall the following result.

**Proposition 3 ([32, Corollary I.2.4])**. The operator $\text{div}: H_0^1(\Omega) \to L^2(\Omega)$ is an isomorphism of $\mathbf{V}^{\perp}$ onto $L^2_{zmv}(\Omega)$. We call $C_{\text{div}}$ the constant such that:

$$ \forall p \in L^2_{zmv}(\Omega), \exists! \mathbf{v} \in \mathbf{V}^{\perp} \text{ such that } \text{div}\mathbf{v} = p \text{ and } \|\mathbf{v}\|_{H_0^1(\Omega)} \le C_{\text{div}} \|p\|_{L^2(\Omega)}. \quad (5) $$

The constant $C_{\text{div}}$ depends only on the domain $\Omega$. Recall that we have:

$$ \|\mathbf{v}\|_{H_0^1(\Omega)}^2 = \|\text{curl } \mathbf{v}\|_{L^2(\Omega)}^2 + \|\text{div } \mathbf{v}\|_{L^2(\Omega)}^2 \ge \|p\|_{L^2(\Omega)}^2, $$

hence $C_{\text{div}} \ge 1$. Using Proposition 3, we can define a lifting in $H^1(\Omega)$ of a function $v \in L^2(\Omega)$:

**Corollary 4.** For all $v \in L^2(\Omega)$, there exists $\mathbf{s} \in H^1(\Omega)$ such that:

$$ \text{div}\mathbf{s} = v \quad \text{and} \quad \|\mathbf{s}\|_{L^2(\Omega)} + h_{\Omega} \|\text{Grad}\mathbf{s}\|_{L^2(\Omega)} \le C_{\Omega} h_{\Omega} \|v\|_{L^2(\Omega)}, \quad (6) $$

where the dimensionless constant $C_{\Omega}$ depends on $C_{\text{div}}$ and $C_{PS}$.

**Proof.** Let $\underline{v} = \int_{\Omega} \underline{v} / |\Omega|$, and $\mathbf{s}_0 \in H_0^1(\Omega)$ be such that $\text{div}\mathbf{s}_0 = v - \underline{v}$ and $\|\text{Grad}\mathbf{s}_0\|_{L^2(\Omega)} \le C_{\text{div}} \|v\|_{L^2(\Omega)}$ (cf. Proposition 3). Let $d' \in \{1, \dots, d\}$. We consider $\mathbf{s} := \mathbf{s}_0 + \underline{v}(x_{d'} - \underline{x}_{d'})\mathbf{e}_{d'} \in H^1(\Omega)$, where $\underline{x}_{d'} = \int_{\Omega} x_{d'} / |\Omega|$. We have: $\text{div}\mathbf{s} = v$ and $\|\text{Grad}\mathbf{s}\|_{L^2(\Omega)} \le \tilde{C}_{\text{div}} \|v\|_{L^2(\Omega)}$, where $\tilde{C}_{\text{div}} = C_{\text{div}} + 1$. Using inequality (4), one can prove that $\|\mathbf{s}\|_{L^2(\Omega)} \le h_{\Omega} \tilde{C}_{\Omega} \|v\|_{L^2(\Omega)}$, where $\tilde{C}_{\Omega} = C_{PS} C_{\text{div}} + 1$. Setting $C_{\Omega} = \tilde{C}_{\Omega} + \tilde{C}_{\text{div}}$, we obtain (6). $\square$

Actually, the constant $C_{\text{div}}$ is such that $C_{\text{div}} = 1/\beta(\Omega)$ where $\beta(\Omega)$, known as the Babuška-Aziz constant, is the inf-sup condition (or Ladyzhenskaya-Babuška-Brezzi condition):

$$ \beta(\Omega) = \inf_{q \in L^2_{zmv}(\Omega) \setminus \{0\}} \sup_{\mathbf{v} \in H_0^1(\Omega) \setminus \{0\}} \frac{(q, \text{div}\mathbf{v})_{L^2(\Omega)}}{\|q\|_{L^2(\Omega)} \|\mathbf{v}\|_{H_0^1(\Omega)}}. \quad (7) $$

Generally, the value of $\beta(\Omega)$ is not known explicitly. In [6], Bernardi et al established results on the discrete approximation of $\beta(\Omega)$ using conforming finite elements. Recently, Gallistl proposed in [30] a numerical scheme with adaptive meshes for computing approximations to $\beta(\Omega)$. In the case of $d=2$, Costabel and Dauge [21] established the following bound:

**Theorem 5 ([21, Theorem 2.3])**. Let $\Omega \subset \mathbb{R}^2$ be a domain contained in a ball of radius $R$, star-shaped with respect to a concentric ball of radius $\rho$. Then

$$ \beta(\Omega) \ge \frac{\rho}{\sqrt{2R}} \left( 1 + \sqrt{1 - \frac{\rho^2}{R^2}} \right)^{-\frac{1}{2}} \ge \frac{\rho}{2R}. \quad (8) $$
Classically, one proves that Problem (2.8) is well-posed using Poincaré-Steklov inequality (2.4) and Prop. 1. Check for instance the proof of [1, Thm. I.5.1].

Let us set $\mathcal{X} = \mathbf{H}_0^1(\Omega) \times L_{zmv}^2(\Omega)$ which is a Hilbert space which we endow with the following norm:

$$ (2.9) \qquad \|(\mathbf{v}, q)\|_{\mathcal{X}} = \left( \|\mathbf{v}\|_{\mathbf{H}_0^1(\Omega)}^2 + \nu^{-2} \|q\|_{L^2(\Omega)}^2 \right)^{\frac{1}{2}}. $$

We consider now the following bilinear symmetric and continuous form:

$$ (2.10) \qquad \left\{ \begin{array}{l} a_S : \mathcal{X} \times \mathcal{X} \to \mathbb{R} \\ (\mathbf{u}', p') \times (\mathbf{v}, q) \mapsto \nu(\mathbf{u}', \mathbf{v})_{\mathbf{H}_0^1(\Omega)} \\ \qquad - (p', \operatorname{div} \mathbf{v})_{L^2(\Omega)} - (q, \operatorname{div} \mathbf{u}')_{L^2(\Omega)} \end{array} \right. . $$

We can write Problem (2.3) in an equivalent way as follows:

$$ (2.11) \quad \text{Find } (\mathbf{u}, p) \in \mathcal{X} \text{ such that } a_S((\mathbf{u}, p), (\mathbf{v}, q)) = \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}_0^1(\Omega)} \quad \forall (\mathbf{v}, q) \in \mathcal{X}. $$

Let us prove that Problem (2.11) is well-posed using the T-coercivity theory.

**Proposition 2.** The bilinear form $a_S(\cdot, \cdot)$ is T-coercive:

$$ (2.12) \qquad \begin{gathered} \exists T \in \mathcal{L}(\mathcal{X}), \text{ bijective, } \exists \alpha_T > 0, \forall (\mathbf{u}', p') \in \mathcal{X}, \\ a_S((\mathbf{u}', p'), T((\mathbf{u}', p'))) \ge \alpha_T \|(\mathbf{u}', p')\|_{\mathcal{X}}^2. \end{gathered} $$

*Proof.* We follow here the proof given in [13, 14, 15]. Let us consider $(\mathbf{u}', p') \in \mathcal{X}$ and let us build $(\mathbf{v}^*, q^*) = T(\mathbf{u}', p') \in \mathcal{X}$ satisfying (2.2) (with $V = \mathcal{X}$). We need three main steps.

1. According to Prop. 1, it exists $\tilde{\mathbf{v}}_{p'} \in \mathbf{H}_0^1(\Omega)$ such that: $\operatorname{div} \tilde{\mathbf{v}}_{p'} = p'$ in $\Omega$ and $\|\tilde{\mathbf{v}}_{p'}\|_{\mathbf{H}_0^1(\Omega)} \le C_{\operatorname{div}} \|p'\|_{L^2(\Omega)}$. Let us set $\mathbf{v}_{p'} = \nu^{-1} \tilde{\mathbf{v}}_{p'}$ so that $\operatorname{div} \mathbf{v}_{p'} = \nu^{-1} p'$ and

$$ (2.13) \qquad \|\mathbf{v}_{p'}\|_{\mathbf{H}_0^1(\Omega)} \le \nu^{-1} C_{\operatorname{div}} \|p'\|_{L^2(\Omega)}. $$

Let us set $(\mathbf{v}^*, q^*) := (\gamma \mathbf{u}' - \mathbf{v}_{p'}, -\gamma p')$, with $\gamma > 0$. We obtain:

$$ (2.14) \qquad a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) = \nu \gamma \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + \nu^{-1} \|p'\|_{L^2(\Omega)}^2 - \nu (\mathbf{u}', \mathbf{v}_{p'})_{\mathbf{H}_0^1(\Omega)}. $$

2. In order to bound the last term of (2.14), we use Young inequality and then inequality (2.13), so that for all $\eta > 0$:

$$ (2.15) \qquad (\mathbf{u}', \mathbf{v}_{p'})_{\mathbf{H}_0^1(\Omega)} \le \frac{\eta}{2} \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + \frac{\eta^{-1}}{2} \left( \frac{C_{\operatorname{div}}}{\nu} \right)^2 \|p'\|_{L^2(\Omega)}^2. $$

3. Using the bound (2.15) in (2.14) and choosing $\eta = \gamma$, we get:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \nu \left( \frac{\gamma}{2} \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + \nu^{-2} \left( 1 - \frac{\gamma^{-1}}{2} (C_{\operatorname{div}})^2 \right) \|p'\|_{L^2(\Omega)}^2 \right). $$

Consider now $\gamma = (C_{\operatorname{div}})^2$. We obtain:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \frac{1}{2} \nu \left( (C_{\operatorname{div}})^2 \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + \nu^{-2} \|p'\|_{L^2(\Omega)}^2 \right). $$

Reminding that $C_{\operatorname{div}} \ge 1$, it comes:

$$ (2.16) \qquad a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \frac{\nu}{2} \|(\mathbf{u}', p')\|_{\mathcal{X}}^2. $$
Let us detail the bound for some remarkable domains. If $\Omega$ is a ball, $\beta(\Omega) \ge \frac{1}{2}$ and if $\Omega$ is a square, $\beta(\Omega) \ge \frac{1}{2\sqrt{2}}$. Suppose now that $\Omega$ is stretched in some direction by a factor $\kappa$, then $\beta(\Omega) \ge \frac{1}{2\kappa}$. Finally, if $\Omega$ is L-shaped (resp. cross-shaped) such that $L = \kappa l$, where $L$ is the largest length and $l$ is the smallest length of an edge, then $\beta(\Omega) \ge \frac{1}{2\sqrt{2\kappa}}$ (resp. $\beta(\Omega) \ge \frac{1}{4\kappa}$).

The variational formulation of Problem (3) reads:

$$ egin{cases} 
\text{Find } (\mathbf{u}, p) \in \mathbf{H}_0^1(\Omega) \times L_{zmv}^2(\Omega) \text{ such that} \\ 
\|\mathbf{v}(\mathbf{u}, \mathbf{v})\|_{\mathbf{H}_0^1(\Omega)} - (p, \text{div}\mathbf{v})_{L^2(\Omega)} = \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}_0^1(\Omega)} \quad \forall \mathbf{v} \in \mathbf{H}_0^1(\Omega); \\ 
(q, \text{div}\mathbf{u})_{L^2(\Omega)} = 0 \quad \forall q \in L_{zmv}^2(\Omega). 
\end{cases} 
\quad (9) $$

Classically, one proves that Problem (9) is well-posed using the Poincaré-Steklov inequality (4) and Proposition 3. Check for instance the proof of [32, Theorem I.5.1].

Let us set $\mathcal{X} = \mathbf{H}_0^1(\Omega) \times L_{zmv}^2(\Omega)$ which is a Hilbert space which we endow with the following norm:

$$ \| (\mathbf{v}, q) \|_{\mathcal{X}} = \left( \| \mathbf{v} \|_{\mathbf{H}_0^1(\Omega)}^2 + v^{-2} \| q \|_{L^2(\Omega)}^2 \right)^{\frac{1}{2}}. \quad (10) $$

We consider now the following bilinear symmetric and continuous form:

$$ \left\{ \begin{array}{l} a_S : \mathcal{X} \times \mathcal{X} \to \mathbb{R} \\ (\mathbf{u}', p') \times (\mathbf{v}, q) \mapsto v(\mathbf{u}', \mathbf{v})_{\mathbf{H}_0^1(\Omega)} - (p', \text{div}\mathbf{v})_{L^2(\Omega)} - (q, \text{div}\mathbf{u}')_{L^2(\Omega)}. \end{array} \right. \quad (11) $$

We can write Problem (3) in an equivalent way as follows:

$$ \text{Find } (\mathbf{u}, p) \in \mathcal{X} \text{ such that } a_S((\mathbf{u}, p), (\mathbf{v}, q)) = \langle \mathbf{f}, \mathbf{v} \rangle_{\mathbf{H}_0^1(\Omega)} \quad \forall (\mathbf{v}, q) \in \mathcal{X}. \quad (12) $$

Let us prove that Problem (12) is well-posed using the T-coercivity theory.

**Proposition 6.** *The bilinear form $a_S(\cdot, \cdot)$ is T-coercive:*

$$ \exists T \in \mathcal{L}(\mathcal{X}), \text{ bijective, } \exists \alpha_T > 0, \forall (\mathbf{u}', p') \in \mathcal{X}, \quad (13) $$

$$ a_S((\mathbf{u}', p'), T((\mathbf{u}', p'))) \ge \alpha_T \| (\mathbf{u}', p') \|_{\mathcal{X}}^2. $$

**Proof.** We follow here the proof given in [5, 9, 19]. Let us consider $(\mathbf{u}', p') \in \mathcal{X}$ and let us build $(\mathbf{v}^*, q^*) = T(\mathbf{u}', p') \in \mathcal{X}$ satisfying (2) (with $V = \mathcal{X}$). We need three main steps.

**Step 1.** According to Proposition 3, there exists $\tilde{\mathbf{v}}_{p'} \in \mathbf{H}_0^1(\Omega)$ such that: $\text{div}\tilde{\mathbf{v}}_{p'} = p'$ in $\Omega$ and $\|\tilde{\mathbf{v}}_{p'}\|_{\mathbf{H}_0^1(\Omega)} \le C_{\text{div}} \|p'\|_{L^2(\Omega)}$. Let us set $\mathbf{v}_{p'} = v^{-1}\tilde{\mathbf{v}}_{p'}$ so that $\text{div}\mathbf{v}_{p'} = v^{-1}p'$ and

$$ \| \mathbf{v}_{p'} \|_{\mathbf{H}_0^1(\Omega)} \le v^{-1} C_{\text{div}} \| p' \|_{L^2(\Omega)}. \quad (14) $$

Let us set $(\mathbf{v}^*, q^*) := (\gamma \mathbf{u}' - \mathbf{v}_{p'}, -\gamma p')$, with $\gamma > 0$. We obtain:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) = v \gamma \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + v^{-1} \|p'\|_{L^2(\Omega)}^2 - v (\mathbf{u}', \mathbf{v}_{p'})_{\mathbf{H}_0^1(\Omega)}. \quad (15) $$

**Step 2.** In order to bound the last term of (15), we use Young's inequality and then inequality (14), so that for all $\eta > 0$:

$$ (\mathbf{u}', \mathbf{v}_{p'})_{\mathbf{H}_0^1(\Omega)} \le \frac{\eta}{2} \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + \frac{\eta^{-1}}{2} \left( \frac{C_{\text{div}}}{v} \right)^2 \|p'\|_{L^2(\Omega)}^2. \quad (16) $$

**Step 3.** Using the bound (16) in (15) and choosing $\eta = \gamma$, we get:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge v \left( \frac{\gamma}{2} \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + v^{-2} \left( 1 - \frac{\gamma^{-1}}{2} (C_{\text{div}})^2 \right) \|p'\|_{L^2(\Omega)}^2 \right). $$

Consider now $\gamma = (C_{\text{div}})^2$. We obtain:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \frac{1}{2} v \left( (C_{\text{div}})^2 \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + v^{-2} \|p'\|_{L^2(\Omega)}^2 \right). $$

Reminding that $C_{\text{div}} \ge 1$, it comes:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \frac{v}{2} \|(\mathbf{u}', p')\|_{\mathcal{X}}^2. \quad (17) $$
We obtain (2.12) with $\alpha_T = \frac{\nu}{2}$.

The operator $T$ such that $T((\mathbf{u}', p')) = (\mathbf{v}^*, q^*)$ is linear and continuous. We have indeed:

$$ \begin{aligned} \|T((\mathbf{u}', p'))\|_{\mathcal{X}}^2 &:= \|\mathbf{v}^*\|_{\mathbf{H}_0^1(\Omega)}^2 + \nu^{-2} \|q^*\|_{L^2(\Omega)}^2 \\ &\le 2\gamma^2 \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + 2\|\mathbf{v}_{p'}\|_{\mathbf{H}_0^1(\Omega)}^2 + \gamma^2 \nu^{-2} \|p'\|_{L^2(\Omega)}^2, \\ &\le 2\gamma^2 \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + (2(C_{\text{div}})^2 + \gamma^2) \nu^{-2} \|p'\|_{L^2(\Omega)}^2. \end{aligned} $$

We deduce that:

$$ (2.17) \qquad \|T((\mathbf{u}', p'))\|_{\mathcal{X}} \le C_{\max} \|(\mathbf{u}', p')\|_{\mathcal{X}}, $$

where $C_{\max} = C_{\text{div}}(\max(2 + (C_{\text{div}})^2, 2(C_{\text{div}})^2))^{\frac{1}{2}}$.

Remark that, given $(\mathbf{v}^*, q^*) \in \mathcal{X}$, choosing $(\mathbf{u}', p') = (\gamma^{-1}\mathbf{v}^* - \gamma^{-2}\mathbf{v}_{q^*}, -\gamma^{-1}q^*)$ yields $T((\mathbf{u}', p')) = (\mathbf{v}^*, q^*)$. Hence, the operator $T \in \mathcal{L}(\mathcal{X})$ is bijective. $\square$

We can now prove the

**Theorem 3.** *Problem (2.11) is well-posed. It admits one and only one solution such that:*

$$ (2.18) \qquad \forall \mathbf{f} \in \mathbf{H}^{-1}(\Omega), \quad \left\{ \begin{array}{l} \|\mathbf{u}\|_{\mathbf{H}_0^1(\Omega)} \le \nu^{-1} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}, \\ \|\mathbf{p}\|_{L^2(\Omega)} \le C_{\text{div}} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}. \end{array} \right. $$

*Proof.* According to Prop. 2, the continuous bilinear form $a_S(\cdot, \cdot)$ is T-coercive. Hence, according to Theorem 1, Problem (2.11) is well-posed. Let us prove (2.18). Consider $(\mathbf{u}, p)$ the unique solution of Problem (2.11). Choosing $\mathbf{v} = 0$, we obtain that $\forall q \in L_{zmv}^2(\Omega)$, $(q, \text{div } \mathbf{u})_{L^2(\Omega)} = 0$, so that $\mathbf{u} \in \mathbf{V}$. Now, choosing $\mathbf{v} = \mathbf{u}$ and using Cauchy-Schwarz inequality, we have: $\nu \|\mathbf{u}\|_{\mathbf{H}_0^1(\Omega)}^2 = \langle \mathbf{f}, \mathbf{u} \rangle_{\mathbf{H}_0^1(\Omega)} \le \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \|\mathbf{u}\|_{\mathbf{H}_0^1(\Omega)}$, so that: $\|\mathbf{u}\|_{\mathbf{H}_0^1(\Omega)} \le \nu^{-1} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}$. Next, we choose in (2.11) $\mathbf{v} = \tilde{\mathbf{v}}_p \in \mathbf{V}^{\perp}$, where $\text{div } \tilde{\mathbf{v}}_p = -p$ (see Prop. 1). Since $\mathbf{u} \in \mathbf{V}$ and $\tilde{\mathbf{v}}_p \in \mathbf{V}^{\perp}$, we have $(\mathbf{u}, \tilde{\mathbf{v}}_p)_{\mathbf{H}_0^1(\Omega)} = 0$. This gives:

$$ \begin{aligned} -(p, \text{div } \tilde{\mathbf{v}}_p)_{L^2(\Omega)} &= \|p\|_{L^2(\Omega)}^2 = \langle \mathbf{f}, \tilde{\mathbf{v}}_p \rangle_{\mathbf{H}_0^1(\Omega)} \\ &\le \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \|\tilde{\mathbf{v}}_p\|_{\mathbf{H}_0^1(\Omega)} \le C_{\text{div}} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \|p\|_{L^2(\Omega)}, \end{aligned} $$

so that: $\|p\|_{L^2(\Omega)} \le C_{\text{div}} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}$. $\square$

2.2. **Comments on the stability constant.** Using (2.17) in (2.16), we have:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \alpha_T (C_{\max})^{-1} \|(\mathbf{u}', p')\|_{\mathcal{X}} \|(\mathbf{v}^*, q^*)\|_{\mathcal{X}} $$

Let us set $\|\|T\|\| := \sup_{((\mathbf{u}', p') \in \mathcal{X} \setminus \{\mathbf{0}, 0\})} \frac{|T((\mathbf{u}', p'))|}{\|(\mathbf{u}', p')\|_{\mathcal{X}}}$.

bound: $\|\|T\|\| \le C_{\max}$.

Hence, we have:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \frac{\alpha_T}{\|\|T\|\|} \|(\mathbf{u}', p')\|_{\mathcal{X}} \|(\mathbf{v}^*, q^*)\|_{\mathcal{X}} $$

We recover the first Banach-Nečas-Babuška condition [8, Thm. 25.9, (BNB1)]. Thus, the T-coercivity approach gives an overestimate of the stability constant $\alpha$.
We obtain (13) with $\alpha_T = \frac{\nu}{2}$.

The operator $T$ such that $T((\mathbf{u}', p')) = (\mathbf{v}^*, q^*)$ is linear and continuous. We have indeed:

$$
\begin{aligned}
\|T((\mathbf{u}', p'))\|_{\mathcal{X}}^2 &:= \|\mathbf{v}^*\|_{\mathbf{H}_0^1(\Omega)}^2 + \nu^{-2} \|q^*\|_{L^2(\Omega)}^2 \\
&\le 2\gamma^2 \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + 2\|\mathbf{v}p'\|_{\mathbf{H}_0^1(\Omega)}^2 + \gamma^2 \nu^{-2} \|p'\|_{L^2(\Omega)}^2, \\
&\le 2\gamma^2 \|\mathbf{u}'\|_{\mathbf{H}_0^1(\Omega)}^2 + (2(C_{\text{div}})^2 + \gamma^2) \nu^{-2} \|p'\|_{L^2(\Omega)}^2.
\end{aligned}
$$

We deduce that:

$$ \|T((\mathbf{u}', p'))\|_{\mathcal{X}} \le C_{\max} \|(\mathbf{u}', p')\|_{\mathcal{X}} \text{ where } C_{\max} = C_{\text{div}} \left(\max(2 + (C_{\text{div}})^2, 2(C_{\text{div}})^2)\right)^{\frac{1}{2}}. \quad (18) $$

Remark that, given $(\mathbf{v}^*, q^*) \in \mathcal{X}$, choosing $(\mathbf{u}', p') = (\gamma^{-1}\mathbf{v}^* - \gamma^{-2}\mathbf{v}_{q^*}, -\gamma^{-1}q^*)$ yields $T((\mathbf{u}', p')) = (\mathbf{v}^*, q^*)$. Hence, the operator $T \in \mathcal{L}(\mathcal{X})$ is bijective. $\square$

We can now prove the following result.

**Theorem 7.** *Problem (12) is well-posed. It admits one and only one solution such that:*

$$ \forall \mathbf{f} \in \mathbf{H}^{-1}(\Omega), \quad \begin{cases} \|\mathbf{u}\|_{\mathbf{H}_0^1(\Omega)} \le \nu^{-1} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}, \\ \|p\|_{L^2(\Omega)} \le C_{\text{div}} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}. \end{cases} \quad (19) $$

**Proof.** According to Proposition 6, the continuous bilinear form $a_S(\cdot, \cdot)$ is $T$-coercive. Hence, according to Theorem 2, Problem (12) is well-posed. Let us prove (19). Consider $(\mathbf{u}, p)$ the unique solution of Problem (12). Choosing $\mathbf{v} = 0$, we obtain that $\forall q \in L_{zmv}^2(\Omega)$, $(q, \text{div}\mathbf{u})_{L^2(\Omega)} = 0$, so that $\mathbf{u} \in \mathbf{V}$. Now, choosing $\mathbf{v} = \mathbf{u}$ and using Cauchy-Schwarz inequality, we have: $\nu \|\mathbf{u}\|_{\mathbf{H}_0^1(\Omega)}^2 = \langle \mathbf{f}, \mathbf{u} \rangle_{\mathbf{H}_0^1(\Omega)} \le \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \|\mathbf{u}\|_{\mathbf{H}_0^1(\Omega)}$, so that: $\|\mathbf{u}\|_{\mathbf{H}_0^1(\Omega)} \le \nu^{-1} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}$. Next, we choose in (12) $\mathbf{v} = \tilde{\mathbf{v}}_p \in \mathbf{V}^{\perp}$, where $\text{div}\tilde{\mathbf{v}}_p = -p$ (see Proposition 3). Since $\mathbf{u} \in \mathbf{V}$ and $\tilde{\mathbf{v}}_p \in \mathbf{V}^{\perp}$, we have $(\mathbf{u}, \tilde{\mathbf{v}}_p)_{\mathbf{H}_0^1(\Omega)} = 0$. This gives:

$$ -(p, \text{div}\tilde{\mathbf{v}}_p)_{L^2(\Omega)} = \|p\|_{L^2(\Omega)}^2 = \langle \mathbf{f}, \tilde{\mathbf{v}}_p \rangle_{\mathbf{H}_0^1(\Omega)} \le \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \|\tilde{\mathbf{v}}_p\|_{\mathbf{H}_0^1(\Omega)} \le C_{\text{div}} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \|p\|_{L^2(\Omega)}, $$

so that: $\|p\|_{L^2(\Omega)} \le C_{\text{div}} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}$. $\square$

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

## 2.2. Comments on the stability constant

Using (18) in (17), we have:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \alpha_T (C_{\max})^{-1} \|(\mathbf{u}', p')\|_{\mathcal{X}} \|(\mathbf{v}^*, q^*)\|_{\mathcal{X}} $$

Let us set $\|T\| := \sup_{((\mathbf{u}', p') \in \mathcal{X} \setminus (\mathbf{0},0))} \frac{\|T((\mathbf{u}', p'))\|_{\mathcal{X}}}{\|(\mathbf{u}', p')\|_{\mathcal{X}}}$. According to (18), we have the bound: $\|T\| \le C_{\max}$.

Hence, we have:

$$ a_S((\mathbf{u}', p'), (\mathbf{v}^*, q^*)) \ge \frac{\alpha_T}{\|T\|} \|(\mathbf{u}', p')\|_{\mathcal{X}} \|(\mathbf{v}^*, q^*)\|_{\mathcal{X}} $$

We recover the first Banach-Nečas-Babuška condition [27, Theorem 25.9, (BNB1)].

Thus, the $T$-coercivity approach gives an overestimate of the stability constant $\alpha$ given below:

$$ \frac{\alpha_T}{\|T\|} \ge \alpha := \inf_{(\mathbf{v}, q) \in \mathcal{X} \setminus (\mathbf{0},0)} \sup_{((\mathbf{u}', p') \in \mathcal{X} \setminus (\mathbf{0},0))} \frac{a_S((\mathbf{u}', p'), (\mathbf{v}, q))}{\|(\mathbf{u}', p')\|_{\mathcal{X}} \|(\mathbf{v}, q)\|_{\mathcal{X}}} $$

It suggests that the stability constant $\alpha$ is proportional to the parameter $\nu$ and depends on the constant $C_{\text{div}}$, therefore on the shape of the domain. More precisely, our estimate gives:

$$ \alpha \le \frac{\nu}{2} \times \left\{ \begin{array}{ll} \left(C_{\text{div}} \sqrt{2 + (C_{\text{div}})^2}\right)^{-1} & \text{if } 1 \le C_{\text{div}} \le \sqrt{2}, \\ \left(\sqrt{2} (C_{\text{div}})^2\right)^{-1} & \text{if } C_{\text{div}} \ge \sqrt{2}. \end{array} \right. $$
given below:

$$ \frac{\alpha_T}{\|T\|} \geq \alpha := \inf_{(\mathbf{v},q) \in \mathcal{X} \setminus (0,0)} \sup_{(\mathbf{u}',\mathbf{p}') \in \mathcal{X} \setminus (0,0)} \frac{a_S((\mathbf{u}',\mathbf{p}'), (\mathbf{v},q))}{\|(\mathbf{u}',\mathbf{p}')\|_{\mathcal{X}} \|(\mathbf{v},q)\|_{\mathcal{X}}} $$

It suggests that the stability constant $\alpha$ is proportional to the parameter $\nu$ and depends on the constant $C_{\text{div}}$, therefore on the shape of the domain. More precisely, our estimate gives:

$$ \alpha \leq \frac{\nu}{2} \times \begin{cases} (C_{\text{div}} \sqrt{2 + (C_{\text{div}})^2})^{-1} & \text{if } 1 \leq C_{\text{div}} \leq \sqrt{2}, \\ (\sqrt{2} (C_{\text{div}})^2)^{-1} & \text{if } C_{\text{div}} \geq \sqrt{2}. \end{cases} $$

In our computations, $\alpha_T$ depends on the choice of the parameters $\eta$ and $\gamma$, so that it could be further optimized to minimize $\frac{\alpha_T}{\|T\|}$. Studying the bilinear form $a_S((\mathbf{u}',\mathbf{p}'), T((\mathbf{v},q)))$, leads to an alternative variational formulation of Stokes problem, as proposed in [15]. It does not depend on the parameters $\eta$ and $\gamma$ because it is coercive. However, the new variational formulation requires a specific treatment of the right-hand side.

### 2.3. Conforming discretization and discrete well-posedness.
If we were using a conforming discretization to solve Problem (2.11) (e.g. Taylor-Hood finite elements [16]), we would use the bilinear form $a_S(\cdot, \cdot)$ to state the discrete variational formulation. Let us call the discrete spaces $X_{c,h} \subset H_0^1(\Omega)$ and $Q_{c,h} \subset L_{zmv}^2(\Omega)$. Then to prove the discrete T-coercivity, we would need to state the discrete counterpart to Proposition 1. To do so, we can build a linear operator $\Pi_c : X \to X_h$, known as Fortin operator, such that (see e.g. [17, §8.4.1]):

$$ (2.19) \quad \exists C_c \forall \mathbf{v} \in H^1(\Omega) \quad \| \mathrm{Grad} \Pi_c \mathbf{v} \|_{L^2(\Omega)} \leq C_c \| \mathrm{Grad} \mathbf{v} \|_{L^2(\Omega)}, $$

$$ (2.20) \quad \forall \mathbf{v} \in H^1(\Omega) \quad (\mathrm{div} \, \Pi_c \mathbf{v}, q_h)_{L^2(\Omega)} = (\mathrm{div} \, \mathbf{v}, q_h)_{L^2(\Omega)}, \quad \forall q_h \in Q_{c,h}. $$

Using a nonconforming discretisation, we will not use the bilinear form $a_S(\cdot, \cdot)$ to exhibit the discrete variational formulation, but we will need a similar operator to (2.19)-(2.20) to prove the discrete T-coercivity, which is stated in Theorem 4.

### 2.4. Discretization notations.
We call $(O, (x_{d'}), d_{d'}=1)$ the Cartesian coordinates system, of orthonormal basis $(e_{d'})_{d'=1}^d$. Consider $(\mathcal{T}_h)_h$ a simplicial triangulation sequence of $\Omega$, where $h$ denotes the mesh size. The triangulations are régular in the sense of Ciarlet. For a triangulation $\mathcal{T}_h$, we use the following index sets:

* $\mathcal{I}_K$ denotes the index set of the elements, such that $\mathcal{T}_h := \bigcup_{\ell \in \mathcal{I}_K} K_\ell$ is the set of elements.

* $\mathcal{I}_F$ denotes the index set of the facets$^2$, such that $\mathcal{F}_h := \bigcup_{f \in \mathcal{I}_F} F_f$ is the set of facets.

Let $\mathcal{I}_F = \mathcal{I}_F^i \cup \mathcal{I}_F^b$, where $\forall f \in \mathcal{I}_F^i, F_f \subset \Omega$ and $\forall f \in \mathcal{I}_F^b, F_f \in \partial\Omega$.

* $\mathcal{I}_S$ denotes the index set of the vertices, such that $(S_j)_{j \in \mathcal{I}_S}$ is the set of vertices.

Let $\mathcal{I}_S = \mathcal{I}_S^i \cup \mathcal{I}_S^b$, where $\forall j \in \mathcal{I}_S^i, S_j \in \Omega$ and $\forall j \in \mathcal{I}_S^b, S_j \in \partial\Omega$.

We also define the following index subsets:

<footnote>$^2$The term facet stands for face (resp. edge) when $d=3$ (resp. $d=2$).</footnote>
In our computations, $\alpha_T$ depends on the choice of the parameters $\eta$ and $\gamma$, so that it could be further optimized to minimize $\frac{\alpha_T}{||T||}$. Studying the bilinear form $a_S((\mathbf{u}', p'), T((\mathbf{v}, q)))$, leads to an alternative variational formulation of Stokes problem, as proposed in [19]. It does not depend on the parameters $\eta$ and $\gamma$ because it is coercive. However, the new variational formulation requires a specific treatment of the right-hand side.

## 2.3. Conforming discretization and discrete well-posedness

If we were using a conforming discretization to solve Problem (12) (e.g. Taylor-Hood finite elements [47]), we would use the bilinear form $a_S(\cdot, \cdot)$ to state the discrete variational formulation. Let us call the discrete spaces $X_{c,h} \subset H_0^1(\Omega)$ and $Q_{c,h} \subset L_{zmv}^2(\Omega)$. Then to prove the discrete T-coercivity, we would need to state the discrete counterpart to Proposition 3. To do so, we can build a linear operator $\Pi_c: X \to X_{c,h}$, known as Fortin operator, such that (see e.g. [8, §8.4.1]):

$$ \exists C_c \text{ such that } \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \quad \| \mathrm{Grad} \Pi_c \mathbf{v} \|_{L^2(\Omega)} \le C_c \| \mathrm{Grad} \mathbf{v} \|_{L^2(\Omega)}, \quad (20) $$

$$ \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \quad (\mathrm{div} \Pi_c \mathbf{v}, q_h)_{L^2(\Omega)} = (\mathrm{div} \mathbf{v}, q_h)_{L^2(\Omega)}, \quad \forall q_h \in Q_{c,h}. \quad (21) $$

Using a nonconforming discretization, we will not use the bilinear form $a_S(\cdot, \cdot)$ to exhibit the discrete variational formulation, but we will need a similar operator to (20)-(21) to prove the discrete T-coercivity, which is stated in Theorem 16.

## 2.4. Discretization notations

We call $(O, (x_{d'})^d_{d'=1})$ the Cartesian coordinates system, of orthonormal basis $(e_{d'})^d_{d'=1}$. Consider $(\mathcal{T}_h)_h$ a simplicial triangulation sequence of $\Omega$, where $h$ denotes the mesh size. The triangulations are regular in the sense of Ciarlet. For a triangulation $\mathcal{T}_h$, we use the following index sets:

* $\mathcal{I}_K$ denotes the index set of the elements, such that $\mathcal{T}_h := \bigcup_{\ell \in \mathcal{I}_K} K_\ell$ is the set of elements.

* $\mathcal{I}_F$ denotes the index set of the facets<fnref n="1" />, such that $\mathcal{T}_h := \bigcup_{f \in \mathcal{I}_F} F_f$ is the set of facets.

Let $\mathcal{I}_F = \mathcal{I}_F^i \cup \mathcal{I}_F^b$, where $\forall f \in \mathcal{I}_F^i, F_f \subset \Omega$ and $\forall f \in \mathcal{I}_F^b, F_f \in \partial\Omega$.

* $\mathcal{I}_S$ denotes the index set of the vertices, such that $(S_j)_{j \in \mathcal{I}_S}$ is the set of vertices.

Let $\mathcal{I}_S = \mathcal{I}_S^i \cup \mathcal{I}_S^b$, where $\forall j \in \mathcal{I}_S^i, S_j \in \Omega$ and $\forall j \in \mathcal{I}_S^b, S_j \in \partial\Omega$.

We also define the following index subsets:

* $\forall \ell \in \mathcal{I}_K, \mathcal{I}_F^{(i,b)} = \{f \in \mathcal{I}_F^{(i,b)} | F_f \in K_\ell\}, \mathcal{I}_S, \ell = \{j \in \mathcal{I}_S | S_j \in K_\ell\}$.

* $\forall j \in \mathcal{I}_S, \mathcal{I}_K, j = \{\ell \in \mathcal{I}_K | S_j \in K_\ell\}$.

For all $\ell \in \mathcal{I}_K$, we call $h_\ell$ the diameter of $K_\ell$ and $\rho_\ell$ the diameter of the sphere inscribed in $K_\ell$, and we let: $\sigma_\ell = \frac{h_\ell}{\rho_\ell}$, $h = \max_{\ell \in \mathcal{I}_K} h_\ell$. When the $(\mathcal{T}_h)_h$ is a shape-regular triangulation sequence (see e.g. [26, Definition 11.2]), there exists a constant $\sigma > 1$, called the shape regularity parameter, such that for all $h$, for all $\ell \in \mathcal{I}_K$, $\sigma_\ell \le \sigma$. For all $f \in \mathcal{I}_F$, $M_f$ denotes the barycenter of $F_f$, and by $\mathbf{n}_f$ its unit normal (outward oriented if $F_f \in \partial\Omega$). For all $j \in \mathcal{I}_S$, for all $\ell \in \mathcal{I}_K, j$, $\lambda_{j,\ell}$ denotes the barycentric coordinate of $S_j$ in $K_\ell$; $F_{j,\ell}$ denotes the face opposite to vertex $S_j$ in element $K_\ell$, and $\mathbf{x}_{j,\ell}$ denotes its barycenter. We call $\mathcal{S}_{j,\ell}$ the outward normal vector of $F_{j,\ell}$ and of norm $|\mathcal{S}_{j,\ell}| = |F_{j,\ell}|$.

<footnote n="1">¹The term facet stands for face (resp. edge) when $d=3$ (resp. $d=2$).</footnote>
* $\forall \ell \in \mathcal{I}_K, \mathcal{I}_{F,\ell}^{(i,b)} = \{f \in \mathcal{I}_F^{(i,b)} | F_f \in K_\ell\}, \mathcal{I}_{S,\ell} = \{j \in \mathcal{I}_S | S_j \in K_\ell\}.$

* $\forall j \in \mathcal{I}_S, \mathcal{I}_{K,j} = \{\ell \in \mathcal{I}_K | S_j \in K_\ell\}, N_j := \text{card}(\mathcal{I}_{K,j}).$

For all $\ell \in \mathcal{I}_K$, we call $h_\ell$ and $\rho_\ell$ the diameters of $K_\ell$ and its inscribed sphere respectively, and we let: $\sigma_\ell = \frac{h_\ell}{\rho_\ell}$, $h = \max_{\ell \in \mathcal{I}_K} h_\ell$. When the $(\mathcal{T}_h)_h$ is a shape-regular triangulation sequence (see e.g. [18, def. 11.2]), it exists a constant $\sigma > 1$, called the shape regularity parameter, such that for all $h$, for all $\ell \in \mathcal{I}_K$, $\sigma_\ell \le \sigma$. For all $f \in \mathcal{I}_F$, $M_f$ denotes the barycentre of $F_f$, and by $\mathbf{n}_f$ its unit normal (outward oriented if $F_f \in \partial\Omega$). For all $j \in \mathcal{I}_S$, for all $\ell \in \mathcal{I}_{K,j}$, $\lambda_{j,\ell}$ denotes the barycentric coordinate of $S_j$ in $K_\ell$; $F_{j,\ell}$ denotes the face opposite to vertex $S_j$ in element $K_\ell$, and $\mathbf{x}_{j,\ell}$ denotes its barycentre. We call $S_{j,\ell}$ the outward normal vector of $F_{j,\ell}$ and of norm $|S_{j,\ell}| = |F_{j,\ell}|$.

Let introduce spaces of piecewise regular elements:

We set $\mathcal{P}_h H^1 = \{v \in L^2(\Omega); \forall \ell \in \mathcal{I}_K, v|_{K_\ell} \in H^1(K_\ell)\}$, endowed with the scalar product :

$$ (v, w)_h := \sum_{\ell \in \mathcal{I}_K} (\mathbf{grad} v, \mathbf{grad} w)_{\mathbf{L}^2(K_\ell)} \quad \|v\|_h^2 = \sum_{\ell \in \mathcal{I}_K} \| \mathbf{grad} v \|_{\mathbf{L}^2(K_\ell)}^2. $$

We set $\mathcal{P}_h H^1 = [\mathcal{P}_h H^1]^d$, endowed with the scalar product :

$$ (\mathbf{v}, \mathbf{w})_h := \sum_{\ell \in \mathcal{I}_K} (\mathbf{Grad} \mathbf{v}, \mathbf{Grad} \mathbf{w})_{\mathbf{L}^2(K_\ell)} \quad \|\mathbf{v}\|_h^2 = \sum_{\ell \in \mathcal{I}_K} \| \mathbf{Grad} \mathbf{v} \|_{\mathbf{L}^2(K_\ell)}^2. $$

Let $f \in \mathcal{I}_F^i$ such that $F_f = \partial K_L \cap \partial K_R$ and $\mathbf{n}_f$ is outward $K_L$ oriented. The jump of $v \in \mathcal{P}_h H^1$ across the facet $F_f$ is defined as follows: $[v]_{F_f} := v|_{K_L} - v|_{K_R}$. For $f \in \mathcal{I}_F^b$, we set: $[v]_{F_f} := v|_{F_f}$.

We set $\mathcal{P}_h \mathbf{H}(\text{div}) = \{\mathbf{v} \in \mathbf{L}^2(\Omega); \forall \ell \in \mathcal{I}_K, \mathbf{v}|_{K_\ell} \in \mathbf{H}(\text{div}; K_\ell)\}$, and we define the operator $\text{div}_h$ such that for all $\mathbf{v} \in \mathcal{P}_h \mathbf{H}(\text{div})$, $\text{div}_h \mathbf{v} \in L^2(\Omega)$ is such that:

$$ \forall q \in L^2(\Omega), \quad (\text{div}_h \mathbf{v}, q)_{L^2(\Omega)} = \sum_{\ell \in \mathcal{I}_K} (\text{div} \mathbf{v}, q)_{L^2(K_\ell)}. $$

We recall classical finite elements estimates [18]. Let $\hat{K}$ be the reference simplex. For $\ell \in \mathcal{I}_K$, we denote by $T_\ell : \hat{K} \to K_\ell$ an affine invertible mapping such that $T_\ell(\hat{K}) = K_\ell$, $T_\ell(\partial\hat{K}) = \partial K_\ell$. We set $T_\ell(\hat{\mathbf{x}}) = \mathbb{B}_\ell \hat{\mathbf{x}} + \mathbf{b}_\ell$, where $\mathbb{B}_\ell \in \mathbb{R}^{d \times d}$ and $\mathbf{b}_\ell \in \mathbb{R}^d$. Let $J_\ell = \det(\mathbb{B}_\ell)$. There holds:

$$ (2.21) \qquad |J_\ell| = d! |K_\ell|, \quad \|\mathbb{B}_\ell\| = \frac{h_\ell}{\rho_{\hat{K}}}, \quad \|\mathbb{B}_\ell^{-1}\| = \frac{h_{\hat{K}}}{\rho_\ell}. $$

Let $f \in \mathcal{I}_{F,\ell}$. According to [2, equation (2.17)], we have:

$$ (2.22) \qquad |F_f| |K_\ell|^{-1} \le (\rho_\ell)^{-1}. $$

For $v \in L^2(K_\ell)$, we set $\hat{v}_\ell = v \circ T_\ell$. Let $v \in \mathcal{P}_h H^1$. By changing the variable, $\mathbf{grad} v|_{K_\ell} = (\mathbb{B}_\ell^{-1})^T \mathbf{grad}_{\hat{\mathbf{x}}} \hat{v}_\ell$, and it holds:

$$ (2.23) \qquad \begin{aligned} (i) \quad & \| \mathbf{grad} v \|_{\mathbf{L}^2(K_\ell)}^2 \le \| \mathbb{B}_\ell^{-1} \|^2 |K_\ell| \| \mathbf{grad}_{\hat{\mathbf{x}}} \hat{v}_\ell \|_{\mathbf{L}^2(\hat{K})}^2, \\ (ii) \quad & \| \mathbf{grad}_{\hat{\mathbf{x}}} \hat{v}_\ell \|_{\mathbf{L}^2(\hat{K})}^2 \le \| \mathbb{B}_\ell \|^2 |K_\ell|^{-1} \| \mathbf{grad} v \|_{\mathbf{L}^2(K_\ell)}^2. \end{aligned} $$

We will use the following notations :

$$ (2.24) \qquad \forall \ell \in \mathcal{I}_K, \forall v \in L^2(K_\ell), \quad \underline{v}_\ell = \int_{K_\ell} v/|K_\ell|, $$

$$ (2.25) \qquad \forall f \in \mathcal{I}_F, \forall v \in L^2(F_f), \quad \text{where } \underline{v}_f = \int_{F_f} v/|F_f|. $$
Let us introduce spaces of piecewise regular elements.

We set $\mathcal{P}_h H^1 = \{v \in L^2(\Omega) \mid \forall \ell \in \mathcal{I}_K, v|_{K_\ell} \in H^1(K_\ell)\}$, endowed with the scalar product :

$$ (v, w)_h := \sum_{\ell \in \mathcal{I}_K} (\mathbf{grad} v, \mathbf{grad} w)_{L^2(K_\ell)} \quad \|v\|_h^2 = \sum_{\ell \in \mathcal{I}_K} \|\mathbf{grad} v\|_{L^2(K_\ell)}^2. $$

We set $\mathcal{P}_h H^1 = (\mathcal{P}_h H^1)^d$, endowed with the scalar product :

$$ (\mathbf{v}, \mathbf{w})_h := \sum_{\ell \in \mathcal{I}_K} (\mathbf{Grad} \mathbf{v}, \mathbf{Grad} \mathbf{w})_{L^2(K_\ell)} \quad \|\mathbf{v}\|_h^2 = \sum_{\ell \in \mathcal{I}_K} \|\mathbf{Grad} \mathbf{v}\|_{L^2(K_\ell)}^2. $$

Let $f \in \mathcal{I}_F^i$ such that $F_f = \partial K_L \cap \partial K_R$ and $n_f$ is outward $K_L$ oriented. The jump of $v \in \mathcal{P}_h H^1$ across the facet $F_f$ is defined as follows: $[v]_{F_f} := v|_{K_L} - v|_{K_R}$. For $f \in \mathcal{I}_F^b$, we set: $[v]_{F_f} := v|_{F_f}$.

We set $\mathcal{P}_h H(\text{div}) = \{\mathbf{v} \in L^2(\Omega) \mid \forall \ell \in \mathcal{I}_K, \mathbf{v}|_{K_\ell} \in \mathbf{H}(\text{div}; K_\ell)\}$, and we define the operator $\text{div}_h$ such that for all $\mathbf{v} \in \mathcal{P}_h H(\text{div})$, $\text{div}_h \mathbf{v} \in L^2(\Omega)$ is such that:

$$ \forall q \in L^2(\Omega), \quad (\text{div}_h \mathbf{v}, q)_{L^2(\Omega)} = \sum_{\ell \in \mathcal{I}_K} (\text{div} \mathbf{v}, q)_{L^2(K_\ell)}. $$

We use the notation $A \lesssim B$ for $A \lesssim CB$ where A and B are scalar quantities and C is a generic positive constant which is independent of the sequence $(\mathcal{I}_h)_h$ and the quantities of interest.

We recall classical finite elements estimates [26]. Let $\hat{K}$ be the reference simplex. For $\ell \in \mathcal{I}_K$, we denote by $T_\ell: \hat{K} \to K_\ell$ an affine invertible mapping such that $T_\ell(\hat{K}) = K_\ell$, $T_\ell(\partial \hat{K}) = \partial K_\ell$. We set $T_\ell(\hat{\mathbf{x}}) = \mathbb{B}_\ell \hat{\mathbf{x}} + \mathbf{b}_\ell$, where $\mathbb{B}_\ell \in \mathbb{R}^{d \times d}$ and $\mathbf{b}_\ell \in \mathbb{R}^d$. Let $J_\ell = \det(\mathbb{B}_\ell)$. There holds:

$$ |J_\ell| = d! |K_\ell|, \quad \|\mathbb{B}_\ell\| = \frac{h_\ell}{\rho_{\hat{K}}}, \quad \|\mathbb{B}_\ell^{-1}\| = \frac{h_{\hat{K}}}{\rho_\ell}. \qquad (22) $$

Let $f \in \mathcal{I}_{F, \ell}$. According to [22, Equation (2.17)], we have:

$$ |F_f| |K_\ell|^{-1} \lesssim (\rho_\ell)^{-1}. \qquad (23) $$

For $v \in L^2(K_\ell)$, we set $\hat{v}_\ell = v \circ T_\ell$. Let $v \in \mathcal{P}_h H^1$. By changing the variable, $\mathbf{grad} v|_{K_\ell} = (\mathbb{B}_\ell^{-1})^T \mathbf{grad}_{\hat{\mathbf{x}}} \hat{v}_\ell$, and it holds:

$$ \| \mathbf{grad} v \|_{L^2(K_\ell)}^2 \lesssim \| \mathbb{B}_\ell^{-1} \|^2 |K_\ell| \| \mathbf{grad}_{\hat{\mathbf{x}}} \hat{v}_\ell \|_{L^2(\hat{K})}^2, \qquad (24) $$

$$ \| \mathbf{grad}_{\hat{\mathbf{x}}} \hat{v}_{\ell} \|_{L^2(\hat{K})}^2 \le \| \mathbb{B}_{\ell} \|^2 |K_{\ell}|^{-1} \| \mathbf{grad} v \|_{L^2(K_{\ell})}^2. \qquad (25) $$

We will use the following notations:

$$ \forall \ell \in \mathcal{I}_K, \forall v \in L^2(K_\ell), \quad \underline{v}_\ell = \int_{K_\ell} v/|K_\ell|, \qquad (26) $$

$$ \forall f \in \mathcal{I}_F, \forall v \in L^2(F_f), \quad \underline{v}_f = \int_{F_f} v/|F_f|. \qquad (27) $$

We recall the Poincaré-Steklov inequality in cells:

**Proposition 8 ( [26, Lemma 12.11]).** For all $\ell \in \mathcal{I}_K$ ($K_\ell$ is a convex set), for all $v \in H^1(K_\ell)$:

$$ \|v - \underline{v}_\ell\|_{L^2(K_\ell)} \le \pi^{-1} h_\ell \|\mathbf{grad} v\|_{L^2(K_\ell)}. \qquad (28) $$

We will need the following Poincaré-Steklov inequality on faces:

**Proposition 9 ( [27, Lemma 36.8]).** For all $\ell \in \mathcal{I}_K$, for all $v \in H^1(K_\ell)$ and for all $f \in \mathcal{I}_{F, \ell}$, we have:

$$ \|v - \underline{v}_f\|_{L^2(F_f)} \lesssim \left( \frac{|F_f|}{|K_\ell|} \right)^{\frac{1}{2}} h_\ell \|\mathbf{grad} v\|_{L^2(K_\ell)} \lesssim (\sigma_\ell)^{\frac{1}{2}} (h_\ell)^{\frac{1}{2}} \|\mathbf{grad} v\|_{L^2(K_\ell)}. \qquad (29) $$

**Proof.** We have: $v - \underline{v}_f = (v - \underline{v}_\ell) - |F_f|^{-1} \int_{F_f} (v - \underline{v}_\ell)$. Hence: $\|v - \underline{v}_f\|_{L^2(F_f)} \le 2\|v - \underline{v}_\ell\|_{L^2(F_f)}$. Changing the variable, using the continuity of the trace operator, we have: $\|v - \underline{v}_\ell\|_{L^2(F_f)} \lesssim |F_f|^{\frac{1}{2}} \|\hat{v}_\ell - \hat{v}_\ell\|_{H^1(\hat{K})}$. Using (28) in $\hat{K}$, we obtain $\|v - \underline{v}_\ell\|_{L^2(F_f)} \lesssim |F_f|^{\frac{1}{2}} \|\mathbf{grad}_{\hat{\mathbf{x}}} \hat{v}_\ell\|_{L^2(\hat{K})}$. Using (25) and applying (23), we get (29). $\square$
We recall the Poincaré-Steklov inequality in cells:

**Proposition 3.** For all $\ell \in \mathcal{I}_K$ ($K_\ell$ is a convex set), for all $v \in H^1(K_\ell)$:

$$ (2.26) \qquad \|v - \underline{v}_\ell\|_{L^2(K_\ell)} \le \pi^{-1} h_\ell \|\text{grad } v\|_{L^2(K_\ell)}, \quad [18, \text{Lemma 12.11}]. $$

We will need the following Poincaré-Steklov inequality on faces [8, Lemma 36.8]:

**Proposition 4.** For all $\ell \in \mathcal{I}_K$, for all $v \in H^1(K_\ell)$ and for all $f \in \mathcal{I}_{F,\ell}$, we have:

$$ \begin{align} \tag{2.27} \|v - \underline{v}_f\|_{L^2(F_f)} &\lesssim \left( \frac{|F_f|}{|K_\ell|} \right)^{\frac{1}{2}} h_\ell \|\text{grad } v\|_{L^2(K_\ell)}, \nonumber \\ &\lesssim (\sigma_\ell)^{\frac{1}{2}} (h_\ell)^{\frac{1}{2}} \|\text{grad } v\|_{L^2(K_\ell)}. \end{align} $$

*Proof.* We have: $v - \underline{v}_f = (v - \underline{v}_\ell) - |F|^{-1} \int_{F_f} (v - \underline{v}_\ell)$. Hence: $\|v - \underline{v}_f\|_{L^2(F_f)} \le 2\|v - \underline{v}_\ell\|_{L^2(F_f)}$. Changing the variable, using the continuity of the trace operator, we have: $\|v - \underline{v}_\ell\|_{L^2(F_f)} \lesssim |F_f|^{\frac{1}{2}} \|\hat{v}_\ell - \hat{v}_\ell\|_{H^1(\hat{K})}$. Using (2.26) in $\hat{K}$, changing the variable again, we get (2.27), applying (2.22). $\square$

For all $D \subset \mathbb{R}^d$, and $k \in \mathbb{N}$, we call $P^k(D)$ the set of order $k$ polynomials on $D$, $\mathbf{P}^k(D) = (P^k(D))^d$, and we consider the broken polynomial space:

$$ P_{\text{disc}}^{k}(\mathcal{T}_{h}) = \{q \in L^{2}(\Omega); \quad \forall \ell \in \mathcal{I}_{K}, q|_{K_{\ell}} \in P^{k}(K_{\ell})\}, \quad \mathbf{P}_{\text{disc}}^{k}(\mathcal{T}_{h}) := (P_{\text{disc}}^{k}(\mathcal{T}_{h}))^{d}. $$

**2.5. Non-conforming discretization and discrete well-posedness.** The non-conforming finite element method was introduced by Crouzeix and Raviart in [2] to solve Stokes problem (2.3). We approximate the vector space $\mathbf{H}^1(\Omega)$ component by component by piecewise polynomials of order $k \in \mathbb{N}^*$.

Let us consider $X_h$ (resp. $X_{0,h}$), the space of nonconforming approximation of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order $k$:

$$ \begin{align} \tag{2.28} X_h &= \left\{ v_h \in P_{\text{disc}}^k(\mathcal{T}_h) ; \forall f \in \mathcal{I}_F^i, \forall q_h \in P^{k-1}(F_f), \int_{F_f} [v_h]_{F_f} q_h = 0 \right\}; \nonumber \\ X_{0,h} &= \left\{ v_h \in X_h ; \forall f \in \mathcal{I}_F^b, \forall q_h \in P^{k-1}(F_f), \int_{F_f} v_h q_h = 0 \right\}. \end{align} $$

The condition on the jumps of $v_h$ on the inner facets is often called the patch-test condition. It allows to prove a discrete Poincaré-Steklov inequality, using the following Lemma:

**Lemma 1.** For all $v \in L^2(\Omega)$, it exists $s \in \mathbf{H}^1(\Omega)$ such that:

$$ (2.29) \qquad \operatorname{div} s = v \quad \text{and} \quad \|s\|_{\mathbf{L}^2(\Omega)} + h_{\Omega} \| \operatorname{Grad} s \|_{\mathbf{L}^2(\Omega)} \le C_{\Omega} h_{\Omega} \|v\|_{\mathbf{L}^2(\Omega)}, $$
where the dimensionless constant $C_{\Omega}$ depends on $C_{\text{div}}$ and $C_{\text{PS}}$.

*Proof.* Let $\underline{v} = \int_{\Omega} v/|\Omega|$, and $s_0 \in \mathbf{H}_0^1(\Omega)$ be such that $\operatorname{div} s_0 = v - \underline{v}$ and $\| \operatorname{Grad} s_0 \|_{\mathbf{L}^2(\Omega)} \le C_{\text{div}} \|v\|_{\mathbf{L}^2(\Omega)}$ (cf. Prop. 1). Let $d' \in \{1, \dots, d\}$. We consider $s := s_0 + \underline{v}(x_{d'} - \underline{x}_{d'}) e_{d'} \in \mathbf{H}^1(\Omega)$, where $\underline{x}_{d'} = \int_{\Omega} x_{d'}/|\Omega|$. We have: $\operatorname{div} s = v$ and $\| \operatorname{Grad} s \|_{\mathbf{L}^2(\Omega)} \le \tilde{C}_{\text{div}} \|v\|_{\mathbf{L}^2(\Omega)}$, where $(\tilde{C}_{\text{div}})^2 = (C_{\text{div}})^2 + 1$. Using inequality (2.4), one can prove that $\|s\|_{\mathbf{L}^2(\Omega)} \le h_{\Omega} \tilde{C}_{\Omega} \|v\|_{\mathbf{L}^2(\Omega)}$, where $(\tilde{C}_{\Omega})^2 = 2 ((C_{\text{PS}} C_{\text{div}})^2 + 1)$. Setting $C_{\Omega} = \tilde{C}_{\Omega} + \tilde{C}_{\text{div}}$, we obtain (2.29). $\square$
For all $D \subset \mathbb{R}^d$, and $k \in \mathbb{N}$, we call $P^k(D)$ the set of order $k$ polynomials on $D$, $\mathbf{P}^k(D) = (P^k(D))^d$, and we consider the broken polynomial space:

$$P_{\mathrm{disc}}^{k}(\mathcal{T}_{h}) = \{q \in L^{2}(\Omega); \quad \forall \ell \in \mathcal{I}_{K}, q|_{K_{\ell}} \in P^{k}(K_{\ell})\}, \quad \mathbf{P}_{\mathrm{disc}}^{k}(\mathcal{T}_{h}) := (P_{\mathrm{disc}}^{k}(\mathcal{T}_{h}))^{d}.$$

## 2.5. Nonconforming discretization and discrete well-posedness

The nonconforming finite element method was introduced by Crouzeix and Raviart in [22] to solve Stokes problem (3). We approximate the vector space $\mathbf{H}^1(\Omega)$ component by component by piecewise polynomials of order $k \in \mathbb{N}^*$. Let us consider $X_h$ (resp. $X_{0,h}$), the space of nonconforming approximation of $\mathbf{H}^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order $k$:

$$X_h = \left\{ v_h \in P_{\mathrm{disc}}^k(\mathcal{T}_h) \mid \forall f \in \mathcal{F}_F^i, \forall q_h \in P^{k-1}(F_f), \int_{F_f} [v_h]_{F_f} q_h = 0 \right\}; \\ X_{0,h} = \left\{ v_h \in X_h \mid \forall f \in \mathcal{F}_F^b, \forall q_h \in P^{k-1}(F_f), \int_{F_f} v_h q_h = 0 \right\}. \tag{30}$$

The condition on the jumps of $v_h$ on the inner facets is often called the patch-test condition. It allows to prove a discrete Poincaré-Steklov inequality, using Corollary 4, and the proof of [40, Theorem D.1]. The proof of [27, Lemma 36.6] is similar, but the vector $\mathbf{s}$ defined in Corollary 4 is constructed in [27, Lemma 36.6] as the gradient of a scalar function, so that it gives a lower estimate when $\Omega$ is nonconvex. Alternative proofs are given in [11, 48].

**Proposition 10.** *The following discrete Poincaré-Steklov inequality holds:*

$$\forall v_h \in X_{0,h}, \quad \|v_h\|_{L^2(\Omega)} \lesssim \sigma C_{\Omega} h_{\Omega} \|v_h\|_{h}. \tag{31}$$

**Proof.** Let $v_h \in X_{0,h}$. According to Corollary 4, there exists $\mathbf{s} \in \mathbf{H}^1(\Omega)$ such that:

$$\mathrm{div} \, \mathbf{s} = v_h \quad \text{and} \quad \|\mathbf{s}\|_{L^2(\Omega)} + h_{\Omega} \|\mathrm{Grad} \, \mathbf{s}\|_{L^2(\Omega)} \le C_{\Omega} h_{\Omega} \|v_h\|_{L^2(\Omega)}. \tag{32}$$

We have, by integration by parts:

$$\|v_h\|_{L^2(\Omega)}^2 = (v_h, \mathrm{div}\,\mathbf{s})_{L^2(\Omega)} = - \sum_{\ell \in \mathcal{I}_K} (\mathrm{grad}\,v_h, \mathbf{s})_{L^2(K_\ell)} + \sum_{\ell \in \mathcal{I}_K} \sum_{f \in \mathcal{F}_{F,\ell}} (v_h, \mathbf{s} \cdot \mathbf{n}_{f,\ell})_{L^2(F_f)}. \tag{33}$$

The first term can be bounded as follows:

$$(\mathrm{grad} \, v_h, \mathbf{s})_{L^2(K_\ell)} \le \|\mathrm{grad} \, v_h\|_{L^2(K_\ell)} \|\mathbf{s}\|_{L^2(K_\ell)}. \tag{34}$$

Due to the patch-test, the second term reads:

$$\sum_{\ell \in \mathcal{I}_K} \sum_{f \in \mathcal{F}_{F,\ell}} (v_h, \mathbf{s} \cdot \mathbf{n}_{f,\ell})_{L^2(F_f)} = \sum_{\ell \in \mathcal{I}_K} \sum_{f \in \mathcal{F}_{F,\ell}} \left( v_h - \underline{v}_{h,f}, (\mathbf{s} - \underline{\mathbf{s}}_f) \cdot \mathbf{n}_{f,\ell} \right)_{L^2(F_f)}, \\ \le \sum_{\ell \in \mathcal{I}_K} \sum_{f \in \mathcal{F}_{F,\ell}} \left\| v_h - \underline{v}_{h,f} \right\|_{L^2(F_f)} \left\| (\mathbf{s} - \underline{\mathbf{s}}_f) \cdot \mathbf{n}_{f,\ell} \right\|_{L^2(F_f)}. \tag{35}$$

Using inequality (29), we have:

$$\|v_h - \underline{v}_{h,f}\|_{L^2(F_f)} \left\| (\mathbf{s} - \underline{\mathbf{s}}_f) \cdot \mathbf{n}_{f,\ell} \right\|_{L^2(F_f)} \lesssim \sigma_{\ell} h_{\ell} \|\mathrm{grad} \, v_h\|_{L^2(K_{\ell})} \|\mathrm{Grads}\|_{L^2(K_{\ell})}. \tag{36}$$

Using (36) in (35), combining the result with (34), inequality (33) now reads:

$$\|v_h\|_{L^2(\Omega)}^2 \lesssim \sum_{\ell \in \mathcal{I}_K} \|\mathrm{grad} \, v_h\|_{L^2(K_\ell)} (\|\mathbf{s}\|_{L^2(K_\ell)} + \sigma_\ell h_\ell \|\mathrm{Grads}\|_{L^2(K_\ell)}) \\ \lesssim \sigma \sum_{\ell \in \mathcal{I}_K} \|\mathrm{grad} \, v_h\|_{L^2(K_\ell)} (\|\mathbf{s}\|_{L^2(K_\ell)} + h_\Omega \|\mathrm{Grads}\|_{L^2(K_\ell)}).$$

We obtain (31) using the discrete Cauchy-Schwarz inequality and (32). □
To prove the next discrete Poincaré-Steklov inequality, we follow the proof of [19, Theorem D.1]. The proof of [8, Lemma 36.6] is similar, but the vector s defined in Lemma 1 is constructed in [8, Lemma 36.6] as the gradient of a scalar function, so that it gives a lower estimate when Ω is non-convex. Alternative proofs are given in [20, 21].

**Proposition 5.** The following discrete Poincaré-Steklov inequality holds:

$$ (2.30) \qquad \forall v_h \in X_{0,h}, \quad \|v_h\|_{L^2(\Omega)} \lesssim \sigma C_{\Omega} h_{\Omega} \|v_h\|_h. $$

*Proof.* Let $v_h \in X_{0,h}$. According to Lemma 1, it exists $s \in H^1(\Omega)$ such that:

$$ (2.31) \qquad \operatorname{div} s = v_h \quad \text{and} \quad \|s\|_{L^2(\Omega)} + h_{\Omega} \|\operatorname{Grad} s\|_{L^2(\Omega)} \le C_{\Omega} h_{\Omega} \|v_h\|_{L^2(\Omega)} $$

We have, by integration by parts:

$$ (2.32) \qquad \begin{aligned} \|v_h\|_{L^2(\Omega)}^2 &= (v_h, \operatorname{div} s)_{L^2(\Omega)}, \\ &= - \sum_{\ell \in \mathcal{I}_K} (\operatorname{grad} v_h, s)_{L^2(\mathcal{K}_\ell)} + \sum_{\ell \in \mathcal{I}_K} \sum_{f \in \mathcal{I}_{F,\ell}} (v_h, s \cdot n_{f,\ell})_{L^2(F_f)}. \end{aligned} $$

The first term can be bounded as follows:

$$ (2.33) \qquad (\operatorname{grad} v_h, s)_{L^2(\mathcal{K}_\ell)} \le \|\operatorname{grad} v_h\|_{L^2(\mathcal{K}_\ell)} \|s\|_{L^2(\mathcal{K}_\ell)}. $$

Due to the patch-test, the second term reads:

$$ (2.34) \qquad \begin{aligned} \sum_{\ell \in \mathcal{I}_K} \sum_{f \in \mathcal{I}_{F,\ell}} (v_h, s \cdot n_{f,\ell})_{L^2(F_f)} &= \sum_{\ell \in \mathcal{I}_K} \sum_{f \in \mathcal{I}_{F,\ell}} (v_h - \underline{v}_{h,f}, (s - \underline{s}_f) \cdot n_{f,\ell})_{L^2(F_f)}, \\ &\le \sum_{\ell \in \mathcal{I}_K} \sum_{f \in \mathcal{I}_{F,\ell}} \|v_h - \underline{v}_{h,f}\|_{L^2(F_f)} \|(s - \underline{s}_f) \cdot n_{f,\ell}\|_{L^2(F_f)}. \end{aligned} $$

Using inequality (2.27), we have:

$$ (2.35) \qquad \begin{aligned} & \|v_h - \underline{v}_{h,f}\|_{L^2(F_f)} \|(s - \underline{s}_f) \cdot n_{f,\ell}\|_{L^2(F_f)} \\ & \lesssim \sigma_{\ell} h_{\ell} \|\operatorname{grad} v_h\|_{L^2(\mathcal{K}_\ell)} \|\operatorname{Grad} s\|_{L^2(\mathcal{K}_\ell)}. \end{aligned} $$

Using (2.35) in (2.34), combining the result with (2.33), inequality (2.32) now reads:

$$ \begin{aligned} \|v_h\|_{L^2(\Omega)}^2 &\lesssim \sum_{\ell \in \mathcal{I}_K} \|\operatorname{grad} v_h\|_{L^2(\mathcal{K}_\ell)} \left( \|s\|_{L^2(\mathcal{K}_\ell)} + \sigma_\ell h_\ell \|\operatorname{Grad} s\|_{L^2(\mathcal{K}_\ell)} \right) \\ &\lesssim \sigma \sum_{\ell \in \mathcal{I}_K} \|\operatorname{grad} v_h\|_{L^2(\mathcal{K}_\ell)} \left( \|s\|_{L^2(\mathcal{K}_\ell)} + h_\Omega \|\operatorname{Grad} s\|_{L^2(\mathcal{K}_\ell)} \right). \end{aligned} $$

We obtain (2.30) using the discrete Cauchy-Schwarz inequality and (2.31). □

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

As a consequence of Proposition 5, we have the

**Proposition 6.** The broken norm $v_h \to \|v_h\|_h$ is a norm over $X_{0,h}$.

The space of nonconforming approximation of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order $k$ is $X_h = (X_h)^d$ (resp. $X_{0,h} = (X_{0,h})^d$). We set $\mathcal{X}_h := X_{0,h} \times Q_h$ where $Q_h = P_{disc}^{k-1}(\mathcal{T}_h) \cap L_{zmv}^2(\Omega)$. We deduce from Proposition 6 the

**Proposition 7.** The broken norm defined below is a norm on $\mathcal{X}_h$:

$$ (2.36) \qquad \| (\cdot, \cdot) \|_{\mathcal{X}_h} : \begin{cases} \mathcal{X}_h \mapsto \mathbb{R} \\ (v_h, q_h) \mapsto \left( \|v_h\|_h^2 + \nu^{-2} \|q_h\|_{L^2(\Omega)}^2 \right)^{\frac{1}{2}}. \end{cases} $$
As a consequence of Proposition 10, we have the following result.

**Proposition 11.** The broken norm $v_h \to \|v_h\|_h$ is a norm over $X_{0,h}$.

The space of nonconforming approximation of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order $k$ is $X_h = (X_h)^d$ (resp. $X_{0,h} = (X_{0,h})^d$). We set $\mathcal{X}_h := X_{0,h} \times Q_h$ where $Q_h = P_{disc}^{k-1}(\mathcal{T}_h) \cap L_{zmv}^2(\Omega)$. We deduce from Proposition 11 the following result.

**Proposition 12.** The broken norm defined below is a norm on $\mathcal{X}_h$:

$$ \| (\cdot, \cdot) \|_{\mathcal{X}_h} : \left\{ \begin{array}{ll} \mathcal{X}_h \mapsto \mathbb{R} & (\mathbf{v}_h, q_h) \to (\|\mathbf{v}_h\|_h^2 + v^{-2} \|q_h\|_{L^2(\Omega)}^2)^{\frac{1}{2}}. \\ & (\mathbf{u}'_h, p'_h) \to (\|\mathbf{u}'_h\|_h^2 + v^{-2} \|p'_h\|_{L^2(\Omega)}^2)^{\frac{1}{2}}. \end{array} \right. $$

Thus, the product space $\mathcal{X}_h$ endowed with the broken norm $\|(\cdot, \cdot)\|_{\mathcal{X}_h}$ is a Hilbert space. We consider the discrete continuous bilinear form $a_{S,h}(\cdot, \cdot)$ such that:

$$ \left\{ \begin{array}{ll} a_{S,h} : \mathcal{X}_h \times \mathcal{X}_h \to \mathbb{R} & (\mathbf{u}'_h, p'_h) \times (\mathbf{v}_h, q_h) \to v(\mathbf{u}'_h, \mathbf{v}_h)_h - (\mathrm{div}_h \mathbf{v}_h, p'_h)_{L^2(\Omega)} - (\mathrm{div}_h \mathbf{u}'_h, q_h)_{L^2(\Omega)}. \\ & (\mathbf{v}_h, q_h) \to v(\mathbf{v}_h)_{F_f} \cdot \mathbf{n}_f = 0. \end{array} \right. $$

Let us set $\mathbf{V}_h$ the discrete space of discrete divergence-free velocities:

$$ \mathbf{V}_h := \{ \mathbf{v}_h \in \mathbf{X}_{0,h} \mid \forall q_h \in Q_h, (\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = 0 \}. \quad (37) $$

We recall that the velocities in $\mathbf{V}_h$ are piecewise divergence-free:

**Proposition 13 ( [13, Lemma 3.1]).** For all $\mathbf{v}_h \in \mathbf{V}_h$, for all $\ell \in \mathcal{I}_K$, $\mathrm{div}\mathbf{v}_h|_{K_\ell} = 0$.

**Proof.** Let $\mathbf{v}_h \in \mathbf{V}_h$. Integrating by parts and using the patch-test, we have:

$$ (\mathrm{div}_h \mathbf{v}_h, 1)_{L^2(\Omega)} = \sum_{\ell \in \mathcal{I}_K} \int_{K_\ell} \mathrm{div} \mathbf{v}_h = \sum_{f \in \mathcal{I}_F} \int_{F_f} [\mathbf{v}_h]_{F_f} \cdot \mathbf{n}_f = 0. $$

Let $q_h \in P_{disc}^{k-1}(\mathcal{T}_h)$ and $q_h = \int_{\Omega} q_h |\Omega|$. Then $q_h - q_h|_{h} \in Q_h$ so that: $(\mathrm{div}_h \mathbf{v}_h, q_h - q_h)_{L^2(\Omega)} = 0$. Hence, we have: $(\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = (\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = 0$. Let $\ell \in \mathcal{I}_K$. Let $q_h \in P_{disc}^{k-1}(\mathcal{T}_h)$ such that $q_h|_{K_\ell} = \mathrm{div}_h \mathbf{v}_h|_{K_\ell}$ and for all $\ell' \in \mathcal{I}_K$, $\ell' \neq \ell$, $q_h|_{K_\ell'} = 0$. We have: $(\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = 0$ and $(\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = (\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(K_\ell)} = \|\mathrm{div}_h \mathbf{v}_h\|_{L^2(K_\ell)}^2$. Hence $\|\mathrm{div}_h \mathbf{v}_h\|_{L^2(K_\ell)}^2 = 0$. $\square$

Let $\mathcal{I}_h : X_{0,h} \to Y_{0,h}$, with $Y_{0,h} = \{v_h \in H_0^1(\Omega) \mid \forall \ell \in \mathcal{I}_K, v_h|_{K_\ell} \in P^k(K_\ell)\}$ be the averaging operator described in [26, §22.4.1]. There exists a constant $C_{\mathcal{I}_h}^{nc} \approx \sigma$ and independent of $h$ such that:

$$ \forall v_h \in X_{0,h}, \quad \|\mathcal{I}_h v_h\|_{H_0^1(\Omega)} \le C_{\mathcal{I}_h}^{nc} \|v_h\|_h. \quad (38) $$

Let $\ell_f \in \mathcal{L}(\mathcal{X}_h, \mathbb{R})$ be such that for all $(\mathbf{v}_h, q_h) \in \mathcal{X}_h$:

$$ \ell_f((\mathbf{v}_h, q_h)) = \begin{cases} (\mathbf{f}, \mathbf{v}_h)_{L^2(\Omega)} & \text{if } \mathbf{f} \in L^2(\Omega), \\ \langle \mathbf{f}, \mathcal{I}_h(\mathbf{v}_h) \rangle_{H_0^1(\Omega)} & \text{if } \mathbf{f} \in H^{-1}(\Omega). \end{cases} $$

The nonconforming discretization of Problem (12) reads:

$$ \text{Find } (\mathbf{u}_h, p_h) \in \mathcal{X}_h \text{ such that } a_{S,h}((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) = \ell_f((\mathbf{v}_h, q_h)) \quad \forall (\mathbf{v}_h, q_h) \in \mathcal{X}_h. \quad (39) $$

To prove that Problem (39) is well-posed, we will also use the T-coercivity theory. We do not need the well-posedness of the continuous problem, i.e. Proposition 6, but we will follow its proof, using a Fortin operator. This operator will be explained later, using the discrete basis functions. We will see that the discrete stability constant depends on this operator (hence polynomial order $k$).
Thus, the product space $\mathcal{X}_h$ endowed with the broken norm $\| \cdot \|_{\mathcal{X}_h}$ is a Hilbert space. We consider the discrete continuous bilinear form $a_{S,h}(\cdot, \cdot)$ such that :

$$ \left\{ \begin{array}{rcl} a_{S,h} : \mathcal{X}_h & \to & \mathbb{R} \\ (\mathbf{u}'_h, p'_h) \times (\mathbf{v}_h, q_h) & \mapsto & \nu(\mathbf{u}'_h, \mathbf{v}_h)_h - (\mathrm{div}_h \mathbf{v}_h, p'_h)_{L^2(\Omega)} - (\mathrm{div}_h \mathbf{u}'_h, q_h)_{L^2(\Omega)} \end{array} \right. . $$

Let us set $\mathbf{V}_h$ the discrete space of discrete divergence-free velocities :

$$ (2.37) \qquad \mathbf{V}_h := \{ \mathbf{v}_h \in \mathbf{X}_{0,h} \mid \forall q_h \in Q_h, (\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = 0 \}. $$

We recall that [22, Lemma 3.1]:

**Proposition 8.** For all $\mathbf{v}_h \in \mathbf{V}_h$, for all $\ell \in \mathcal{I}_K$, $\mathrm{div} \, \mathbf{v}_{h|K_\ell} = 0$.

*Proof.* Let $\mathbf{v}_h \in \mathbf{V}_h$. Integrating by parts and using the patch-test, we have:

$$ (\mathrm{div}_h \mathbf{v}_h, 1)_{L^2(\Omega)} = \sum_{\ell \in \mathcal{I}_K} \int_{K_\ell} \mathrm{div} \, \mathbf{v}_h = \sum_{f \in \mathcal{I}_F} \int_{F_f} [\mathbf{v}_h]_{F_f} \cdot \mathbf{n}_f = 0. $$

Let $q_h \in P_{disc}^{k-1}(\mathcal{T}_h)$ and $q_h = \int_{\Omega} q_h/|\Omega|$. Then $q_h - q_h \in Q_h$ so that: $(\mathrm{div}_h \mathbf{v}_h, q_h - q_h)_{L^2(\Omega)} = 0$. Hence, we have: $(\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = (\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = 0$. Let $\ell \in \mathcal{I}_K$. Let $q_h \in P_{disc}^{k-1}(\mathcal{T}_h)$ such that $q_h|_{K_\ell} = \mathrm{div} \, \mathbf{v}_{h|K_\ell}$ and for all $\ell' \in \mathcal{I}_K$, $\ell' \neq \ell$, $q_h|_{K_{\ell'}} = 0$. We have: $(\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = 0$ and $(\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)} = (\mathrm{div} \, \mathbf{v}_h, q_h)_{L^2(K_\ell)} = \| \mathrm{div} \, \mathbf{v}_h \|_{L^2(K_\ell)}^2$. Hence $\| \mathrm{div} \, \mathbf{v}_h \|_{L^2(K_\ell)}^2 = 0$. $\square$

Let $\mathcal{I}_h: X_{0,h} \to Y_{0,h}$, with $Y_{0,h} = \{v_h \in \mathbf{H}_0^1(\Omega); \forall \ell \in \mathcal{I}_K, v_{h|K_\ell} \in P^k(K_\ell)\}$ be the averaging operator described in [18, §22.4.1]. It exists a constant $C_{\mathcal{I}_h}^{nc} > 0$ independent of $\mathcal{T}_h$ such that :

$$ (2.38) \qquad \forall v_h \in X_{0,h}, \quad \|\mathcal{I}_h v_h\|_{\mathbf{H}_0^1(\Omega)} \le C_{\mathcal{I}_h}^{nc} \|v_h\|_h. $$

Let $\ell_f \in \mathcal{L}(\mathcal{X}_h, \mathbb{R})$ be such that for all $(\mathbf{v}_h, q_h) \in \mathcal{X}_h$:

If $\mathbf{f} \in \mathbf{L}^2(\Omega): \quad \ell_f((\mathbf{v}_h, q_h)) = (\mathbf{f}, \mathbf{v}_h)_{\mathbf{L}^2(\Omega)}$.

If $\mathbf{f} \in \mathbf{H}^{-1}(\Omega): \quad \ell_f((\mathbf{v}_h, q_h)) = \langle \mathbf{f}, \mathcal{I}_h(\mathbf{v}_h) \rangle_{\mathbf{H}_0^1(\Omega)}$.

The nonconforming discretization of Problem (2.11) reads:

Find $(\mathbf{u}_h, p_h) \in \mathcal{X}_h$ such that

$$ (2.39) \qquad a_{S,h}((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) = \ell_f((\mathbf{v}_h, q_h)) \quad \forall (\mathbf{v}_h, q_h) \in \mathcal{X}_h. $$

To prove that Problem (2.39) is well-posed, we will also use the T-coercivity theory. We do not need the well-posedness of the continuous problem, i.e. Prop. 2, but we will follow its proof, using a Fortin operator. This operator will be explained latter, using the discrete basis functions. We will see that the discrete stability constant depends on this operator (hence polynomial order $k$).

**Proposition 9.** Suppose that it exists a Fortin operator $\Pi_{nc}: \mathbf{H}^1(\Omega) \to \mathbf{X}_h$ such that

$$ (2.40) \qquad \exists C_{nc} \mid \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \quad \|\Pi_{nc}\mathbf{v}\|_h \le C_{nc} \|\mathrm{Grad} \, \mathbf{v}\|_{\mathbf{L}^2(\Omega)}, $$

$$ (2.41) \qquad \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \quad (\mathrm{div}_h \Pi_{nc} \mathbf{v}, q_h)_{L^2(\Omega)} = (\mathrm{div} \, \mathbf{v}, q_h)_{L^2(\Omega)}, \quad \forall q \in Q_h, $$

where the constant $C_{nc}$ does not depend on $h$. Then, the bilinear form $a_{S,h}(\cdot, \cdot)$ is T-coercive:

$$ (2.42) \qquad \exists T_h \in \mathcal{L}(\mathcal{X}_h), \text{ bijective, } \exists \alpha_{T_h} > 0, \forall (\mathbf{u}'_h, p'_h) \in \mathcal{X}_h, \\ a_{S,h}((\mathbf{u}'_h, p'_h), T_h((\mathbf{u}'_h, p'_h))) \ge \alpha_{T_h} \|(\mathbf{u}'_h, p'_h)\|_{\mathcal{X}_h}^2. $$
**Proposition 14.** Suppose that there exists a Fortin operator $\Pi_{nc} : H^1(\Omega) \to X_h$ such that

$$ \exists C_{nc} \text{ such that } \forall \mathbf{v} \in H^1(\Omega) \quad \|\Pi_{nc}\mathbf{v}\|_{h} \le C_{nc} \|\mathbf{Grad}\mathbf{v}\|_{L^2(\Omega)}, \quad (40) $$

$$ \forall \mathbf{v} \in H^1(\Omega) \quad (\mathrm{div}_h \Pi_{nc} \mathbf{v}, q_h)_{L^2(\Omega)} = (\mathrm{div} \mathbf{v}, q_h)_{L^2(\Omega)}, \quad \forall q \in Q_h, \quad (41) $$

where the constant $C_{nc}$ does not depend on $h$. Then, the bilinear form $a_{S,h}(\cdot, \cdot)$ is T-coercive:

$$ \exists T_h \in \mathcal{L}(\mathcal{X}_h), \text{ bijective, } \exists \alpha_{T_h} > 0, \forall (\mathbf{u}_h', \mathbf{p}_h') \in \mathcal{X}_h, \quad (42) $$

$$ a_{S,h}((\mathbf{u}_h', \mathbf{p}_h'), T_h((\mathbf{u}_h', \mathbf{p}_h'))) \ge \alpha_{T_h} \| (\mathbf{u}_h', \mathbf{p}_h') \|_{\mathcal{X}_h}^2. $$

We will exhibit $C_{nc}$ for $k=1$, $d=2$, $3$ in Section 3.1, obtaining $C_{nc} = 1$; and then for $k=2$, $d=2$ in Section 3.3, obtaining $C_{nc} = \sigma^2 + 1$. Let us prove Proposition 14.

**Proof.** We follow the proof of Proposition 6. Let us consider $(\mathbf{u}_h', \mathbf{p}_h') \in \mathcal{X}_h$ and let us build $(\mathbf{v}_h^*, q_h^*) \in \mathcal{X}_h$ satisfying (2) (with $V = \mathcal{X}_h$). We need three main steps.

**Step 1.** According to Proposition 3, there exists $\tilde{\mathbf{v}}_{p_h'} \in \mathbf{V}^{\perp}$ such that $\mathrm{div} \tilde{\mathbf{v}}_{p_h'} = \mathbf{p}_h'$ in $\Omega$ and $\|\tilde{\mathbf{v}}_{p_h'}\|_{H_0^1(\Omega)} \le C_{\mathrm{div}} \|\mathbf{p}_h'\|_{L^2(\Omega)}$. Let us set $\mathbf{v}_{p_h'} = v^{-1} \tilde{\mathbf{v}}_{p_h'}$. Consider $\mathbf{v}_{h,p_h'} = \Pi_{nc} \mathbf{v}_{p_h'}$, for all $q_h \in Q_h$ we have:

$$ (\mathrm{div}_h \mathbf{v}_{h,p_h'}, q_h)_{L^2(\Omega)} = v^{-1} (\mathbf{p}_h', q_h)_{L^2(\Omega)} \quad \text{and} \quad \|\mathbf{v}_{h,p_h'}\|_{h} \le v^{-1} C_{\mathrm{div}}^{nc} \|\mathbf{p}_h'\|_{L^2(\Omega)}, \quad (43) $$

where $C_{\mathrm{div}}^{nc} = C_{nc} C_{\mathrm{div}}$. Let us set $(\mathbf{v}_h^*, q_h^*) := (\gamma_{nc} \mathbf{u}_h' - \mathbf{v}_{h,p_h'}, -\gamma_{nc} \mathbf{p}_h')$, with $\gamma_{nc} > 0$. We obtain:

$$ a_{S,h}((\mathbf{u}_h', \mathbf{p}_h'), (\mathbf{v}_h^*, q_h^*)) = v \gamma_{nc} \|\mathbf{u}_h'\|_{h}^{2} + v^{-1} \|\mathbf{p}_h'\|_{L^2(\Omega)}^{2} - v(\mathbf{u}_h', \mathbf{v}_{h,p_h'})_h. \quad (44) $$

**Step 2.** In order to bound the last term of (44), we use Young's inequality and then inequality (43) so that for all $\eta_{nc} > 0$:

$$ (\mathbf{u}_h', \mathbf{v}_{h,p_h'})_h \le \frac{\eta_{nc}}{2} \|\mathbf{u}_h'\|_{h}^{2} + \frac{\eta_{nc}^{-1}}{2} \left( \frac{C_{\mathrm{div}}^{nc}}{v} \right)^{2} \|\mathbf{p}_h'\|_{L^2(\Omega)}^{2}. \quad (45) $$

**Step 3.** Using the bound (45) in (44) and choosing $\eta_{nc} = \gamma_{nc}$, we get:

$$ a_{S,h}((\mathbf{u}_h', \mathbf{p}_h'), (\mathbf{v}_h^*, q_h^*)) \ge v \left( \frac{\gamma_{nc}}{2} \|\mathbf{u}_h'\|_{h}^{2} + v^{-2} \left( 1 - \frac{(\gamma_{nc})^{-1}}{2} (C_{\mathrm{div}}^{nc})^{2} \right) \|\mathbf{p}_h'\|_{L^2(\Omega)}^{2} \right). $$

Consider now $\gamma_{nc} = (C_{\mathrm{div}}^{nc})^2$. We obtain:

$$ a_{S,h}((\mathbf{u}_h', \mathbf{p}_h'), (\mathbf{v}_h^*, q_h^*)) \ge \alpha_T C_{\mathrm{min}}^{nc} \| (\mathbf{u}_h', \mathbf{p}_h') \|_{\mathcal{X}_h}^2, \text{ where } C_{\mathrm{min}}^{nc} = \min((C_{\mathrm{div}}^{nc})^2, 1). $$

We obtain (42) with $\alpha_{T_h} = \alpha_T C_{\mathrm{min}}^{nc}$. Suppose that $C_{nc} \ge 1$. Then $\alpha_{T_h} = \alpha_T = \frac{\gamma}{2}$.

The operator $T_h$ such that $T_h(\mathbf{u}_h', \mathbf{p}_h') = (\mathbf{v}_h^*, \mathbf{p}_h^*)$ is linear and continuous. We have indeed:

$$ \| T_h(\mathbf{u}_h', \mathbf{p}_h') \|_{\mathcal{X}_h}^2 = \| \mathbf{v}_h^* \|_{h}^{2} + v^{-2} \| q_h^* \|_{L^2(\Omega)}^{2} \le (C_{\mathrm{max}}^{nc})^{2} \| (\mathbf{u}_h', \mathbf{p}_h') \|_{\mathcal{X}_h}^2 $$

$$ \text{where } C_{\mathrm{max}}^{nc} = C_{\mathrm{div}}^{nc} \left( \max \left( 2 + (C_{\mathrm{div}}^{nc})^2, 2(C_{\mathrm{div}}^{nc})^2 \right) \right)^{\frac{1}{2}}. $$

Remark that, given $(\mathbf{v}_h^*, q_h^*) \in \mathcal{X}_h$, choosing $(\mathbf{u}_h', \mathbf{p}_h') = (\gamma_{nc}^{-1} \mathbf{v}_h^* - \gamma_{nc}^{-2} \mathbf{v}_{h,q_h^*}, -\gamma_{nc}^{-1} q_h^*)$ yields $T_h((\mathbf{u}_h', \mathbf{p}_h')) = (\mathbf{v}_h^*, q_h^*)$. Hence, the operator $T_h \in \mathcal{L}(\mathcal{X}_h)$ is bijective. $\square$

**Remark 15.** We recover the first Banach-Nečas-Babuška condition [27, Theorem 25.9, (BNB1)]:

$$ a_{S,h}((\mathbf{u}_h', \mathbf{p}_h'), (\mathbf{v}_h^*, q_h^*)) \ge \alpha_{T_h} (C_{\mathrm{max}}^{nc})^{-1} \|(\mathbf{u}_h', \mathbf{p}_h')\|_{\mathcal{X}_h} \|(\mathbf{v}_h^*, q_h^*)\|_{\mathcal{X}_h}. $$

We can now prove the discrete counterpart of Theorem 7.
We will exhibit $C_{nc}$ for $k=1, d=2, 3$ in §3.1, obtaining $C_{nc} = 1$; and then for $k=2, d=2$ in §3.3, obtaining $C_{nc} = \sigma^2 + 1$. Let us prove Proposition 9.

*Proof.* We follow the proof of Proposition 2. Let us consider $(\mathbf{u}'_h, p'_h) \in \mathcal{X}_h$ and let us build $(\mathbf{v}^*_h, q^*_h) \in \mathcal{X}_h$ satisfying (2.2) (with $V = \mathcal{X}_h$). We need three main steps.

1. According to Prop. 1, it exists $\tilde{\mathbf{v}}_{p'_h} \in \mathbf{V}^{\perp}$ such that $\mathrm{div} \, \tilde{\mathbf{v}}_{p'_h} = p'_h$ in $\Omega$ and $\|\tilde{\mathbf{v}}_{p'_h}\|_{\mathbf{H}_0^1(\Omega)} \le C_{\mathrm{div}} \|p'_h\|_{L^2(\Omega)}$. Let us set $\mathbf{v}_{p'_h} = \nu^{-1} \tilde{\mathbf{v}}_{p'_h}$. Consider $\mathbf{v}_{h,p'_h} = \Pi_{nc} \mathbf{v}_{p'_h}$, for all $q_h \in Q_h$, we have: $(\mathrm{div}_h \mathbf{v}_{h,p'_h}, q_h)_{L^2(\Omega)} = \nu^{-1} (p'_h, q_h)_{L^2(\Omega)}$ and

$$ (2.43) \qquad \|\mathbf{v}_{h,p'_h}\|_h \le \nu^{-1} C_{\mathrm{div}}^{nc} \|p'_h\|_{L^2(\Omega)} \quad \text{where } C_{\mathrm{div}}^{nc} = C_{nc} C_{\mathrm{div}}. $$

Let us set $(\mathbf{v}^*_h, q^*_h) := (\gamma_{nc}\mathbf{u}'_h - \mathbf{v}_{h,p'_h}, -\gamma_{nc}p'_h)$, with $\gamma_{nc} > 0$. We obtain:

$$ (2.44) \qquad a_{S,h}((\mathbf{u}'_h, p'_h), (\mathbf{v}^*_h, q^*_h)) = \nu \gamma_{nc} \|\mathbf{u}'_h\|_h^2 + \nu^{-1} \|p'_h\|_{L^2(\Omega)}^2 - \nu(\mathbf{u}'_h, \mathbf{v}_{h,p'_h})_h. $$

2. In order to bound the last term of (2.44), we use Young inequality and then inequality (2.43) so that for all $\eta_{nc} > 0$:

$$ (2.45) \qquad (\mathbf{u}'_h, \mathbf{v}_{h,p'_h})_h \le \frac{\eta_{nc}}{2} \|\mathbf{u}'_h\|_h^2 + \frac{\eta_{nc}^{-1}}{2} \left( \frac{C_{\mathrm{div}}^{nc}}{\nu} \right)^2 \|p'_h\|_{L^2(\Omega)}^2. $$

3. Using the bound (2.45) in (2.44) and choosing $\eta_{nc} = \gamma_{nc}$, we get:

$$ a_{S,h}((\mathbf{u}'_h, p'_h), (\mathbf{v}^*_h, q^*_h)) \ge \nu \left( \frac{\gamma_{nc}}{2} \|\mathbf{u}'_h\|_h^2 + \nu^{-2} \left( 1 - \frac{(\gamma_{nc})^{-1}}{2} (C_{\mathrm{div}}^{nc})^2 \right) \|p'_h\|_{L^2(\Omega)}^2 \right). $$

Consider now $\gamma_{nc} = (C_{\mathrm{div}}^{nc})^2$. We obtain:

$$ a_{S,h}((\mathbf{u}'_h, p'_h), (\mathbf{v}^*_h, q^*_h)) \ge \alpha_T C_{\mathrm{min}}^{nc} \|(\mathbf{u}'_h, p'_h)\|_{\mathcal{X}_h}^2, \text{ where } C_{\mathrm{min}}^{nc} = \min((C_{\mathrm{div}}^{nc})^2, 1). $$

We obtain (2.42) with $\alpha_{T_h} = \alpha_T C_{\mathrm{min}}^{nc}$. Suppose that $C_{nc} \ge 1$. Then $\alpha_{T_h} = \alpha_T = \frac{\nu}{2}$.

The operator $T_h$ such that $T_h(\mathbf{u}'_h, p'_h) = (\mathbf{v}^*_h, p^*_h)$ is linear and continuous. We have indeed:

$$ \|T_h(\mathbf{u}'_h, p'_h)\|_{\mathcal{X}_h}^2 = \|\mathbf{v}^*_h\|_h^2 + \nu^{-2} \|q^*_h\|_{L^2(\Omega)} \le (C_{\mathrm{max}}^{nc})^2 \|(\mathbf{u}'_h, p'_h)\|_{\mathcal{X}_h}^2 $$

$$ \text{where } C_{\mathrm{max}}^{nc} = C_{\mathrm{div}}^{nc} \left( \max \left( 2 + (C_{\mathrm{div}}^{nc})^2, 2(C_{\mathrm{div}}^{nc})^2 \right) \right)^{\frac{1}{2}}. $$

Remark that, given $(\mathbf{v}^*_h, q^*_h) \in \mathcal{X}_h$, choosing $(\mathbf{u}'_h, p'_h) = (\gamma_{nc}^{-1} \mathbf{v}^*_h - \gamma_{nc}^{-2} \mathbf{v}_{h,q^*_h}, -\gamma_{nc}^{-1} q^*_h)$ yields $T_h((\mathbf{u}'_h, p'_h)) = (\mathbf{v}^*_h, q^*_h)$. Hence, the operator $T_h \in \mathcal{L}(\mathcal{X}_h)$ is bijective. $\square$

**Remark 1.** We recover the first Banach-Nečas-Babuška condition [8, Thm. 25.9, (BNB1)]:

$$ a_{S,h}((\mathbf{u}'_h, p'_h), (\mathbf{v}^*_h, q^*_h)) \ge \alpha_{T_h} (C_{\mathrm{max}}^{nc})^{-1} \|(\mathbf{u}'_h, p'_h)\|_{\mathcal{X}_h} \|(\mathbf{v}^*_h, q^*_h)\|_{\mathcal{X}_h}. $$

We can now prove the discrete counterpart of Theorem 3.

**Theorem 4.** Suppose that it exists a Fortin operator $\Pi_{nc}: \mathbf{H}^1(\Omega) \to \mathbf{X}_h$ satisfying (2.40)-(2.41). Then Problem (2.39) is well-posed. It admits one and only one solution $(\mathbf{u}_h, p_h)$ such that:

$$ (2.46) \qquad \begin{cases} \|\mathbf{u}_h\|_h & \lesssim \nu^{-1} C_0^{nc} \|\mathbf{f}\|_{L^2(\Omega)} \\ \|p_h\|_{L^2(\Omega)} & \lesssim 2 C_0^{nc} C_{\mathrm{div}}^{nc} \|\mathbf{f}\|_{L^2(\Omega)} \\ \\ \|\mathbf{u}_h\|_h & \lesssim \nu^{-1} C_{-1}^{nc} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \\ \|p_h\|_{L^2(\Omega)} & \lesssim 2 C_{-1}^{nc} C_{\mathrm{div}}^{nc} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \end{cases}, $$
**Theorem 16.** Suppose that there exists a Fortin operator $\Pi_{nc}: \mathbf{H}^1(\Omega) \to \mathbf{X}_h$ satisfying (40)-(41). Then Problem (39) is well-posed. It admits one and only one solution $(\mathbf{u}_h, p_h)$ such that:

$$ \begin{array}{ll} \text{if } \mathbf{f} \in \mathbf{L}^2(\Omega): & \|\mathbf{u}_h\|_h \lesssim v^{-1} C_0^{nc} \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)}, \quad \|p_h\|_{L^2(\Omega)} \lesssim 2 C_0^{nc} C_{\text{div}}^{nc} \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)} \\ \text{if } \mathbf{f} \in \mathbf{H}^{-1}(\Omega): & \|\mathbf{u}_h\|_h \lesssim v^{-1} C_{\mathcal{F}_h}^{nc} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)}, \quad \|p_h\|_{L^2(\Omega)} \lesssim 2 C_{\mathcal{F}_h}^{nc} C_{\text{div}}^{nc} \|\mathbf{f}\|_{\mathbf{H}^{-1}(\Omega)} \end{array} \quad (46) $$

where $C_0^{nc} = \sigma C_0 h_\Omega$.

**Proof.** Consider $(\mathbf{u}_h, p_h)$ the unique solution of Problem (39). Choosing $v_h = 0$, we obtain that $\text{div}_h \mathbf{u}_h = 0$. Let $\mathbf{f} \in \mathbf{L}^2(\Omega)$. Now, choosing $v_h = \mathbf{u}_h$ in (39), using Cauchy-Schwarz inequality, we get that: $\|\mathbf{u}_h\|_h \le v^{-1} \sigma C_0 h_\Omega \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)}$ using inequality (31). Consider $(v_h, q_h) = (v_{h,p_h}, 0)$ in (39), where $v_{h,p_h} = \Pi_{nc}v_{p_h}$ is built as $v_{h,p'_h}$ in point 1, setting $p'_h = p_h$. Suppose that $\mathbf{f} \in \mathbf{L}^2(\Omega)$. Notice that $v^{-1} \|p_h\|_{L^2(\Omega)}^2 = v(\mathbf{u}_h, v_{h,p_h})_h - (\mathbf{f}, v_{h,p_h})_{L^2(\Omega)}$. Using Cauchy-Schwarz inequality, we have: $v^{-1} \|p_h\|_{L^2(\Omega)}^2 \le v \|\mathbf{u}_h\|_h \|v_{h,p_h}\|_h + \|\mathbf{f}\|_{L^2(\Omega)} \|v_{h,p_h}\|_{L^2(\Omega)}$. Using Poincaré-Steklov inequality (31), hypothesis (40), and the previous estimate on $\|\mathbf{u}_h\|_h$, we have:

$$ \|p_h\|_{L^2(\Omega)}^2 \lesssim 2\sigma C_0 h_\Omega \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)} \|v_{h,p_h}\|_h \lesssim 2\sigma C_0 h_\Omega C_{\text{div}}^{nc} \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)} \|p_h\|_{L^2(\Omega)}. $$

Let $\mathbf{f} \in \mathbf{H}^{-1}(\Omega)$. We apply the same reasoning, using inequality (38). $\square$

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

As a corollary of Theorem 16, the following a priori error estimates hold:

**Corollary 17 ([22, Theorems 3, 4, 6], [29, Equation (47)]).** Under the assumption of Theorem 16, suppose that $(\mathbf{u}, p) \in (\mathbf{H}^{1+k}(\Omega) \cap \mathbf{H}_0^1(\Omega)) \times (H^k(\Omega) \cap L_{zmv}^2(\Omega))$, we have the estimates:

$$ \|\mathbf{u} - \mathbf{u}_h\|_h \lesssim \sigma^l h^k (|\mathbf{u}|_{\mathbf{H}^{k+1}(\Omega)} + v^{-1} |p|_{H^k(\Omega)}), \quad (47) $$

$$ v^{-1} \|p - p_h\|_{L^2(\Omega)} \lesssim \sigma^l h^k (|\mathbf{u}|_{\mathbf{H}^{k+1}(\Omega)} + v^{-1} |p|_{H^k(\Omega)}). \quad (48) $$

Suppose moreover that the domain $\Omega$ is convex. Then we have:

$$ \|\mathbf{u} - \mathbf{u}_h\|_{\mathbf{L}^2(\Omega)} \lesssim \sigma^{2l} h^{k+1} (|\mathbf{u}|_{\mathbf{H}^{k+1}(\Omega)} + v^{-1} |p|_{H^k(\Omega)}). \quad (49) $$

The hidden constants depend on $k$ but they don't depend on the mesh. The parameter $\sigma$ is the shape regularity parameter and the exponent $l \in \mathbb{N}^*$ depends on $k$. When $k=1, d=2, 3$, we have: $l=1$ and when $k=2, d=2$, we have: $l=2$.

When $\Omega$ is not convex, the exponent on $h$ in Equation (49) is equal to $k+s$ where $s \in ]0, 1[$ depends on $\Omega$ (cf. [27, Theorem 31.33]).

The main issue with nonconforming mixed finite elements is the construction of the basis functions. In a recent paper, Sauter explains such a construction in two dimensions [41, Theorem 1.3], and gives a bound to the discrete counterpart $\beta_{\mathcal{F}}(\Omega)$ of $\beta(\Omega)$ defined in (7):

$$ \beta_{\mathcal{F}}(\Omega) = \inf_{q_h \in Q_h \setminus \{0\}} \sup_{v_h \in X_{0,h} \setminus \{0\}} \frac{(\text{div}_h v_h, q_h)_{L^2(\Omega)}}{\|q_h\|_{L^2(\Omega)} \|v_h\|_h} \ge c_{\mathcal{F}} (\log(k+1))^{-\alpha}, \quad (50) $$

where the parameter $\alpha$ is explicit and depends on $k$ and on the mesh topology; and the constant $c_{\mathcal{F}}$ depends only on the shape-regularity of the mesh.
where $C_0^{nc} = \sigma C_\Omega h_\Omega$ and $C_{-1}^{nc} = C_{\mathcal{I}_h}^{nc} C_{PS} h_\Omega$.

*Proof.* Consider $(\mathbf{u}_h, p_h)$ the unique solution of Problem (2.39). Choosing $\mathbf{v}_h = 0$, we obtain that $\mathrm{div}_h \mathbf{u}_h = 0$. Let $\mathbf{f} \in \mathbf{L}^2(\Omega)$. Now, choosing $\mathbf{v}_h = \mathbf{u}_h$ in (2.39), using Cauchy-Schwarz inequality, we get that: $\|\mathbf{u}_h\|_h \le \nu^{-1} \sigma C_\Omega h_\Omega \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)}$ using inequality (2.30). Consider $(\mathbf{v}_h, q_h) = (\mathbf{v}_{h,p_h}, 0)$ in (2.39), where $\mathbf{v}_{h,p_h} = \Pi_{nc}\mathbf{v}_{p_h}$ is built as $\mathbf{v}_{h,p'_h}$ in point 1, setting $p'_h = p_h$. Suppose that $\mathbf{f} \in \mathbf{L}^2(\Omega)$. Notice that $\nu^{-1} \|p_h\|_{L^2(\Omega)}^2 = \nu(\mathbf{u}_h, \mathbf{v}_{h,p_h})_h - (\mathbf{f}, \mathbf{v}_{h,p_h})_{\mathbf{L}^2(\Omega)}$. Using Cauchy-Schwarz inequality, we have: $\nu^{-1} \|p_h\|_{L^2(\Omega)}^2 \le \nu \|\mathbf{u}_h\|_h \|\mathbf{v}_{h,p_h}\|_h + \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)} \|\mathbf{v}_{h,p_h}\|_{\mathbf{L}^2(\Omega)}$. Using Poincaré-Steklov inequality (2.30), *Hypothesis (2.40)*, and the previous estimate on $\|\mathbf{u}_h\|_h$, we have:

$$ \|p_h\|_{L^2(\Omega)}^2 \lesssim 2 \sigma C_\Omega h_\Omega \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)} \|\mathbf{v}_{h,p_h}\|_h \lesssim 2 \sigma C_\Omega h_\Omega C_{\mathrm{div}}^{nc} \|\mathbf{f}\|_{\mathbf{L}^2(\Omega)} \|p_h\|_{L^2(\Omega)}. $$

Let $\mathbf{f} \in \mathbf{H}^{-1}(\Omega)$. We apply the same reasoning, using inequalities (2.38) and (2.30). $\square$

As a corollary of Theorem 4, the following a priori error estimate follows [2, Theorems 3, 4, 6] and [3, Eq. (47)]:

**Corollary 1.** Under the assumption of Theorem 4, suppose that $(\mathbf{u}, p) \in (\mathbf{H}^{1+k}(\Omega) \cap \mathbf{H}_0^1(\Omega)) \times (H^k(\Omega) \cap L_{zmv}^2(\Omega))$, we have the estimates:

$$ (2.47) \qquad \|\mathbf{u} - \mathbf{u}_h\|_h \lesssim \sigma^l h^k (|\mathbf{u}|_{\mathbf{H}^{k+1}(\Omega)} + \nu^{-1} |p|_{H^k(\Omega)}), $$

$$ (2.48) \qquad \nu^{-1} \|p - p_h\|_{L^2(\Omega)} \lesssim \sigma^l h^k (|\mathbf{u}|_{\mathbf{H}^{k+1}(\Omega)} + \nu^{-1} |p|_{H^k(\Omega)}), $$

Suppose moreover that the domain $\Omega$ is convex. Then we have:

$$ (2.49) \qquad \|\mathbf{u} - \mathbf{u}_h\|_{L^2(\Omega)} \lesssim \sigma^{2l} h^{k+1} (|\mathbf{u}|_{\mathbf{H}^{k+1}(\Omega)} + \nu^{-1} |p|_{H^k(\Omega)}), $$

The hidden constants depend on $k$ but they don't dependent on mesh. The parameter $\sigma$ is the shape regularity parameter and the exponent $l \in \mathbb{N}^*$ depends on $k$. When $k=1, d=2, 3$, we have: $l=1$ and when $k=2, d=2$, we have: $l=2$.

When $\Omega$ is not convex, the exponent on $h$ in Equation (2.49) is equal to $k+s$ where $s \in ]0, 1[$ depends on $\Omega$ (cf. [8, Theorem 31.33]).

The main issue with nonconforming mixed finite elements is the construction the basis functions. In a recent paper, Sauter explains such a construction in two dimensions [23, Theorem 1.3], and gives a bound to the discrete counterpart $\beta_T(\Omega)$ of $\beta(\Omega)$ defined in (2.6):

$$ (2.50) \qquad \beta_T(\Omega) = \inf_{q_h \in Q_h \setminus \{0\}} \sup_{\mathbf{v}_h \in \mathbf{X}_{0,h} \setminus \{0\}} \frac{(\mathrm{div}_h \mathbf{v}_h, q_h)_{L^2(\Omega)}}{\|q_h\|_{L^2(\Omega)} \|\mathbf{v}_h\|_h} \ge c_T (\log(k+1))^{-\alpha}, $$

where the parameter $\alpha$ is explicit and depends on $k$ and on the mesh topology; and the constant $c_T$ depends only on the shape-regularity of the mesh.

### 3. EXAMPLES OF NON CONFORMING DISCRETIZATION FOR STOKES PROBLEM

#### 3.1. Nonconforming Crouzeix-Raviart mixed finite elements for $k=1$.
We study the lowest order nonconforming Crouzeix-Raviart mixed finite elements
### 3. Examples of nonconforming discretization for Stokes problem

#### 3.1. Nonconforming Crouzeix-Raviart mixed finite elements for $k = 1$

We study the lowest order nonconforming Crouzeix-Raviart mixed finite elements [22]. Let us consider $X_{CR}$ (resp. $X_{0,CR}$), the space of nonconforming approximation of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order 1:

$$X_{CR} = \left\{ v_h \in P_{disc}^1(\mathcal{T}_h) \mid \forall f \in \mathcal{J}_F^i, \int_{F_f} [v_h]_{F_f} = 0 \right\}; \qquad (51)$$

$$X_{0,CR} = \left\{ v_h \in X_{CR} \mid \forall f \in \mathcal{J}_F^b, \int_{F_f} v_h = 0 \right\}.$$

The space of nonconforming approximation of of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order 1 is $X_{CR} = (X_{CR})^d$ (resp. $X_{0,CR} = (X_{0,CR})^d$). We set $\mathcal{X}_{CR} := X_{0,CR} \times Q_{CR}$ where $Q_{CR} = P_{disc}^0(\mathcal{T}_h) \cap L_{zmv}^2(\Omega)$.

We can endow $X_{CR}$ with the basis $(\psi_f)_{f \in \mathcal{J}_F}$ such that:

$$\forall \ell \in \mathcal{J}_K, \psi_f|_{K_\ell} = \begin{cases} 1 - d\lambda_{i,\ell} & \text{if } f \in \mathcal{J}_{F,\ell}, \\ 0 & \text{otherwise,} \end{cases}$$

where $i \in \mathcal{J}_S$ is such that $S_i$ is the vertex opposite to $F_f$ in $K_\ell$. We then have $\psi_{f|F_f} = 1$, so that $[\psi_f]_{F_f} = 0$ if $f \in \mathcal{J}_F^i$, and for all $f, f' \in \mathcal{J}_F$, $f' \neq f, \int_{F_f'} \psi_f = 0$.

We have: $X_{CR} = \text{vect}((\psi_f)_{f \in \mathcal{J}_F})$ and $X_{0,CR} = \text{vect}((\psi_f)_{f \in \mathcal{J}_F^i})$.

The Crouzeix-Raviart interpolation operator $\pi_{CR}$ for scalar functions is defined by:

$$\pi_{CR} : \left\{ \begin{array}{l} H^1(\Omega) \to X_{CR} \\ v \mapsto \sum_{f \in \mathcal{J}_F} \pi_f v \psi_f, \text{ where } \pi_f v = \frac{1}{|F_f|} \int_{F_f} v. \end{array} \right.$$

Notice that $\forall f \in \mathcal{J}_F$, $\int_{F_f} \pi_{CR} v = \int_{F_f} v$. Moreover, the Crouzeix-Raviart interpolation operator preserves the constants, so that $\pi_{CR} v_{|\Omega} = v_{|\Omega}$ where $v_{|\Omega} = \int_{\Omega} v/|\Omega|$. We recall that for $k = 1$, the coefficient $C_{nc}$ in (40) is equal to 1:

**Lemma 18 ([3, Lemma 2]).** The Crouzeix-Raviart interpolation operator $\pi_{CR}$ is such that:

$$\forall v \in H^1(\Omega), \quad \|\pi_{CR} v\|_h \le \|\text{grad } v\|_{L^2(\Omega)}. \qquad (52)$$

**Proof.** We have, integrating by parts twice and using Cauchy-Schwarz inequality:

$$\begin{align*} \text{grad } \pi_{CR} v|_{K_\ell} &= |K_\ell|^{-1} \int_{K_\ell} \text{grad } \pi_{CR} v, \\ &= |K_\ell|^{-1} \sum_{f \in \mathcal{J}_{F,\ell}} \int_{F_f} \pi_{CR} v \mathbf{n}_f, \\ &= |K_\ell|^{-1} \sum_{f \in \mathcal{J}_{F,\ell}} \int_{F_f} v \mathbf{n}_f, \\ &= |K_\ell|^{-1} \int_{K_\ell} \text{grad } v, \end{align*}$$

$$|\text{grad } \pi_{CR} v|_{K_\ell} \le |K_\ell|^{-\frac{1}{2}} \|\text{grad } v\|_{L^2(K_\ell)}$$

$$\Rightarrow \quad \|\text{grad } \pi_{CR} v\|_{L^2(K_\ell)}^2 \le \|\text{grad } v\|_{L^2(K_\ell)}^2.$$

Summing these local estimates over $\ell \in \mathcal{J}_K$, we obtain (52). □

For a vector $v \in H^1(\Omega)$ of components $(v_{d'})^d_{d'=1}$, the Crouzeix-Raviart interpolation operator is such that: $\Pi_{CR} v = (\pi_{CR} v_{d'})^d_{d'=1}$. Let us set $\Pi_f v = (\pi_f v_{d'})^d_{d'=1}$.
[2]. Let us consider $X_{CR}$ (resp. $X_{0,CR}$), the space of nonconforming approximation of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order 1:

$$ (3.1) \qquad \begin{aligned} X_{CR} &= \left\{ v_h \in P_{disc}^1(\mathcal{T}_h); \quad \forall f \in \mathcal{I}_F^i, \int_{F_f} [v_h]_{F_f} = 0 \right\}; \\ X_{0,CR} &= \left\{ v_h \in X_{CR}; \quad \forall f \in \mathcal{I}_F^b, \int_{F_f} v_h = 0 \right\}. \end{aligned} $$

The space of nonconforming approximation of of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order 1 is $X_{CR} = (X_{CR})^d$ (resp. $X_{0,CR} = (X_{0,CR})^d$). We set $\mathcal{X}_{CR} := X_{0,CR} \times Q_{CR}$ where $Q_{CR} = P_{disc}^0(\mathcal{T}_h) \cap L_{zmv}^2(\Omega)$.

We can endow $X_{CR}$ with the basis $(\psi_f)_{f \in \mathcal{I}_F}$ such that: $\forall \ell \in \mathcal{I}_K$,

$$ \psi_{f|K_{\ell}} = \begin{cases} 1 - d\lambda_{i,\ell} & \text{if } f \in \mathcal{I}_{F,\ell}, \\ 0 & \text{otherwise,} \end{cases} $$

where $S_i$ is the vertex opposite to $F_f$ in $K_\ell$. We then have $\psi_{f|F_f} = 1$, so that $[\psi_f]_{F_f} = 0$ if $f \in \mathcal{I}_F^i$, and for all $f, f' \in \mathcal{I}_F$, $f' \neq f$, $\int_{F_{f'}} \psi_f = 0$.

We have: $X_{CR} = \text{vec}((\psi_f)_{f \in \mathcal{I}_F})$ and $X_{0,CR} = \text{vec}((\psi_f)_{f \in \mathcal{I}_F^b})$.

The Crouzeix-Raviart interpolation operator $\pi_{CR}$ for scalar functions is defined by:

$$ \pi_{CR} : \left\{ \begin{array}{ll} H^1(\Omega) & \rightarrow X_{CR} \\ v & \mapsto \sum_{f \in \mathcal{I}_F} \pi_f v \psi_f, \quad \text{where } \pi_f v = \frac{1}{|F_f|} \int_{F_f} v. \end{array} \right. $$

Notice that $\forall f \in \mathcal{I}_F$, $\int_{F_f} \pi_{CR}v = \int_{F_f} v$. Moreover, the Crouzeix-Raviart interpolation operator preserves the constants, so that $\pi_{CR}v_{\Omega} = v_{\Omega}$ where $v_{\Omega} = \int_{\Omega} v/|\Omega|$. We recall that for $k=1$, the coefficient $C_{nc}$ in (2.40) is equal to 1, as proven in [24, Lemma 2]:

**Lemma 2.** The Crouzeix-Raviart interpolation operator $\pi_{CR}$ is such that:

$$ (3.2) \qquad \forall v \in H^1(\Omega), \quad \|\pi_{CR}v\|_h \le \|\text{grad } v\|_{L^2(\Omega)}. $$

*Proof.* We have, integrating by parts twice and using Cauchy-Schwarz inequality:

$$ \begin{aligned} \text{grad } \pi_{CR} v|_{K_{\ell}} &= |K_{\ell}|^{-1} \int_{K_{\ell}} \text{grad } \pi_{CR} v = |K_{\ell}|^{-1} \sum_{f \in \mathcal{I}_{F,\ell}} \int_{F_f} \pi_{CR} v \mathbf{n}_f, \\ &= |K_{\ell}|^{-1} \sum_{f \in \mathcal{I}_{F,\ell}} \int_{F_f} v \mathbf{n}_f = |K_{\ell}|^{-1} \int_{K_{\ell}} \text{grad } v, \\ |\text{grad } \pi_{CR} v|_{K_{\ell}}| &\le |K_{\ell}|^{-\frac{1}{2}} \|\text{grad } v\|_{L^2(K_{\ell})} \\ \Rightarrow \quad \|\text{grad } \pi_{CR} v\|_{L^2(K_{\ell})}^2 &\le \|\text{grad } v\|_{L^2(K_{\ell})}^2. \end{aligned} $$

Summing these local estimates over $\ell \in \mathcal{I}_K$, we obtain (3.2). □

For a vector $v \in H^1(\Omega)$ of components $(v'_{d'=1})^d_{d'=1}$, the Crouzeix-Raviart interpolation operator is such that: $\Pi_{CR}v = (\pi_{CR}v'_{d'=1})^d_{d'=1}$. Let us set $\Pi_f v = (\pi_f v'_{d'=1})^d_{d'=1}$.
**Lemma 19.** The Crouzeix-Raviart interpolation operator $\Pi_{CR}$ can play the role of the Fortin operator:

$$ \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \quad \|\Pi_{CR}\mathbf{v}\|_h \le \|\mathbf{Grad}\mathbf{v}\|_{\ell^2(\Omega)}, \quad (53) $$

$$ \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \quad (\mathrm{div}_h \Pi_{CR}\mathbf{v}, q_h)_{L^2(\Omega)} = (\mathrm{div}\mathbf{v}, q_h)_{L^2(\Omega)}, \quad \forall q \in Q_h, \quad (54) $$

Moreover, for all $\mathbf{v} \in \mathbf{P}^1(\Omega)$, $\Pi_{CR}\mathbf{v} = \mathbf{v}$.

**Proof.** We obtain (53) applying Lemma 18 component by component. By integrating by parts, we have $\forall \mathbf{v} \in \mathbf{H}^1(\Omega)$, $\forall \ell \in \mathcal{I}_K$:

$$ \int_{K_{\ell}} \mathrm{div} \Pi_{CR} \mathbf{v} = \sum_{f \in \mathcal{F}_{F, \ell}} \int_{F_f} \Pi_{CR} \mathbf{v} \cdot \mathbf{n}_f = \sum_{f \in \mathcal{F}_{F, \ell}} \int_{F_f} \Pi_f \mathbf{v} \cdot \mathbf{n}_f = \sum_{f \in \mathcal{F}_{F, \ell}} \int_{F_f} \mathbf{v} \cdot \mathbf{n}_f = \int_{K_{\ell}} \mathrm{div} \mathbf{v}, $$

so that (54) is satisfied.

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

We can apply the T-coercivity theory to show the following result:

**Theorem 20.** Let $\mathcal{X}_h = \mathcal{X}_{CR}$. Then the continuous bilinear form $a_{S,h}(\cdot, \cdot)$ is $T_h$-coercive and Problem (39) is well-posed.

**Proof.** Using estimates (53) and (31), we apply the proof of Theorem 16.

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

Since the constant of the interpolation operator $\Pi_{CR}$ is equal to 1, we have $C_{\min}^{nc} = \min((C_{\mathrm{div}})^2, 1) = 1$ and $C_{\max}^{nc} = C_{\max}$: the stability constant of the nonconforming Crouzeix-Raviart mixed finite elements is independent of the mesh. This is not the case for higher order (see [14, Theorem 2.2]).

## 3.2. Comments on higher-order methods

For higher order, we cannot build the interpolation operator component by component, since higher-order divergence moments must be preserved. Thus, for $k > 1$, we must build $\Pi_{nc}$ so that for all $\mathbf{v} \in \mathbf{H}^1(\Omega)$, for all $\ell \in \mathcal{I}_K$, for all $q \in P^{k-1}(K_\ell)$:

$$ \int_{K_{\ell}} q \mathrm{div} \Pi_{nc} \mathbf{v} = \int_{K_{\ell}} q \mathrm{div} \mathbf{v}. $$

We recall that by integration by parts, we have:

$$ \int_{K_{\ell}} q \mathrm{div} \Pi_{nc} \mathbf{v} + \int_{K_{\ell}} \mathbf{grad} q \cdot \Pi_{nc} \mathbf{v} = \int_{\partial K_{\ell}} q \Pi_{nc} \mathbf{v} \cdot \mathbf{n}|_{\partial K_{\ell}}. \quad (55) $$

Hence, to obtain a local estimate of $\|\mathbf{Grad} \Pi_{nc} \mathbf{v}\|_{\ell^2(K_\ell)}$, we will need the following lemma:

**Lemma 21.** Let $\mathbf{v} \in \mathbf{H}^1(K_\ell)$ and $q \in P^{k-1}(K_\ell) \cap L_{zmv}^2(K_\ell)$. We have:

$$ \left| \int_{\partial K_{\ell}} q (\mathbf{v} - \underline{\mathbf{v}}_{\ell}) \cdot \mathbf{n}|_{\partial K_{\ell}} \right| \le (\sqrt{d} + 1) \pi^{-1} h_{\ell} \| \mathbf{grad} q \|_{L^2(K_{\ell})} \| \mathbf{Grad} \mathbf{v} \|_{\ell^2(K_{\ell})} \quad (56) $$

**Proof.** We have by integration by parts, and then using Cauchy-Schwarz inequality:

$$ \begin{aligned} \left| \int_{\partial K_{\ell}} q (\mathbf{v} - \underline{\mathbf{v}}_{\ell}) \cdot \mathbf{n}|_{\partial K_{\ell}} \right| & \le \left| \int_{K_{\ell}} q \mathrm{div}(\mathbf{v} - \underline{\mathbf{v}}_{\ell}) \right| + \left| \int_{K_{\ell}} \mathbf{grad} q \cdot (\mathbf{v} - \underline{\mathbf{v}}_{\ell}) \right|, \\ & \le \sqrt{d} \|q\|_{L^2(K_{\ell})} \|\mathbf{Grad}(\mathbf{v} - \underline{\mathbf{v}}_{\ell})\|_{\ell^2(K_{\ell})} + \|\mathbf{grad} q\|_{L^2(K_{\ell})} \|(\mathbf{v} - \underline{\mathbf{v}}_{\ell})\|_{L^2(K_{\ell})}, \\ & \le (\sqrt{d} + 1) \pi^{-1} h_{\ell} \|\mathbf{grad} q\|_{L^2(K_{\ell})} \|\mathbf{Grad} \mathbf{v}\|_{\ell^2(K_{\ell})}, \text{ using (28) twice.} \end{aligned} $$

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

Noticing that for $\mathbf{v} \in \mathbf{V}_h$, for all $\ell \in \mathcal{I}_K$, $\mathrm{div}\mathbf{v}|_{K_\ell} = 0$, we obtain in the way the following result.
**Lemma 3.** The Crouzeix-Raviart interpolation operator $\Pi_{CR}$ can play the role of the Fortin operator:

$$ (3.3) \qquad \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \quad \|\Pi_{CR}\mathbf{v}\|_h \le \|\mathbf{Grad}\,\mathbf{v}\|_{L^2(\Omega)}, $$

$$ (3.4) \qquad \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \quad (\mathrm{div}_h \Pi_{CR}\mathbf{v}, q_h)_{L^2(\Omega)} = (\mathrm{div}\,\mathbf{v}, q_h)_{L^2(\Omega)}, \quad \forall q \in Q_h, $$

Moreover, for all $\mathbf{v} \in \mathbf{P}^1(\Omega)$, $\Pi_{CR}\mathbf{v} = \mathbf{v}$.

*Proof.* We obtain (3.3) applying Lemma 2 component by component. By integrating by parts, we have $\forall \mathbf{v} \in \mathbf{H}^1(\Omega), \forall \ell \in \mathcal{I}_K$:

$$ \begin{aligned} \int_{K_{\ell}} \mathrm{div} \, \Pi_{CR} \mathbf{v} &= \sum_{f \in \mathcal{I}_{F, \ell}} \int_{F_f} \Pi_{CR} \mathbf{v} \cdot \mathbf{n}_f = \sum_{f \in \mathcal{I}_{F, \ell}} \int_{F_f} \Pi_f \mathbf{v} \cdot \mathbf{n}_f \\ &= \sum_{f \in \mathcal{I}_{F, \ell}} \int_{F_f} \mathbf{v} \cdot \mathbf{n}_f = \int_{K_{\ell}} \mathrm{div} \, \mathbf{v} \end{aligned} $$

so that (3.4) is satisfied.

□

We can apply the T-coercivity theory to show the following result:

**Theorem 5.** Let $\mathcal{X}_h = \mathcal{X}_{CR}$. Then the continuous bilinear form $a_{S,h}(\cdot, \cdot)$ is $T_h$-coercive and Problem (2.39) is well-posed.

*Proof.* Using estimates (3.3) and (2.30), we apply the proof of Theorem 4. □

Since the constant of the interpolation operator $\Pi_{CR}$ is equal to 1, we have $C_{\min}^{nc} = \min((C_{\mathrm{div}})^2, 1) = 1$ and $C_{\max}^{nc} = C_{\max}$: the stability constant of the non-conforming Crouzeix-Raviart mixed finite elements is independent of the mesh. This is not the case for higher-order (see [25, Theorem 2.2]).

**3.2. Comments on higher-order methods.** For higher-order, we cannot build the interpolation operator component by component, since higher-order divergence moments must be preserved. Thus, for $k > 1$, we must build $\Pi_{nc}$ so that for all $\mathbf{v} \in \mathbf{H}^1(\Omega)$, for all $\ell \in \mathcal{I}_K$, for all $q \in P^{k-1}(K_\ell)$:

$$ \int_{K_{\ell}} q \mathrm{div} \, \Pi_{nc} \mathbf{v} = \int_{K_{\ell}} q \mathrm{div} \, \mathbf{v}. $$

We recall that by integration by parts, we have:

$$ (3.5) \qquad \int_{K_{\ell}} q \mathrm{div} \, \Pi_{nc} \mathbf{v} + \int_{K_{\ell}} \mathbf{grad} \, q \cdot \Pi_{nc} \mathbf{v} = \int_{\partial K_{\ell}} q \Pi_{nc} \mathbf{v} \cdot \mathbf{n}|_{\partial K_{\ell}}. $$

Hence, to obtain a local estimate of $\|\mathbf{Grad} \, \Pi_{nc} \mathbf{v}\|_{L^2(K_\ell)}$, we will need the following Lemma:

**Lemma 4.** Let $\mathbf{v} \in \mathbf{H}^1(K_\ell)$ and $q \in P^{k-1}(K_\ell) \cap L_{zmv}^2(K_\ell)$. We have:

$$ (3.6) \qquad \left| \int_{\partial K_{\ell}} q (\mathbf{v} - \underline{\mathbf{v}}_{\ell}) \cdot \mathbf{n}|_{\partial K_{\ell}} \right| \le (\sqrt{d} + 1) \pi^{-1} h_{\ell} \| \mathbf{grad} \, q \|_{L^2(K_{\ell})} \| \mathbf{Grad} \, \mathbf{v} \|_{L^2(K_{\ell})} $$
**Corollary 22.** Let $\mathbf{v}_h \in \mathbf{V}_h$ and $q \in H^1(K_\ell)$. We have:

$$ \left| \int_{\partial K_\ell} q (\mathbf{v}_h - \mathbf{v}_{h,\ell}) \cdot \mathbf{n}|_{\partial K_\ell} \right| \le \pi^{-1} h_\ell \| \mathbf{grad} q \|_{L^2(K_\ell)} \| \mathbf{Grad} \mathbf{v}_h \|_{L^2(K_\ell)} \quad (57) $$

In the next section, we will see that for $k=2$, $d=2$, we will need Lemma 21. For $k \ge 3$, it could be necessary to bound the tangential components of $\mathbf{v} - \mathbf{v}_\ell$. To do so, we would need to preserve curl integrals on $K_\ell$. Indeed, by integration by parts, we have:

$$ \text{For } d=2, \mathbf{v} \in \mathbf{H}^1(\Omega), q \in P^{k-1}(K_\ell): \quad \int_{K_\ell} q (\mathbf{curl} q \cdot \mathbf{v} - \mathbf{curl} \mathbf{v} q) = \int_{\partial K_\ell} q \mathbf{v} \times \mathbf{n}|_{\partial K_\ell}. $$

$$ \text{For } d=3, \mathbf{v} \in \mathbf{H}^1(\Omega), \mathbf{w} \in \mathbf{P}^{k-1}(K_\ell): \quad \int_{K_\ell} (\mathbf{w} \cdot \mathbf{curl}\mathbf{v} - \mathbf{curl}\mathbf{w} \cdot \mathbf{v}) = \int_{\partial K_\ell} (\mathbf{n}|_{\partial K_\ell} \times \mathbf{v} \times \mathbf{n}|_{\partial K_\ell}) \cdot (\mathbf{w} \times \mathbf{n}|_{\partial K_\ell}). $$

### 3.3. Fortin-Soulie mixed finite elements

We consider here the case $k=2$, $d=2$ and we study the so-called Fortin-Soulie mixed finite elements [29]. We consider a shape-regular triangulation sequence $(\mathcal{T}_h)_h$.

Let us consider $X_{FS}$ (resp. $X_{0,FS}$), the space of nonconforming approximation of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order 2:

$$ X_{FS} = \left\{ v_h \in P_{disc}^2(\mathcal{T}_h) \mid \forall f \in \mathcal{F}_F^i, \forall q_h \in P^1(F_f), \int_{F_f} [v_h]_{F_f} q_h = 0 \right\}; \quad (58) $$

$$ X_{0,FS} = \left\{ v_h \in X_{FS} \mid \forall f \in \mathcal{F}_F^b, \forall q_h \in P^1(F_f), \int_{F_f} v_h q_h = 0 \right\}. $$

The space of nonconforming approximation of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order 2 is $X_{FS} = (X_{FS})^2$ (resp. $X_{0,FS} = (X_{0,FS})^2$). We set $\mathcal{X}_{FS} = X_{0,FS} \times Q_{FS}$ where $Q_{FS} := P_{disc}(\mathcal{T}_h) \cap L_{zmv}^2(\Omega)$.

The building of a basis for $X_{0,FS}$ is more involved than for $X_{0,CR}$ since we cannot use two points per facet as degrees of freedom. Indeed, for all $\ell \in \mathcal{I}_K$, there exists a polynomial of order 2 vanishing on the Gauss-Legendre points of the facets of the boundary $\partial K_\ell$. Let $f \in \mathcal{I}_F$. The barycentric coordinates of the two Gauss-Legendre points $(p_{+,f}, p_{-,f})$ on $F_f$ are such that:

$$ p_{+,f} = (c_+, c_-), \quad p_{-,f} = (c_-, c_+), \quad \text{where } c_\pm = (1 \pm 1/\sqrt{3})/2. $$

These points can be used to integrate exactly order three polynomials:

$$ \forall g \in P^3(F_f), \int_{F_f} g = \frac{|F_f|}{2} (g(p_{+,f}) + g(p_{-,f})). $$

For all $\ell \in \mathcal{I}_K$, we define the quadratic function $\phi_{K_\ell}$ that vanishes on the six Gauss-Legendre points of the facets of $K_\ell$ (see Figure 1):

$$ \phi_{K_\ell} := 2 - 3 \sum_{i \in \mathcal{I}_{S,\ell}} \lambda_{i,\ell}^2 \text{ such that } \forall f \in \mathcal{F}_{F,\ell}, \forall q \in P^1(F_f), \quad \int_{F_f} \phi_{K_\ell} q = 0. \quad (59) $$

We consider the set of the elliptic functions $\phi_{K_\ell}$:

$$ \Phi_h := \{\phi_h \in L^2(\Omega) \mid \forall \ell \in \mathcal{I}_K, \phi_h|_{K_\ell} = v_{K_\ell} \phi_{K_\ell}, v_{K_\ell} \in \mathbb{R}\}. \quad (60) $$

We also define the spaces of $P^2$-Lagrange functions:

$$ X_{LG} := \{v_h \in H^1(\Omega) \mid \forall \ell \in \mathcal{I}_K, v_h|_{K_\ell} \in P^2(K_\ell)\}, \quad X_{0,LG} := \{v_h \in X_{LG} \mid v_h|_{\partial\Omega} = 0\}. $$

The Proposition below allows to build a basis for $X_{0,FS}$:

**Proposition 23 ( [29, Proposition 1])**. We have the following decomposition: $X_{FS} = X_{LG} + \Phi_h$ with $\dim(X_{LG} \cap \Phi_h) = 1$. Any function of $X_{FS}$ can be written as the sum of a function of $X_{LG}$ and a function of $\Phi_h$. This representation can be made unique by specifying one degree of freedom.
Proof. We have by integration by parts, and then using Cauchy-Schwarz inequality:

$$
\begin{align*}
\left| \int_{\partial K_{\ell}} q (\mathbf{v} - \mathbf{v}_{\ell}) \cdot \mathbf{n}_{|\partial K_{\ell}} \right| &\le \left| \int_{K_{\ell}} q \operatorname{div}(\mathbf{v} - \mathbf{v}_{\ell}) \right| + \left| \int_{K_{\ell}} \operatorname{grad} q \cdot (\mathbf{v} - \mathbf{v}_{\ell}) \right|, \\
&\le \sqrt{d} \|q\|_{L^2(K_{\ell})} \| \operatorname{Grad}(\mathbf{v} - \mathbf{v}_{\ell}) \|_{L^2(K_{\ell})} + \| \operatorname{grad} q \|_{L^2(K_{\ell})} \|(\mathbf{v} - \mathbf{v}_{\ell})\|_{L^2(K_{\ell})}, \\
&\le (\sqrt{d} + 1) \pi^{-1} h_{\ell} \| \operatorname{grad} q \|_{L^2(K_{\ell})} \| \operatorname{Grad} \mathbf{v} \|_{L^2(K_{\ell})}, \text{ using (2.26) twice.}
\end{align*}
$$

□

In the next section, we will see that for $k = 2$, $d = 2$, we will need Lemma 4.
For $k \ge 3$, it could be necessary to bound the tangential components of $\mathbf{v} - \mathbf{v}_\ell$. To
do so, we would need to preserve curl integrals on $K_\ell$. Indeed, by integration by
parts, we have:

• For $d = 2$, $\mathbf{v} \in \mathbf{H}^1(\Omega)$ and $q \in P^{k-1}(K_\ell)$:

$$
\int_{K_{\ell}} q (\mathbf{curl} q \cdot \mathbf{v} - \mathbf{curl} \mathbf{v} q) = \int_{\partial K_{\ell}} q \mathbf{v} \times \mathbf{n}_{|\partial K_{\ell}}.
$$

• For $d = 3$, $\mathbf{v} \in \mathbf{H}^1(\Omega)$ and $\mathbf{w} \in \mathbf{P}^{k-1}(K_\ell)$:

$$
\int_{K_{\ell}} (\mathbf{w} \cdot \mathrm{curl} \mathbf{v} - \mathrm{curl} \mathbf{w} \cdot \mathbf{v}) = \int_{\partial K_{\ell}} (\mathbf{n}_{|\partial K_{\ell}} \times \mathbf{v} \times \mathbf{n}_{|\partial K_{\ell}}) \cdot (\mathbf{w} \times \mathbf{n}_{|\partial K_{\ell}}).
$$

3.3. **Fortin-Soulie mixed finite elements.** We consider here the case $k = 2$, $d = 2$ and we study the so-called Fortin-Soulie mixed finite elements [3]. We consider a shape-regular triangulation sequence $(\mathcal{T}_h)_h$.

Let us consider $X_{FS}$ (resp. $X_{0,FS}$), the space of nonconforming approximation of $H^1(\Omega)$ (resp. $H_0^1(\Omega)$) of order 2:

$$
\begin{equation} \tag{3.7}
\begin{aligned}
X_{FS} &= \left\{ v_h \in P_{disc}^2(\mathcal{T}_h) ; \quad \forall f \in \mathcal{I}_F^i, \forall q_h \in P^1(F_f), \int_{F_f} [v_h]_{F_f} q_h = 0 \right\}; \\
X_{0,FS} &= \left\{ v_h \in X_{FS} ; \quad \forall f \in \mathcal{I}_F^b, \forall q_h \in P^1(F_f), \int_{F_f} v_h q_h = 0 \right\}.
\end{aligned}
\end{equation}
$$

The space of nonconforming approximation of $\mathbf{H}^1(\Omega)$ (resp. $\mathbf{H}_0^1(\Omega)$) of order 2 is
$\mathbf{X}_{\text{FS}} = (X_{\text{FS}})^2$ (resp. $\mathbf{X}_{0,\text{FS}} = (X_{0,\text{FS}})^2$). We set $\mathcal{X}_{\text{FS}} = \mathbf{X}_{0,\text{FS}} \times Q_{\text{FS}}$ where
$Q_{\text{FS}} := P_{\text{disc}}^1(\mathcal{T}_h) \cap L_{zmv}^2(\Omega)$.

The building of a basis for $X_{0,FS}$ is more involved than for $X_{0,CR}$ since we cannot
use two points per facet as degrees of freedom. Indeed, for all $\ell \in K_\ell$, it exists a
polynomial of order 2 vanishing on the Gauss-Legendre points of the facets of the
boundary $\partial K_\ell$. Let $f \in \mathcal{I}_F$. The barycentric coordinates of the two Gauss-Legendre
points $(p_+,f, p_-,f)$ on $F_f$ are such that:

$$
p_{+,f} = (c_+, c_-), p_{-,f} = (c_-, c_+), \text{ where } c_{\pm} = (1 \pm 1/\sqrt{3})/2.
$$

These points can be used to integrate exactly order three polynomials:

$$
\forall g \in P^3(F_f), \int_{F_f} g = \frac{|F_f|}{2} (g(p_{+,f}) + g(p_{-,f})).
$$
<figure><img src="image_8.png" /><figcaption>**Figure 1.** The six Gauss-Legendre points of an element $K_\ell$ and the elliptic function $\phi_{K_\ell}$.</figcaption></figure>

Notice that $\Phi_h \cap X_{LG} = \text{vec}(v_\Phi)$, where for all $\ell \in \mathcal{I}_K$, $v_{\Phi|_{K_\ell}} = \phi_{K_\ell}$. Then, counting the degrees of freedom, one can show that $\dim(X_{FS}) = \dim(X_{LG}) + \dim(\Phi_h) + 1$. For problems involving Dirichlet boundary conditions we can thus prove that for $X_{0,FS}$ the representation is unique and $X_{0,FS} = X_{0,LG} \oplus \Phi_h$. We have $X_{LG} = \text{vec}((\phi_{S_i})_{i \in \mathcal{I}_S}, (\phi_{F_f})_{f \in \mathcal{I}_F})$ where the basis functions are such that:

$$ \forall i, j \in \mathcal{I}_S, \forall f, f' \in \mathcal{I}_F: \quad \phi_{S_i}(S_j) = \delta_{i,j}, \phi_{S_i}(M_f) = 0, \quad \phi_{M_f}(M_{f'}) = \delta_{f,f'}, \phi_{M_f}(S_i) = 0. $$

For all $\ell \in \mathcal{I}_K$, we will denote by $(\phi_{\ell,j})_{j=1}^6$ the local nodal basis such that:

$$ (\phi_{\ell,j})_{j=1}^3 = (\phi_{S_i|K_\ell})_{i \in \mathcal{I}_{S,\ell}} \quad \text{and} \quad (\phi_{\ell,j})_{j=4}^6 = (\phi_{F_f|K_\ell})_{f \in \mathcal{I}_{F,\ell}}. $$

The spaces $X_{FS}$ and $X_{0,FS}$ are such that:

$$ \begin{aligned} X_{FS} &= \text{vec} \left( (\phi_{S_i})_{i \in \mathcal{I}_S}, (\phi_{F_f})_{f \in \mathcal{I}_F}, (\phi_{K_\ell})_{\ell \in \mathcal{I}_K} \right), \\ X_{0,FS} &= \text{vec} \left( (\phi_{S_i})_{i \in \mathcal{I}_S^i}, (\phi_{F_f})_{f \in \mathcal{I}_F^i}, (\phi_{K_\ell})_{\ell \in \mathcal{I}_K} \right). \end{aligned} \quad (61) $$

We propose here an alternative definition of the Fortin interpolation operator proposed in [29]. Let us first recall the Scott-Zhang interpolation operator [18, 44]. For all $i \in \mathcal{I}_S$, we choose some $\ell_i \in \mathcal{I}_{K,i}$, and we build the $L^2(K_\ell)$-dual basis $(\tilde{\phi}_{\ell,i})_{j=1}^6$ of the local nodal basis such that:

$$ \forall j, j' \in \{1, \dots, 6\}, \quad \int_{K_{\ell_i}} \tilde{\phi}_{\ell_i, j} \phi_{\ell_i, j'} = \delta_{j, j'}. $$

Let us define the Fortin-Soulie interpolation operator for scalar functions by:

$$ \pi_{FS}: \left\{ \begin{array}{l} H^1(\Omega) \to X_{FS} \\ v \mapsto \tilde{\pi}v + \sum_{\ell \in \mathcal{I}_K} v_{K_\ell} \phi_{K_\ell} \quad \text{with} \quad \tilde{\pi}v = \sum_{i \in \mathcal{I}_S} v_{S_i} \phi_{S_i} + \sum_{f \in \mathcal{I}_F} v_{F_f} \phi_{F_f}. \end{array} \right. \quad (62) $$

* The coefficients $(v_{S_i})_{i \in \mathcal{I}_S}$ are fixed so that: $\forall i \in \mathcal{I}_S$, $v_{S_i} = \int_{K_{\ell_i}} v \tilde{\phi}_{\ell_i, j_i}$, where $j_i$ is the index such that $\int_{K_{\ell_i}} \tilde{\phi}_{\ell_i, j_i} \phi_{S_i|K_{\ell_i}} = 1$.

* The coefficients $(v_{F_f})_{f \in \mathcal{I}_F}$ are fixed so that: $\forall f \in \mathcal{I}_F$, $\int_{F_f} \tilde{\pi}v = \int_{F_f} v$.

* The coefficients $v_{K_\ell}$ are fixed so that: $\int_{K_\ell} \pi_{FS}v = \int_{K_\ell} v$.

The definition (62) is more general than the one given in [29], which holds for $v \in H^2(\Omega)$.

We set $\mathbf{v}_{S_i} := (\tilde{\pi}v_1(S_i), \tilde{\pi}v_2(S_i))^T$ and $\mathbf{v}_{F_f} := (\tilde{\pi}v_1(F_f), \tilde{\pi}v_2(F_f))^T$.

We can define two different Fortin-Soulie interpolation operators for vector functions. First, let

$$ \tilde{\Pi}_{FS}: \left\{ \begin{array}{l} H^1(\Omega) \to X_{FS} \\ \mathbf{v} \mapsto \sum_{i \in \mathcal{I}_S} \mathbf{v}_{S_i} \phi_{S_i} + \sum_{f \in \mathcal{I}_F} \mathbf{v}_{F_f} \phi_{F_f} + \sum_{\ell \in \mathcal{I}_K} \tilde{\mathbf{v}}_{K_\ell} \phi_{K_\ell} \end{array} \right. $$
For all $\ell \in \mathcal{I}_K$, we define the quadratic function $\phi_{K\ell}$ that vanishes on the six Gauss-Legendre points of the facets of $K_\ell$ (see Fig. 1):

$$ (3.8) \quad \phi_{K\ell} := 2 - 3 \sum_{i \in \mathcal{I}_{S,\ell}} \lambda_{i,\ell}^2 \text{ such that } \forall f \in \mathcal{I}_{F,\ell}, \forall q \in P^1(F_f), \int_{F_f} \phi_{K\ell} q = 0. $$

<figure><img src="image_3.png" /><figcaption>FIGURE 1. The six Gauss-Legendre points of an element $K_\ell$ and the elliptic function $\phi_{K\ell}$.</figcaption></figure>

We consider the set of the elliptic functions $\phi_{K\ell}$:

$$ (3.9) \quad \Phi_h := \{ \phi_h \in L^2(\Omega); \forall \ell \in \mathcal{I}_K, \phi_{h|K_\ell} = v_{K_\ell} \phi_{K_\ell}, v_{K_\ell} \in \mathbb{R} \}. $$

We also define the spaces of $P^2$-Lagrange functions:

$$ \begin{aligned} X_{LG} &:= \{ v_h \in H^1(\Omega); \forall \ell \in \mathcal{I}_K, v_{h|K_\ell} \in P^2(K_\ell) \}, \\ X_{0, LG} &:= \{ v_h \in X_{LG}; v_{h|\partial\Omega} = 0 \}. \end{aligned} $$

The Proposition below proved in [3, Prop. 1] allows to build a basis for $X_{0,FS}$:

**Proposition 10.** We have the following decomposition: $X_{FS} = X_{LG} + \Phi_h$ with $\dim(X_{LG} \cap \Phi_h) = 1$. Any function of $X_{FS}$ can be written as the sum of a function of $X_{LG}$ and a function of $\Phi_h$. This representation can be made unique by specifying one degree of freedom.

Notice that $\Phi_h \cap X_{LG} = \text{vect}(v_\Phi)$, where for all $\ell \in \mathcal{I}_K$, $v_{\Phi|K_\ell} = \phi_{K_\ell}$. Then, counting the degrees of freedom, one can show that $\dim(X_{FS}) = \dim(X_{LG}) + \dim(\Phi_h) + 1$. For problems involving Dirichlet boundary conditions we can prove thus that for $X_{0,FS}$ the representation is unique and $X_{0,FS} = X_{0, LG} \oplus \Phi_h$. We have $X_{LG} = \text{vect}((\phi_{S_i})_{i \in \mathcal{I}_S}, (\phi_{F_f})_{f \in \mathcal{I}_F})$ where the basis functions are such that:

$$ \forall i, j \in \mathcal{I}_S, \forall f, f' \in \mathcal{I}_F \left\{ \begin{array}{l} \phi_{S_i}(S_j) = \delta_{i,j}, \quad \phi_{S_i}(M_f) = 0 \\ \phi_{M_f}(M_{f'}) = \delta_{f,f'}, \quad \phi_{M_f}(S_i) = 0 \end{array} \right. . $$

For all $\ell \in \mathcal{I}_K$, we will denote by $(\phi_{\ell,j})_{j=1}^6$ the local nodal basis such that:

$$ (\phi_{\ell,j})_{j=1}^3 = (\phi_{S_i|K_\ell})_{i \in \mathcal{I}_{S,\ell}} \quad \text{and} \quad (\phi_{\ell,j})_{j=4}^6 = (\phi_{F_f|K_\ell})_{f \in \mathcal{I}_{F,\ell}}. $$

The spaces $X_{FS}$ and $X_{0,FS}$ are such that:

$$ (3.10) \quad \left\{ \begin{array}{l} X_{FS} = \text{vect}((\phi_{S_i})_{i \in \mathcal{I}_S}, (\phi_{F_f})_{f \in \mathcal{I}_F}, (\phi_{K_\ell})_{\ell \in \mathcal{I}_K}) \\ X_{0,FS} = \text{vect}((\phi_{S_i})_{i \in \mathcal{I}_S^i}, (\phi_{F_f})_{f \in \mathcal{I}_F^i}, (\phi_{K_\ell})_{\ell \in \mathcal{I}_K}) \end{array} \right. . $$
where the coefficients $(\tilde{\mathbf{v}}_{K_\ell})_{\ell \in \mathcal{I}_K}$ are such that:

$$ \forall \ell \in \mathcal{I}_K, \quad \int_{K_\ell} \tilde{\Pi}_{FS} \mathbf{v} = \int_{K_\ell} \mathbf{v}. \qquad (63) $$

The interpolation operator $\tilde{\Pi}_{FS}$ preserves the local averages, but it doesn't preserve the divergence. We then define a second interpolation operator which preserves the divergence in a weak sense:

$$ \Pi_{FS}: \left\{ \begin{array}{l} \mathbf{H}^1(\Omega) \to \mathbf{X}_{FS} \\ \mathbf{v} \mapsto \sum_{i \in \mathcal{I}_S} \mathbf{v}_{S_i} \phi_{S_i} + \sum_{f \in \mathcal{I}_F} \mathbf{v}_{F_f} \phi_{F_f} + \sum_{\ell \in \mathcal{I}_K} \mathbf{v}_{K_\ell} \phi_{K_\ell} \end{array} \right. $$

For all $\ell \in \mathcal{I}_K$, the vector coefficient $\mathbf{v}_{K_\ell} \in \mathbb{R}^2$ is now fixed so that condition (41) is satisfied. We can impose for example that the projection $\Pi_{FS}\mathbf{v}$ satisfies:

$$ \int_{K_\ell} T_\ell^{-1}(\mathbf{x}) \operatorname{div} \Pi_{FS} \mathbf{v} = \int_{K_\ell} T_\ell^{-1}(\mathbf{x}) \operatorname{div} \mathbf{v}. \qquad (64) $$

Notice that due to (59), the patch-test condition is still satisfied.

**Proposition 24.** The Fortin-Soulie interpolation operator $\Pi_{FS}$ is such for all $0 \le s \le 1$, for all $\mathbf{v} \in \mathbf{H}^{1+s}(\Omega)$, we have:

$$ \forall \ell \in \mathcal{I}_K, \quad \| \mathrm{Grad}(\Pi_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)} \lesssim (\sigma \ell)^2 (h_\ell)^s |\mathbf{v}|_{\mathbf{H}^{1+s}(K_\ell)}, \qquad (65) $$

$$ \| \Pi_{FS} \mathbf{v} - \mathbf{v} \|_h \lesssim \sigma^2 h^s |\mathbf{v}|_{\mathbf{H}^{1+s}(\Omega)}. \qquad (66) $$

**Remark 25.** Albeit we are inspired by the proof of [22, Lemma 4], we changed the transition from Equation (4.27) to (4.29) there by using only the properties related to the normal component of the velocity, cf (56). In the original proof, one needs a stronger assumption on the regularity of $\mathbf{v}$ (namely, $\mathbf{v} \in \bigcap_{0<s<s_\Omega} \mathbf{H}^{1+s}(\Omega)$ with $s_\Omega > \frac{1}{2}$). Finally, because we do not split the integral over the boundaries of elements into the sum of $d+1$ integrals over the facets, we obtain purely local estimates, which appear to be new for the Fortin-Soulie element in the case of low-regularity fields $\mathbf{v}$.

**Proof.** Let $\mathbf{v} \in \mathbf{H}^1(\Omega)$. We have:

$$ \| \mathrm{Grad} (\Pi_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)} \le \| \mathrm{Grad} (\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) \|_{L^2(K_\ell)} + \| \mathrm{Grad} (\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)}. \qquad (67) $$

Notice that for all $\ell \in \mathcal{I}_K$, $(\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v})|_{K_\ell} = (\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}) \phi_{K_\ell}$. Using (24), we obtain that:

$$ \| \mathrm{Grad} (\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) \|_{L^2(K_\ell)} \lesssim |\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}| \| \mathrm{grad} \phi_{K_\ell} \|_{L^2(K_\ell)} \lesssim \| \mathbb{B} \ell^{-1} \| |K_\ell|^{\frac{1}{2}} |\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}|. \qquad (68) $$

Let us estimate $|\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}|$. On the one hand, we have<fnref n="2" />:

$$ \begin{align*} \int_{K_\ell} (\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) &= \int_{K_\ell} (\Pi_{FS} \mathbf{v} - \mathbf{v}) && \text{from (63),} \\ &= \int_{\partial K_\ell} \mathbf{x} (\Pi_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n} |_{\partial K_\ell} && \text{by IBP and using (64),} \\ &= \int_{\partial K_\ell} (\mathbf{x} - \bar{\mathbf{x}}) (\Pi_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n} |_{\partial K_\ell} && \text{since } \int_{\partial K_\ell} (\Pi_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n} |_{\partial K_\ell} = 0, \\ &= \int_{\partial K_\ell} (\mathbf{x} - \bar{\mathbf{x}}) (\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n} |_{\partial K_\ell} && \text{from (59).} \end{align*} $$

<footnote n="2">²Let $\mathbf{w} = \Pi_{FS}\mathbf{v} - \mathbf{v} = (w_1, w_2)^T$ and $\mathbf{x} := (x_1, x_2)^T$. By integration by parts, we have for $d' = 1, 2: \int_{\partial K_\ell} x_{d'} \mathbf{w} \cdot \mathbf{n} |_{\partial K_\ell} = \int_{K_\ell} x_{d'} \operatorname{div} \mathbf{w} + \int_{K_\ell} \mathbf{w} \cdot \mathbf{e}_{d'}$. Due to (64), $\int_{K_\ell} x_{d'} \operatorname{div} \mathbf{w} = 0$, so that for $d' = 1, 2: \int_{K_\ell} w_{d'} = \int_{\partial K_\ell} x_{d'} \mathbf{w} \cdot \mathbf{n} |_{\partial K_\ell}$.</footnote>
We propose here an alternative definition of the Fortin interpolation operator proposed in [3]. Let us first recall the Scott-Zhang interpolation operator [26, 27]. For all $i \in \mathcal{I}_S$, we choose some $\ell_i \in \mathcal{I}_{K,i}$, and we build the $L^2(K_{\ell_i})$-dual basis $(\tilde{\phi}_{\ell_i,j})_{j=1}^6$ of the local nodal basis such that:

$$ \forall j, j' \in \{1, \dots, 6\}, \quad \int_{K_{\ell_i}} \tilde{\phi}_{\ell_i, j} \phi_{\ell_i, j'} = \delta_{j, j'}. $$

Let us define the Fortin-Soulie interpolation operator for scalar functions by:

$$ (3.11) \qquad \pi_{FS} : \left\{ \begin{array}{rcl} H^1(\Omega) & \to & X_{FS} \\ v & \mapsto & \sum_{i \in \mathcal{I}_S} v_{S_i} \phi_{S_i} + \sum_{f \in \mathcal{I}_F} v_{F_f} \phi_{F_f} + \sum_{\ell \in \mathcal{I}_K} v_{K_\ell} \phi_{K_\ell} \end{array} \right. $$

* The coefficients $(v_{S_i})_{i \in \mathcal{I}_S}$ are fixed so that: $\forall i \in \mathcal{I}_S$, $v_{S_i} = \int_{K_{\ell,i}} v \tilde{\phi}_{\ell_i, j_i}$, where $j_i$ is the index such that $\int_{K_{\ell_i}} \tilde{\phi}_{\ell_i, j_i} \phi_{S_i|K_{\ell_i}} = 1$.

* The coefficients $(v_{F_f})_{f \in \mathcal{I}_F}$ are fixed so that: $\forall f \in \mathcal{I}_F$, $\int_{F_f} \tilde{\pi} v = \int_{F_f} v$.

* The coefficients $v_{K_\ell}$ are fixed so that: $\int_{K_\ell} \pi_{FS} v = \int_{K_\ell} v$.

The definition (3.11) is more general than the one given in [3], which holds for $v \in H^2(\Omega)$.

We set $v_{S_i} := (\tilde{\pi}v_1(S_i), \tilde{\pi}v_2(S_i))^T$ and $v_{F_f} := (\tilde{\pi}v_1(F_f), \tilde{\pi}v_2(F_f))^T$.

We can define two different Fortin-Soulie interpolation operators for vector functions. First, let

$$ \tilde{\Pi}_{FS} : \left\{ \begin{array}{rcl} \mathbf{H}^1(\Omega) & \to & \mathbf{X}_{FS} \\ \mathbf{v} & \mapsto & \sum_{i \in \mathcal{I}_S} \mathbf{v}_{S_i} \phi_{S_i} + \sum_{f \in \mathcal{I}_F} \mathbf{v}_{F_f} \phi_{F_f} + \sum_{\ell \in \mathcal{I}_K} \tilde{\mathbf{v}}_{K_\ell} \phi_{K_\ell} \end{array} \right. $$

where the coefficients $(\tilde{\mathbf{v}}_{K_\ell})_{\ell \in \mathcal{I}_K}$ are such that:

$$ (3.12) \qquad \forall \ell \in \mathcal{I}_K, \quad \int_{K_\ell} \tilde{\Pi}_{FS} \mathbf{v} = \int_{K_\ell} \mathbf{v}. $$

The interpolation operator $\tilde{\Pi}_{FS}$ preserves the local averages, but it doesn't preserve the divergence. We then define a second interpolation operator which preserves the divergence in a weak sense:

$$ \Pi_{FS} : \left\{ \begin{array}{rcl} \mathbf{H}^1(\Omega) & \to & \mathbf{X}_{FS} \\ \mathbf{v} & \mapsto & \sum_{i \in \mathcal{I}_S} \mathbf{v}_{S_i} \phi_{S_i} + \sum_{f \in \mathcal{I}_F} \mathbf{v}_{F_f} \phi_{F_f} + \sum_{\ell \in \mathcal{I}_K} \mathbf{v}_{K_\ell} \phi_{K_\ell} \end{array} \right. $$

For all $\ell \in \mathcal{I}_K$, the vector coefficient $\mathbf{v}_{K_\ell} \in \mathbb{R}^2$ is now fixed so that condition (2.41) is satisfied. We can impose for example that the projection $\Pi_{FS} \mathbf{v}$ satisfies:

$$ (3.13) \qquad \int_{K_\ell} T_\ell^{-1}(\mathbf{x}) \operatorname{div} \Pi_{FS} \mathbf{v} = \int_{K_\ell} T_\ell^{-1}(\mathbf{x}) \operatorname{div} \mathbf{v}. $$

Notice that due to (3.8), the patch-test condition is still satisfied.
On the other hand, it holds: $$ \int_{K_{\ell}} (\Pi_{FS}\mathbf{v} - \tilde{\Pi}_{FS}\mathbf{v}) = (\mathbf{v}_{K_{\ell}} - \tilde{\mathbf{v}}_{K_{\ell}}) \int_{K_{\ell}} \phi_{K_{\ell}} = \frac{|K_{\ell}|}{4} (\mathbf{v}_{K_{\ell}} - \tilde{\mathbf{v}}_{K_{\ell}}). $$

Hence we have:

$$ |\mathbf{v}_{K_{\ell}} - \tilde{\mathbf{v}}_{K_{\ell}}| \le 4|K_{\ell}|^{-1} \left| \int_{\partial K_{\ell}} (\mathbf{x} - \underline{\mathbf{x}}) (\tilde{\Pi}_{FS}\mathbf{v} - \mathbf{v}) \cdot \mathbf{n} |\partial K_{\ell}| \right|. \quad (69) $$

In order to bound the right-hand side of (69) component by component, we can use Lemma 21, with $q = x'_{d'} - \int_{K_{\ell}} x'_{d'} |K_{\ell}| (d' = 1, 2)$, so that $\| \mathbf{grad} q \|_{\mathbb{L}^2(K_{\ell})} = |K_{\ell}|^{\frac{1}{2}}$. We obtain:

$$ |\mathbf{v}_{K_{\ell}} - \tilde{\mathbf{v}}_{K_{\ell}}| \le 4 d \times (\sqrt{d} + 1) |K_{\ell}|^{-\frac{1}{2}} \pi^{-1} h_{\ell} \| \mathbf{Grad}(\tilde{\Pi}_{FS}\mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_{\ell})}. \quad (70) $$

Combining (68) and (70), we have:

$$ \| \mathbf{Grad} (\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) \|_{\mathbb{L}^2(K_{\ell})} \lesssim \| \mathbb{B}^{-1} \| h_{\ell} \| \mathbf{Grad} (\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_{\ell})} \lesssim \sigma_{\ell} \| \mathbf{Grad} (\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_{\ell})}. $$

For all $\mathbf{v} \in \mathbf{P}^2(K_{\ell})$ we have $\tilde{\Pi}_{FS}(\mathbf{v}) = \mathbf{v}$ and $\hat{\Pi}_{FS}\hat{\mathbf{v}}_{\ell} = \hat{\mathbf{v}}_{\ell}$. Hence, using Bramble-Hilbert/Deny-Lions Lemma [26, Lemma 11.9], we have for $m = 0, 1$, for all $\mathbf{v} \in \mathbf{H}^{1+m}(\Omega)$:

$$ \forall \ell \in \mathcal{I}_K, \quad \| \mathbf{Grad}(\tilde{\Pi}_{FS}\mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_{\ell})} \lesssim \sigma_{\ell} (h_{\ell})^m |\mathbf{v}|_{\mathbf{H}^{1+m}(K_{\ell})}. $$

We deduce that for $m = 0, 1$, for all $\mathbf{v} \in \mathbf{H}^{1+m}(\Omega)$:

$$ \begin{align} \forall \ell \in \mathcal{I}_K, \quad & \| \mathbf{Grad}(\Pi_{FS}\mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_{\ell})} \lesssim (\sigma_{\ell})^2 (h_{\ell})^m |\mathbf{v}|_{\mathbf{H}^{1+m}(K_{\ell})}, \tag{71} \\ & \| \mathbf{Grad}(\Pi_{FS}\mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(\Omega)} \lesssim \sigma_{\ell}^2 h_{\ell}^m |\mathbf{v}|_{\mathbf{H}^{1+m}(\Omega)} \quad \text{by summation.} \nonumber \end{align} $$

Using interpolation property [46, Lemma 22.2], we obtain (65) and (66). □

Hence, using the triangular inequality, we have:

$$ \| \mathbf{Grad} \Pi_{FS} \mathbf{v} \|_{\mathbb{L}^2(K_{\ell})} \le \| \mathbf{Grad} (\Pi_{FS} \mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_{\ell})} + \| \mathbf{Grad} \mathbf{v} \|_{\mathbb{L}^2(K_{\ell})} \lesssim ((\sigma_{\ell})^2 + 1) \| \mathbf{v} \|_{\mathbf{H}^{1}(K_{\ell})}. $$

By summation over $\ell$, we deduce that the coefficient $C_{nc}$ in (40) is here equal to $\sigma^2 + 1$. We recall that the discrete Poincaré–Steklov inequality (31) holds.

**Theorem 26.** Let $X_h = X_{FS}$. Then the continuous bilinear form $a_{S,h}(\cdot, \cdot)$ is $T_h$-coercive and Problem (39) is well-posed.

**Proof.** According to Proposition 24, the Fortin-Soulie interpolation operator $\Pi_{FS}$ satisfies (40)-(41), so that we can apply the proof of Theorem 16. □

Notice that in the recent paper [42], the inf-sup condition of the mixed Fortin-Soulie finite element is proved directly on a triangle and then using the macro-element technique [45], but it seems difficult to use this technique to build a Fortin operator, which is needed to compute error estimates.

The study can be extended to higher orders for $d = 2$ using the following papers: [4] for $k \ge 4$, $k$ even, [15] for $k = 3$ and [14] for $k \ge 5$, $k$ odd. In [24], the authors propose a local Fortin operator for the lowest order Taylor-Hood finite element [47] for $d = 3$.

# 4. Numerical results improving consistency

## 4.1. $\mathbf{H}(\mathrm{div})$-conforming velocity reconstruction

Consider Problem (3) with data $\mathbf{f} = -\mathbf{grad}\phi$, where $\phi \in H^1(\Omega) \cap L^2_{zmv}(\Omega)$. The unique solution is then $(\mathbf{u}, p) := (0, \phi)$. By integrating by parts, the source term in (9) reads:

$$ \forall \mathbf{v} \in \mathbf{H}_0^1(\Omega), \quad \int_{\Omega} \mathbf{f} \cdot \mathbf{v} = \int_{\Omega} \phi \, \mathrm{div} \mathbf{v}. \quad (72) $$
**Proposition 11.** The Fortin-Soulie interpolation operator $\Pi_{FS}$ is such for all $0 \le s \le 1$, for all $\mathbf{v} \in \mathbf{H}^{1+s}(\Omega)$, we have:

$$ (3.14) \quad \forall \ell \in \mathcal{I}_K, \quad \| \mathrm{Grad}(\Pi_{FS} \mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_\ell)} \lesssim (\sigma_\ell)^2 (h_\ell)^s |\mathbf{v}|_{\mathbf{H}^{1+s}(K_\ell)} $$

$$ (3.15) \quad \| \Pi_{FS} \mathbf{v} - \mathbf{v} \|_h \lesssim \sigma^2 h^s |\mathbf{v}|_{\mathbf{H}^{1+s}(\Omega)} $$

**Remark 2.** Albeit we are inspired by the proof of [2, Lemma 4], we changed the transition from equation (4.27) to (4.29) there by using only the properties related to the normal component of the velocity, cf (3.6). In the original proof, one needs a stronger assumption on the regularity of $\mathbf{v}$ (namely, $\mathbf{v} \in \bigcap_{0<s<s_\Omega} \mathbf{H}^{1+s}(\Omega)$ with $s_\Omega > \frac{1}{2}$). Finally, because we do not split the integral over the boundaries of elements into the sum of $d+1$ integrals over the facets, we obtain purely local estimates, which appear to be new for the Fortin-Soulie element in the case of low-regularity fields $\mathbf{v}$.

**Proof.** Let $\mathbf{v} \in \mathbf{H}^1(\Omega)$. We have:

$$ (3.16) \quad \| \mathrm{Grad} (\Pi_{FS} \mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_\ell)} \le \| \mathrm{Grad} (\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) \|_{\mathbb{L}^2(K_\ell)} + \| \mathrm{Grad} (\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{\mathbb{L}^2(K_\ell)} $$

Notice that for all $\ell \in \mathcal{I}_K$, $(\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v})|_{K_\ell} = (\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}) \phi_{K_\ell}$.

Using (2.23)-(i), we obtain that:

$$ (3.17) \quad \| \mathrm{Grad} (\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) \|_{\mathbb{L}^2(K_\ell)} \lesssim |\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}| \| \mathrm{grad} \, \phi_{K_\ell} \|_{\mathbb{L}^2(K_\ell)} \lesssim \| \mathbb{B}_\ell^{-1} \| |K_\ell|^{\frac{1}{2}} |\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}|. $$

Let us estimate $|\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}|$. On the one hand, we have<fnref n="3" />:

$$ \begin{align*} \int_{K_\ell} (\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) &= \int_{K_\ell} (\Pi_{FS} \mathbf{v} - \mathbf{v}) \text{ from (3.12)}, \\ &= \int_{\partial K_\ell} \mathbf{x} (\Pi_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n}_{|\partial K_\ell} \text{ by IBP and using (3.13)}, \\ &= \int_{\partial K_\ell} (\mathbf{x} - \underline{\mathbf{x}}) (\Pi_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n}_{|\partial K_\ell} \\ &\qquad \text{since } \int_{\partial K_\ell} (\Pi_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n}_{|\partial K_\ell} = 0, \\ &= \int_{\partial K_\ell} (\mathbf{x} - \underline{\mathbf{x}}) (\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n}_{|\partial K_\ell} \text{ from (3.8)}. \end{align*} $$

On the other hand, it holds:

$$ \int_{K_\ell} (\Pi_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) = (\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}) \int_{K_\ell} \phi_{K_\ell} = \frac{|K_\ell|}{4} (\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}). $$

Hence we have:

$$ (3.18) \quad |\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}| \le 4 |K_\ell|^{-1} \left| \int_{\partial K_\ell} (\mathbf{x} - \underline{\mathbf{x}}) (\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \cdot \mathbf{n}_{|\partial K_\ell} \right|. $$

<footnote n="3">³Let set $\mathbf{w} = \Pi_{FS} \mathbf{v} - \mathbf{v} = (w_1, w_2)^T$ and $\mathbf{x} := (x_1, x_2)^T$. By IBP, we have for $d' = 1, 2$: $\int_{\partial K_\ell} x_{d'} \mathbf{w} \cdot \mathbf{n}_{|\partial K_\ell} = \int_{K_\ell} x_{d'} \mathrm{div} \mathbf{w} + \int_{K_\ell} \mathbf{w} \cdot \mathbf{e}_{d'}$. Due to (3.13), $\int_{K_\ell} x_{d'} \mathrm{div} \mathbf{w} = 0$, so that for $d' = 1, 2$: $\int_{K_\ell} w_{d'} = \int_{\partial K_\ell} x_{d'} \mathbf{w} \cdot \mathbf{n}_{|\partial K_\ell}$.</footnote>
Recall that the nonconforming space $\mathbf{X}_h$ defined in (30) is a subset of $\mathcal{P}_h\mathbf{H}^1$: using a nonconforming finite element method, the integration by parts must be done on each element of the triangulation, and we have:

$$ \forall \mathbf{v}_h \in \mathbf{V}_h, \quad \int_{\Omega} \mathbf{f} \cdot \mathbf{v}_h = \sum_{f \in \mathcal{F}_F} \int_{F_f} [\mathbf{v}_h - \mathbf{v}_{h, \ell}]_{F_f} \cdot \mathbf{n}_f \phi. \quad (73) $$

Using Corollary 22, we have: $\sum_{f \in \mathcal{F}_F} \int_{F_f} [\mathbf{v}_h - \mathbf{v}_{h, \ell}]_{F_f} \cdot \mathbf{n}_f \phi \lesssim h \|\mathbf{v}_h\|_h \|\mathbf{grad} \phi\|_{L^2(\Omega)}$. Applying (73) to the right-hand side of (39) and choosing $\mathbf{v}_h = \mathbf{u}_h$, it holds: $v \|\mathbf{u}_h\|_h \lesssim h \|\mathbf{grad} \phi\|_{L^2(\Omega)}$ (as expected by (47)). Hence, the term with the jumps acts as a numerical source for the discrete velocity, whose numerical influence is proportional to $h/v$. Thus, we cannot obtain exactly $\mathbf{u}_h = 0$. Linke proposed in [38] to project the test function $\mathbf{v}_h \in \mathbf{X}_h$ on a discrete subspace of $\mathbf{H}(\text{div}; \Omega)$, like Raviart-Thomas or Brezzi-Douglas-Marini finite elements (see [12,39], or the monograph [8]). Let $\Pi_{\text{div}}: \mathbf{X}_{0,h} \to \mathbf{P}_{\text{disc}}^k(\mathcal{T}_h) \cap \mathbf{H}_0(\text{div}; \Omega)$ be some interpolation operator built so that for all $\mathbf{v}_h \in \mathbf{X}_{0,h}$, for all $\ell \in \mathcal{I}_K$, $(\text{div}\Pi_{\text{div}}\mathbf{v}_h)|_{K_\ell} = \text{div}\mathbf{v}_h|_{K_\ell}$. Integrating by parts, we have for all $\mathbf{v}_h \in \mathbf{X}_{0,h}$:

$$ \int_{\Omega} \mathbf{f} \cdot \Pi_{\text{div}} \mathbf{v}_h = \int_{\Omega} \phi \, \text{div} \Pi_{\text{div}} \mathbf{v}_h = \sum_{\ell \in K_\ell} \int_{K_\ell} \phi \, \text{div} \Pi_{\text{div}} \mathbf{v}_h = \sum_{\ell \in K_\ell} \int_{K_\ell} \phi \, \text{div} \mathbf{v}_h = (\text{div}_h \mathbf{v}_h, \phi)_{L^2(\Omega)}. $$

The projection $\Pi_{\text{div}}$ allows to eliminate the terms of the integrals of the jumps in (73). Let us write Problem (39) as:

$$ \text{Find } (\mathbf{u}_h, p_h) \in \mathcal{X}_h \text{ such that } a_{S,h}((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) = \ell_f((\Pi_{\text{div}}\mathbf{v}_h, q_h)) \quad \forall (\mathbf{v}_h, q_h) \in \mathcal{X}_h. \quad (74) $$

In the case of $\mathcal{X}_h = \mathcal{X}_{CR}$ and a projection on Brezzi-Douglas-Marini finite elements, the following error estimate holds if $(\mathbf{u}, p) \in \mathbf{H}^2(\Omega) \times \mathbf{H}^1(\Omega)$:

$$ \|\mathbf{u} - \mathbf{u}_h\|_{L^2(\Omega)} \le \tilde{C} h^2 |\mathbf{u}|_{\mathbf{H}^2(\Omega)}, \quad (75) $$

where the constant $\tilde{C}$ is independent of $h$. The proof is detailed in [10] for shape-regular meshes and [2] for anisotropic meshes. We remark that the error doesn't depend on the norm of the pressure nor on the $v$ parameter. We will provide some numerical results to illustrate the effectiveness of this formulation, even with a projection on the Raviart-Thomas finite elements, which, for a fixed polynomial order, are less precise than the Brezzi-Douglas-Marini finite elements. For $k \in \mathbb{N}^*$, the space of Raviart-Thomas finite elements can be defined as:

$$ \mathbf{X}_{RT_k} := \{ \mathbf{v} \in \mathbf{H}(\text{div}; \Omega) | \forall \ell \in \mathcal{I}_k, \mathbf{v}|_{K_\ell} = \mathbf{a}_\ell + b_\ell \mathbf{x}, (\mathbf{a}_\ell, b_\ell) \in \mathbf{P}^k(K_\ell) \times \mathbf{P}^k(K_\ell) \}. $$

Let $k \le 1$. The Raviart-Thomas interpolation operator $\Pi_{RT_k}: \mathbf{H}^1(\Omega) \cup \mathbf{X}_h \to \mathbf{X}_{RT_k}$ is defined by:

$$ \forall \mathbf{v} \in \mathbf{H}^1(\Omega) \cup \mathbf{X}_h, \quad \left\{ \begin{array}{ll} \forall f \in \mathcal{F}_F, & \displaystyle \int_{F_f} \Pi_{RT_k} \mathbf{v} \cdot \mathbf{n}_f q = \int_{F_f} \mathbf{v} \cdot \mathbf{n}_f q, & \forall q \in \mathbf{P}^k(F_f) \\ \\ \text{for } k=1, \forall \ell \in \mathcal{I}_K, & \displaystyle \int_{K_\ell} \Pi_{RT_1} \mathbf{v} = \int_{K_\ell} \mathbf{v} & . \end{array} \right. \quad (76) $$

Note that the Raviart-Thomas interpolation operator preserves the constants. Let $\mathbf{v}_h \in \mathbf{X}_h$. In order to compute the left-hand side of (73), we must evaluate $(\Pi_{RT_k} \mathbf{v}_h)|_{K_\ell}$ for all $\ell \in \mathcal{I}_K$. Calculations are performed using the proposition below:

**Proposition 27 ( [31, Lemma 3.11] ).** Let $k \le 1$. Let $\hat{\Pi}_{RT_k}: \mathbf{H}^1(\hat{K}) \to \mathbf{P}^k(\hat{K})$ be the Raviart-Thomas interpolation operator restricted to the reference element, so that:

$$ \forall \hat{\mathbf{v}} \in \mathbf{H}^1(\hat{K}), \quad \left\{ \begin{array}{ll} \forall \hat{F} \in \partial \hat{K}, & \displaystyle \int_{\hat{F}} \hat{\Pi}_{RT_k} \hat{\mathbf{v}} \cdot \mathbf{n}_{\hat{F}} \hat{q} = \int_{\hat{F}} \hat{\mathbf{v}} \cdot \mathbf{n}_{\hat{F}} \hat{q}, & \forall \hat{q} \in \mathbf{P}^k(\hat{F}) \\ \\ \text{for } k=1, & \displaystyle \int_{\hat{K}} \hat{\Pi}_{RT_k} \hat{\mathbf{v}} = \int_{\hat{K}} \hat{\mathbf{v}} & . \end{array} \right. \quad (77) $$
In order to bound the right-hand-side of (3.18) component by component, we can use Lemma 4, with $q = x_d' - \int_{K_\ell} x_{d'}/|K_\ell|$ ($d' = 1, 2$), so that $||\text{grad } q||_{L^2(K_\ell)} = |K_\ell|^{1/2}$. We obtain:

$$ (3.19) \quad |\mathbf{v}_{K_\ell} - \tilde{\mathbf{v}}_{K_\ell}| \le 4d \times (\sqrt{d} + 1) |K_\ell|^{-1/2} \pi^{-1} h_\ell \| \text{Grad}(\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)}. $$

Combining (3.17) and (3.19), we have:

$$ \begin{aligned} \| \text{Grad}(\tilde{\Pi}_{FS} \mathbf{v} - \tilde{\Pi}_{FS} \mathbf{v}) \|_{L^2(K_\ell)} &\lesssim \| \mathbb{B}_\ell^{-1} \| h_\ell \| \text{Grad}(\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)}, \\ &\lesssim \sigma_\ell \| \text{Grad}(\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)}. \end{aligned} $$

For all $\mathbf{v} \in \mathbf{P}^2(K_\ell)$ we have $\tilde{\Pi}_{FS}(\mathbf{v}) = \mathbf{v}$ and $\tilde{\Pi}_{FS}\hat{\mathbf{v}}_\ell = \hat{\mathbf{v}}_\ell$. Hence, using Bramble-Hilbert/Deny-Lions Lemma [18, Lemma 11.9], we have for $m = 0, 1$, for all $\mathbf{v} \in \mathbf{H}^{1+m}(\Omega)$:

$$ \forall \ell \in \mathcal{I}_K, \quad \| \text{Grad}(\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)} \lesssim \sigma_\ell (h_\ell)^m |\mathbf{v}|_{\mathbf{H}^{1+m}(K_\ell)}. $$

We deduce that for $m = 0, 1$, for all $\mathbf{v} \in \mathbf{H}^{1+m}(\Omega)$, for all $\ell \in \mathcal{I}_K$:

$$ \| \text{Grad}(\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)} \lesssim (\sigma_\ell)^2 (h_\ell)^m |\mathbf{v}|_{\mathbf{H}^{1+m}(K_\ell)}. $$

Hence, by summation, we get that for $m = 0, 1$, for all $\mathbf{v} \in \mathbf{H}^{1+m}(\Omega)$:

$$ \| \text{Grad}(\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(\Omega)} \lesssim \sigma^2 h^m |\mathbf{v}|_{\mathbf{H}^{1+m}(\Omega)}. $$

Using interpolation property [28, Lemma 22.2], we obtain (3.14) and (3.15). $\square$

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

Hence, using the triangular inequality, we have:

$$ \begin{aligned} \| \text{Grad} \tilde{\Pi}_{FS} \mathbf{v} \|_{L^2(K_\ell)} &\le \| \text{Grad}(\tilde{\Pi}_{FS} \mathbf{v} - \mathbf{v}) \|_{L^2(K_\ell)} + \| \text{Grad} \mathbf{v} \|_{L^2(K_\ell)}, \\ &\lesssim ((\sigma_\ell)^2 + 1) \| \mathbf{v} \|_{\mathbf{H}^{1+m}(K_\ell)}. \end{aligned} $$

By summation over $\ell$, we deduce that the coefficient $C_{nc}$ in (2.40) is here equal to $\sigma^2 + 1$. We recall that the discrete Poincaré-Steklov inequality (2.30) holds.

**Theorem 6.** Let $\mathcal{X}_h = \mathcal{X}_{FS}$. Then the continuous bilinear form $a_{S,h}(\cdot, \cdot)$ is $T_h$-coercive and Problem (2.39) is well-posed.

*Proof.* According to Proposition 11, the Fortin-Soulie interpolation operator $\Pi_{FS}$ satisfies (2.40)-(2.41), so that we can apply the proof of Theorem 4. $\square$

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

Notice that in the recent paper [29], the inf-sup condition of the mixed Fortin-Soulie finite element is proven directly on a triangle and then using the macro-element technique [30], but it seems difficult to use this technique to build a Fortin operator, which is needed to compute error estimates.

The study can be extended to higher orders for $d = 2$ using the following papers: [31] for $k \ge 4$, $k$ even, [32] for $k = 3$ and [25] for $k \ge 5$, $k$ odd. In [33], the authors propose a local Fortin operator for the lowest order Taylor-Hood finite element [16] for $d = 3$.

## 4. NUMERICAL RESULTS IMPROVING CONSISTENCY

4.1. **H(div)-conforming velocity reconstruction.** Consider Problem (2.3) with data $\mathbf{f} = -\text{grad} \phi$, where $\phi \in H^1(\Omega) \cap L^2_{zmv}(\Omega)$. The unique solution is then $(\mathbf{u}, p) := (0, \phi)$. By integrating by parts, the source term in (2.8) reads:

$$ (4.1) \qquad \forall \mathbf{v} \in \mathbf{H}_0^1(\Omega), \quad \int_{\Omega} \mathbf{f} \cdot \mathbf{v} = \int_{\Omega} \phi \, \mathrm{div} \, \mathbf{v}. $$
We then have: $\forall \ell \in \mathcal{I}_K$,

$$ (\Pi_{RT_k} \mathbf{v})|_{K_\ell}(\mathbf{x}) = \mathbb{B}_\ell (\hat{\Pi}_{RT_k} \mathbb{B}_\ell^{-1} \hat{\mathbf{v}}_\ell) \circ T_\ell^{-1}(\mathbf{x}) \quad \text{where} \quad \hat{\mathbf{v}}_\ell = \mathbf{v} \circ T_\ell(\hat{\mathbf{x}}). \qquad (78) $$

The proof is based on the equality of the $\hat{F}$ and $\hat{K}$-moments of $(\Pi_{RT_k} \mathbf{v})|_{K_\ell} \circ T_\ell(\hat{\mathbf{x}})$ and $\mathbb{B}_\ell(\hat{\Pi}_{RT_k} \mathbb{B}_\ell^{-1} \hat{\mathbf{v}}_\ell)(\hat{\mathbf{x}})$. For $k=0$, setting for $d' \in \{1, \dots, d\}$: $\psi_{f,d'} := \psi_f e_{d'}$, we obtain that:

$$ \forall \ell \in \mathcal{I}_K, \forall f \in \mathcal{I}_{F,\ell}, \quad (\Pi_{RT_0} \psi_{f,d'})|_{K_\ell} = (d|K_\ell|)^{-1} (\mathbf{x} - \bar{O}S_{f,\ell}) \mathcal{I}_{f,\ell} \cdot \mathbf{e}_{d'}, \qquad (79) $$

where $S_{f,\ell}$ is the vertex opposite to $F_f$ in $K_\ell$.

For $k=1$, the vector $(\Pi_{RT_1} \mathbf{v}_h)|_{K_\ell}$ is described by eight unknowns:

$$ (\Pi_{RT_1} \mathbf{v}_h)|_{K_\ell} = \mathbb{A}_\ell \mathbf{x} + (\mathbf{b}_\ell \cdot \mathbf{x}) \mathbf{x} + \mathbf{d}_\ell, \quad \text{where} \quad \mathbb{A}_\ell \in \mathbb{R}^{2 \times 2}, \mathbf{b}_\ell \in \mathbb{R}^2, \mathbf{d}_\ell \in \mathbb{R}^2. $$

We compute only once the inverse of the matrix of the linear system (77), in $\mathbb{R}^{8 \times 8}$.

## 4.2. Application with manufactured solutions

In the Tables 1, 2 and 3, we call $\varepsilon_0^\nu(\mathbf{u}) = \|\mathbf{u} - \mathbf{u}_h\|_{L^2(\Omega)}/\|\mathbf{u}, \mathbf{p}\|_{\mathcal{X}}$ the velocity error in $L^2(\Omega)$-norm, where $\mathbf{u}_h$ is the solution to Problem (39) (columns $\mathbf{X}_{CR}$ and $\mathbf{X}_{FS}$) or (74) (columns $\mathbf{X}_{CR} + \Pi_{RT_0}$ and $\mathbf{X}_{FS} + \Pi_{RT_1}$) and $\mathbf{h}$ is the mesh size.

We first consider the Stokes problem (3) in $\Omega = (0, 1)^2$ with $\mathbf{u} = 0$, $\mathbf{p} = (x_1)^3 + (x_2)^3 - 0.5$, $\mathbf{f} = \mathbf{grad} \mathbf{p} = 3((x_1)^2, (x_2)^2)^T$. We report in Table 1 the velocity error $\varepsilon_0^\nu(\mathbf{u}) := \nu \|\mathbf{u}_h\|_{L^2(\Omega)} / \|\mathbf{p}\|_{L^2(\Omega)}$, for $h = 5.00 \times 10^{-2}$ and for different values of $\mathbf{v}$.

<table><caption>**Table 1.** Values of $\varepsilon_0^\nu(\mathbf{u})$ for $h = 5.00 \times 10^{-2}$.</caption><thead><tr><th>ν</th><th>X<sub>CR</sub></th><th>X<sub>CR</sub> + Π<sub>RT<sub>0</sub></sub></th><th>X<sub>FS</sub></th><th>X<sub>FS</sub> + Π<sub>RT<sub>1</sub></sub></th></tr></thead><tbody><tr><td>1.00 × 10<sup>-4</sup></td><td>7.96 × 10<sup>-4</sup></td><td>4.59 × 10<sup>-17</sup></td><td>8.81 × 10<sup>-7</sup></td><td>1.54 × 10<sup>-16</sup></td></tr><tr><td>1.00 × 10<sup>-5</sup></td><td>7.96 × 10<sup>-4</sup></td><td>4.59 × 10<sup>-17</sup></td><td>8.81 × 10<sup>-7</sup></td><td>1.54 × 10<sup>-16</sup></td></tr><tr><td>1.00 × 10<sup>-6</sup></td><td>7.96 × 10<sup>-4</sup></td><td>4.59 × 10<sup>-17</sup></td><td>8.81 × 10<sup>-7</sup></td><td>1.54 × 10<sup>-16</sup></td></tr></tbody></table>

The $L^2(\Omega)$-norm of the discrete velocity $\|\mathbf{u}_h\|_{L^2(\Omega)}$ is proportional to $\nu^{-1}$. Using the projection, we obtain $\varepsilon_0^\nu(\mathbf{u}) = 0$ close to machine precision.

We now consider the Stokes problem (3) in $\Omega = (0, 1)^2$ with:

$$ \mathbf{u} = \begin{pmatrix} (1 - \cos(2\pi x_1)) \sin(2\pi x_2) \\ (\cos(2\pi x_2) - 1) \sin(2\pi x_1) \end{pmatrix}, \quad \mathbf{p} = \sin(2\pi x_1) \text{ and } \mathbf{f} = -\nu \Delta \mathbf{u} + \mathbf{grad} \mathbf{p}. $$

We report in Table 2 (resp. 3) the values of $\varepsilon_0^\nu(\mathbf{u})$ in the case $\nu = 1.00 \times 10^{-3}$ (resp. $\nu = 1.00 \times 10^{-4}$) for different mesh sizes. We observe that when there is no projection, $\varepsilon_0^\nu(\mathbf{u})$ is independent of $\mathbf{v}$, whereas using the projection, $\varepsilon_0^\nu(\mathbf{u})$ is proportional to $\mathbf{v}$.

<table><caption>**Table 2.** Values of $\varepsilon_0^\nu(\mathbf{u})$ for $\nu = 1.00 \times 10^{-3}$.</caption><thead><tr><th>h</th><th>X<sub>CR</sub></th><th>X<sub>CR</sub> + Π<sub>RT<sub>0</sub></sub></th><th>X<sub>FS</sub></th><th>X<sub>FS</sub> + Π<sub>RT<sub>1</sub></sub></th></tr></thead><tbody><tr><td>5.00 × 10<sup>-2</sup></td><td>1.32 × 10<sup>-3</sup></td><td>2.74 × 10<sup>-5</sup></td><td>4.73 × 10<sup>-6</sup></td><td>5.05 × 10<sup>-7</sup></td></tr><tr><td>2.50 × 10<sup>-2</sup></td><td>3.30 × 10<sup>-4</sup></td><td>6.93 × 10<sup>-6</sup></td><td>5.06 × 10<sup>-7</sup></td><td>6.42 × 10<sup>-8</sup></td></tr><tr><td>1.25 × 10<sup>-2</sup></td><td>8.25 × 10<sup>-5</sup></td><td>1.74 × 10<sup>-6</sup></td><td>6.31 × 10<sup>-8</sup></td><td>8.10 × 10<sup>-9</sup></td></tr><tr><td>6.25 × 10<sup>-3</sup></td><td>2.04 × 10<sup>-5</sup></td><td>4.35 × 10<sup>-7</sup></td><td>7.44 × 10<sup>-9</sup></td><td>1.03 × 10<sup>-9</sup></td></tr><tr><td>Rate</td><td>h<sup>2.00</sup></td><td>h<sup>1.99</sup></td><td>h<sup>3.08</sup></td><td>h<sup>2.97</sup></td></tr></tbody></table>
Recall that the nonconforming space $\mathbf{X}_h$ defined in (2.28) is a subset of $\mathcal{P}_h\mathbf{H}^1$: using a nonconforming finite element method, the integration by parts must be done on each element of the triangulation, and we have:

$$ (4.2) \quad \forall \mathbf{v} \in \mathcal{P}_h \mathbf{H}^1, \quad \int_{\Omega} \mathbf{f} \cdot \mathbf{v} = (\mathrm{div}_h \mathbf{v}, \phi)_{L^2(\Omega)} + \sum_{f \in \mathcal{I}_F} \int_{F_f} [\mathbf{v}]_{F_f} \cdot \mathbf{n}_f \phi. $$

Using Lemma 4, we have: $\sum_{f \in \mathcal{I}_F} \int_{F_f} [\mathbf{v}]_{F_f} \cdot \mathbf{n}_f \phi \lesssim h \|\mathbf{v}_h\|_h \|\mathrm{grad} \phi\|_{L^2(\Omega)}$. Applying (4.2) to the right-hand-side of (2.39) and choosing $\mathbf{v}_h = \mathbf{u}_h$, it holds: $\nu \|\mathbf{u}_h\|_h \lesssim h \|\mathrm{grad} \phi\|_{L^2(\Omega)}$ (as expected by (2.47)). Hence, the term with the jumps acts as a numerical source for the discrete velocity, which numerical influence is proportional to $h/\nu$. Thus, we cannot obtain exactly $\mathbf{u}_h = 0$. Linke proposed in [34] to project the test function $\mathbf{v}_h \in \mathbf{X}_h$ on a discrete subspace of $\mathbf{H}(\mathrm{div}; \Omega)$, like Raviart-Thomas or Brezzi-Douglas-Marini finite elements (see [35, 36], or the monograph [17]). Let $\Pi_{\mathrm{div}}: \mathbf{X}_{0,h} \to P_{\mathrm{disc}}^k(\mathcal{T}_h) \cap \mathbf{H}_0(\mathrm{div}; \Omega)$ be some interpolation operator built so that for all $\mathbf{v}_h \in \mathbf{X}_{0,h}$, for all $\ell \in \mathcal{I}_K$, $(\mathrm{div} \Pi_{\mathrm{div}} \mathbf{v}_h)|_{K_\ell} = \mathrm{div} \mathbf{v}_{h|K_\ell}$. Integrating by parts, we have for all $\mathbf{v}_h \in \mathbf{X}_{0,h}$:

$$ \begin{aligned} \int_{\Omega} \mathbf{f} \cdot \Pi_{\mathrm{div}} \mathbf{v}_h &= \int_{\Omega} \phi \mathrm{div} \Pi_{\mathrm{div}} \mathbf{v}_h = \sum_{\ell \in K_{\ell}} \int_{K_{\ell}} \phi \mathrm{div} \Pi_{\mathrm{div}} \mathbf{v}_h, \\ &= \sum_{\ell \in K_{\ell}} \int_{K_{\ell}} \phi \mathrm{div} \mathbf{v}_h = (\mathrm{div}_h \mathbf{v}_h, \phi)_{L^2(\Omega)}. \end{aligned} $$

The projection $\Pi_{\mathrm{div}}$ allows to eliminate in Equation (4.2) the last term of the right-hand-side.

Let us write Problem (2.39) as: Find $(\mathbf{u}_h, p_h) \in \mathcal{X}_h$ such that

$$ (4.3) \quad a_{S,h}((\mathbf{u}_h, p_h), (\mathbf{v}_h, q_h)) = \ell_{\mathbf{f}} ((\Pi_{\mathrm{div}}\mathbf{v}_h, q_h)) \quad \forall (\mathbf{v}_h, q_h) \in \mathcal{X}_h. $$

In the case of $\mathcal{X}_h = \mathcal{X}_{CR}$ and a projection on Brezzi-Douglas-Marini finite elements, the following error estimate holds if $(\mathbf{u}, p) \in \mathbf{H}^2(\Omega) \times \mathbf{H}^1(\Omega)$:

$$ (4.4) \quad \|\mathbf{u} - \mathbf{u}_h\|_{\mathbf{L}^2(\Omega)} \le \tilde{C} h^2 |\mathbf{u}|_{\mathbf{H}^2(\Omega)}, $$

where the constant $\tilde{C}$ if independent of $h$. The proof is detailed in [37] for shape-regular meshes and [38] for anisotropic meshes. We remark that the error doesn't depend on the norm of the pressure nor on the $\nu$ parameter. We will provide some numerical results to illustrate the effectiveness of this formulation, even with a projection on the Raviart-Thomas finite elements, which, for a fixed polynomial order, are less precise than the Brezzi-Douglas-Marini finite elements.

For all $\ell \in \mathcal{I}_K$, we let $P_H^k(K_\ell)$ be the set of homogeneous polynomials of order $k$ on $K_\ell$.

For $k \in \mathbb{N}^*$, the space of Raviart-Thomas finite elements can be defined as:

$$ \mathbf{X}_{RT_k} := \{ \mathbf{v} \in \mathbf{H}(\mathrm{div}; \Omega); \forall \ell \in \mathcal{I}_k, \mathbf{v}|_{K_\ell} = \mathbf{a}_\ell + b_\ell \mathbf{x} \mid (\mathbf{a}_\ell, b_\ell) \in P^k(K_\ell)^d \times P_H^k(K_\ell) \}. $$

Let $k \le 1$. The Raviart-Thomas interpolation operator $\Pi_{RT_k}: \mathbf{H}^1(\Omega) \cup \mathbf{X}_h \to \mathbf{X}_{RT_k}$ is defined by: $\forall \mathbf{v} \in \mathbf{H}^1(\Omega) \cup \mathbf{X}_h$,

$$ (4.5) \quad \left\{ \begin{array}{ll} \forall f \in \mathcal{I}_F, & \displaystyle \int_{F_f} \Pi_{RT_k} \mathbf{v} \cdot \mathbf{n}_f q = \int_{F_f} \mathbf{v} \cdot \mathbf{n}_f q, & \forall q \in P^k(F_f) \\ \text{for } k=1, \forall \ell \in \mathcal{I}_K, & \displaystyle \int_{K_\ell} \Pi_{RT_1} \mathbf{v} = \int_{K_\ell} \mathbf{v} & \end{array} \right. $$
<table><caption>**Table 3.** Values of $\epsilon_0^\nu(\mathbf{u})$ for $\nu = 1.00 \times 10^{-4}$.</caption><thead><tr><th>$h$</th><th>$\mathbf{X}_{CR}$</th><th>$\mathbf{X}_{CR} + \Pi_{RT_0}$</th><th>$\mathbf{X}_{FS}$</th><th>$\mathbf{X}_{FS} + \Pi_{RT_1}$</th></tr></thead><tbody><tr><td>$5.00 \times 10^{-2}$</td><td>$1.32 \times 10^{-3}$</td><td>$2.74 \times 10^{-6}$</td><td>$4.70 \times 10^{-6}$</td><td>$5.05 \times 10^{-8}$</td></tr><tr><td>$2.50 \times 10^{-2}$</td><td>$3.30 \times 10^{-4}$</td><td>$6.93 \times 10^{-7}$</td><td>$5.10 \times 10^{-7}$</td><td>$6.43 \times 10^{-9}$</td></tr><tr><td>$1.25 \times 10^{-2}$</td><td>$8.25 \times 10^{-5}$</td><td>$1.74 \times 10^{-7}$</td><td>$6.37 \times 10^{-8}$</td><td>$8.11 \times 10^{-10}$</td></tr><tr><td>$6.25 \times 10^{-3}$</td><td>$2.04 \times 10^{-5}$</td><td>$4.36 \times 10^{-8}$</td><td>$7.51 \times 10^{-9}$</td><td>$9.77 \times 10^{-11}$</td></tr><tr><td>Rate</td><td>$h^{2.00}$</td><td>$h^{1.99}$</td><td>$h^{3.08}$</td><td>$h^{2.99}$</td></tr></tbody></table>

Let us consider the Stokes problem (3) with a low-regular velocity. Let $\Omega = (0, 1)^2$, $S_0 = (0.5, 0.5)$, and $(r, \theta)$ be the polar coordinates centered on $S_0$. We set:

$$ \mathbf{u} = r^{\alpha} \mathbf{e}_{\theta}, \quad p = r - |\Omega|^{-1} \int_{\Omega} r \text{ so that } \mathbf{f} := -\nabla \Delta \mathbf{u} + \mathbf{grad} \, p = \nu(1 - \alpha^2) r^{\alpha-2} \mathbf{e}_{\theta} + \mathbf{e}_r. $$

We report in Table 4 the values of $\epsilon_0^\nu(\mathbf{u})$ for $\nu = 1.00 \times 10^{-4}$, and for different mesh sizes, with $\alpha = 1$ and $\alpha = 0.45$. For $\alpha = 1$, $\mathbf{u} = (-y, x)^T \in \mathbf{H}^2(\Omega)$. For $\alpha = 0.45$, $\mathbf{u} \in \bigcap_{0<s<\alpha} \mathbf{H}^{1+s}(\Omega)$, hence $\mathbf{u} \notin \mathbf{H}^2(\Omega)$. It seems that the Raviart-Thomas projection is less efficient in that last case.

<table><caption>**Table 4.** Values of $\epsilon_0^\nu(\mathbf{u})$, regular and low-regular velocity, $\nu = 1.00 \times 10^{-4}$.</caption><thead><tr><th></th><th colspan="2">α = 1</th><th colspan="2">α = 0.45</th></tr><tr><th>$h$</th><th>$\mathbf{X}_{FS}$</th><th>$\mathbf{X}_{FS} + \Pi_{RT_1}$</th><th>$\mathbf{X}_{FS}$</th><th>$\mathbf{X}_{FS} + \Pi_{RT_1}$</th></tr></thead><tbody><tr><td>$1.00 \times 10^{-1}$</td><td>$3.03 \times 10^{-5}$</td><td>$2.81 \times 10^{-6}$</td><td>$3.05 \times 10^{-5}$</td><td>$3.94 \times 10^{-6}$</td></tr><tr><td>$5.00 \times 10^{-2}$</td><td>$4.34 \times 10^{-6}$</td><td>$1.54 \times 10^{-6}$</td><td>$4.57 \times 10^{-6}$</td><td>$2.15 \times 10^{-6}$</td></tr><tr><td>$2.50 \times 10^{-2}$</td><td>$4.72 \times 10^{-7}$</td><td>$2.41 \times 10^{-8}$</td><td>$9.70 \times 10^{-7}$</td><td>$8.52 \times 10^{-7}$</td></tr><tr><td>Rate</td><td>$h^{3.00}$</td><td>$h^{3.43}$</td><td>$h^{2.48}$</td><td>$h^{1.11}$</td></tr></tbody></table>

In order to enhance the numerical results, one could also use a posteriori error estimators to adapt the mesh near point $S_0$ (see [23, 25] for $k = 1$ and [1] for $k = 2$).

Alternatively, using the nonconforming Crouzeix-Raviart mixed finite element method, one can build a divergence-free basis, as described in [34] for $k = 1$. When $k = 1$, following the initial work of [7], one can also add $P^1$-Lagrange basis functions to the space of the discrete pressures as explained in [37]. The consistency of the discrete velocity is then improved. Notice that using conforming finite elements, the Scott-Vogelius finite elements [28, 43, 49] produce velocity approximations that are exactly divergence free.

The code used to get the numerical results can be downloaded on GitHub [35].

## 5. Conclusion

We analysed the discretization of Stokes problem with nonconforming finite elements in light of the T-coercivity theory. Furthermore, we obtained local stability estimates for order 1 in 2 or 3 dimensions without mesh regularity assumption; and for order 2 in 2 dimensions in the case of a shape-regular triangulation sequence. This local approach, splitting the normal and the tangential components could help to generalize our results to order $k \ge 3$ (using maybe also other internal moment conservation). This is ongoing work. We then provided numerical results to illustrate the importance of using $\mathbf{H}(\text{div})$-conforming projection. Further, we intend to extend the study to other mixed finite element methods.
Note that the Raviart-Thomas interpolation operator preserves the constants. Let $\mathbf{v}_h \in \mathbf{X}_h$. In order to compute the left-hand-side of (4.2), we must evaluate $(\Pi_{RT_k} \mathbf{v}_h)|_{K_\ell}$ for all $\ell \in \mathcal{I}_K$. Calculations are performed using the proposition below, which corresponds to [39, Lemma 3.11]:

**Proposition 12.** Let $k \le 1$. Let $\hat{\Pi}_{RT_k} : \mathbf{H}^1(\hat{K}) \to \mathbf{P}^k(\hat{K})$ be the Raviart-Thomas interpolation operator restricted to the reference element, so that: $\forall \hat{\mathbf{v}} \in \mathbf{H}^1(\hat{K})$,

$$ (4.6) \qquad \left\{ \begin{array}{ll} \forall \hat{F} \in \partial \hat{K}, & \displaystyle \int_{\hat{F}} \hat{\Pi}_{RT_k} \hat{\mathbf{v}} \cdot \mathbf{n}_{\hat{F}} \hat{q} = \int_{\hat{F}} \hat{\mathbf{v}} \cdot \mathbf{n}_{\hat{F}} \hat{q}, & \forall \hat{q} \in P^k(\hat{F}) \\ \text{for } k=1, & \displaystyle \int_{\hat{K}} \hat{\Pi}_{RT_k} \hat{\mathbf{v}} = \int_{\hat{K}} \hat{\mathbf{v}} & . \end{array} \right. $$

We then have: $\forall \ell \in \mathcal{I}_K$,

$$ (4.7) \qquad (\Pi_{RT_k} \mathbf{v})|_{K_\ell}(\mathbf{x}) = \mathbb{B}_\ell (\hat{\Pi}_{RT_k} \mathbb{B}_\ell^{-1} \hat{\mathbf{v}}_\ell) \circ T_\ell^{-1}(\mathbf{x}) \quad \text{where} \quad \hat{\mathbf{v}}_\ell = \mathbf{v} \circ T_\ell(\hat{\mathbf{x}}). $$

The proof is based on the equality of the $\hat{F}$ and $\hat{K}$-moments of $(\Pi_{RT_k} \mathbf{v})|_{K_\ell} \circ T_\ell(\hat{\mathbf{x}})$ and $\mathbb{B}_\ell (\hat{\Pi}_{RT_k} \mathbb{B}_\ell^{-1} \hat{\mathbf{v}}_\ell)(\hat{\mathbf{x}})$. For $k=0$, setting for $d' \in \{1, \dots, d\}$: $\psi_{f,d'} := \psi_f e_{d'}$, we obtain that:

$$ (4.8) \quad \forall \ell \in \mathcal{I}_K, \forall f \in \mathcal{I}_{F,\ell}, \quad (\Pi_{RT_0} \psi_{f,d'})|_{K_\ell} = (d|K_\ell|)^{-1} (\mathbf{x} - \vec{O}S_{f,\ell}) S_{f,\ell} \cdot \mathbf{e}_{d'}, $$

where $S_{f,\ell}$ is the vertex opposite to $F_f$ in $K_\ell$.

For $k=1$, the vector $(\Pi_{RT_1} \mathbf{v}_h)|_{K_\ell}$ is described by eight unknowns:

$$ (\Pi_{RT_1} \mathbf{v}_h)|_{K_\ell} = \mathbb{A}_\ell \mathbf{x} + (\mathbf{b}_\ell \cdot \mathbf{x}) \mathbf{x} + \mathbf{d}_\ell, \text{ where } \mathbb{A}_\ell \in \mathbb{R}^{2 \times 2}, \mathbf{b}_\ell \in \mathbb{R}^2, \mathbf{d}_\ell \in \mathbb{R}^2. $$

We compute only once the inverse of the matrix of the linear system (4.6), in $\mathbb{R}^{8 \times 8}$.

4.2. **Application with manufactured solutions.** In the Tables 1, 2 and 3, we call $\varepsilon_0'(\mathbf{u}) = \|\mathbf{u} - \mathbf{u}_h\|_{\mathbf{L}^2(\Omega)} / \|(\mathbf{u}, p)\|_\chi$ the velocity error in $\mathbf{L}^2(\Omega)$-norm, where $\mathbf{u}_h$ is the solution to Problem (2.39) (columns $\mathbf{X}_{CR}$ and $\mathbf{X}_{FS}$) or (4.3) (columns $\mathbf{X}_{CR} + \Pi_{RT_0}$ and $\mathbf{X}_{FS} + \Pi_{RT_1}$) and $h$ is the mesh size.

We first consider Stokes problem (2.3) in $\Omega = (0,1)^2$ with $\mathbf{u} = 0, p = (x_1)^3 + (x_2)^3 - 0.5, \mathbf{f} = \mathbf{grad} p = 3 ((x_1)^2, (x_2)^2)^T$. We report in Table 1 $\varepsilon_0'(\mathbf{u})$ for $h = 5.00 e - 2$ and for different values of $\nu$.

<table><caption>TABLE 1. Values of $\varepsilon_0'(\mathbf{u})$ for $h = 5.00 e - 2$</caption><thead><tr><th>ν</th><th>X<sub>CR</sub></th><th>X<sub>CR</sub> + Π<sub>RT<sub>0</sub></sub></th><th>X<sub>FS</sub></th><th>X<sub>FS</sub> + Π<sub>RT<sub>1</sub></sub></th></tr></thead><tbody><tr><td>1.00 e - 4</td><td>7.96 e - 4</td><td>4.59 e - 17</td><td>8.81 e - 7</td><td>1.54 e - 16</td></tr><tr><td>1.00 e - 5</td><td>7.96 e - 4</td><td>4.59 e - 17</td><td>8.81 e - 7</td><td>1.54 e - 16</td></tr><tr><td>1.00 e - 6</td><td>7.96 e - 4</td><td>4.59 e - 17</td><td>8.81 e - 7</td><td>1.54 e - 16</td></tr></tbody></table>

Here we have $\|(\mathbf{u}, p)\|_\chi = \nu \|p\|_{L^2(\Omega)}$. Hence, the $\mathbf{L}^2(\Omega)$-norm of the discrete velocity $\|\mathbf{u}_h\|_{\mathbf{L}^2(\Omega)}$ is proportional to $\nu^{-1}$. Using the projection, we obtain $\varepsilon_0'(\mathbf{u}) = 0$ close to machine precision.

We now consider Stokes problem (2.3) in $\Omega = (0,1)^2$ with:

$$ \mathbf{u} = \begin{pmatrix} (1 - \cos(2 \pi x_1)) \sin(2 \pi x_2) \\ (\cos(2 \pi x_2) - 1) \sin(2 \pi x_1) \end{pmatrix}, \quad p = \sin(2 \pi x_1) \text{ and } \mathbf{f} = -\nu \Delta \mathbf{u} + \mathbf{grad} p. $$
Acknowledgements

The author thanks Mahran Rihani and Albéric Lefort, and she is grateful to the reviewers for their
careful assessment.

References

[1] M. Ainsworth, A. Allendes, G. R. Barrenechea, R. Rankin, "Computable error bounds for nonconforming Fortin-Soulie finite element approximation of the Stokes problem", *IMA J. Numer. Anal.* **32** (2011), no. 2, p. 417-447.

[2] T. Apel, V. Kempf, A. Linke, C. Merdon, "A nonconforming pressure-robust finite element method for the Stokes equations on anisotropic meshes", *IMA J. Numer. Anal.* **42** (2022), no. 1, p. 392-416.

[3] T. Apel, S. Nicaise, J. Schöberl, "Crouzeix-Raviart type finite elements on anisotropic meshes", *Numer. Math.* **89** (2001), no. 2, p. 193-223.

[4] À. Baran, G. Stoyan, "Gauss-Legendre elements: a stable, higher order non-conforming finite element family", *Computing* **79** (2007), no. 1, p. 1-21.

[5] M. Barré, P. Ciarlet Jr., "The T-coercivity approach for mixed problems", *Comptes Rendus. Mathématique* **362** (2024), p. 1051-1088.

[6] C. Bernardi, M. Costabel, M. Dauge, V. Girault, "Continuity properties of the inf-sup constant for the divergence", *SIAM J. Math. Anal.* **48** (2016), no. 2, p. 1250-1271.

[7] C. Bernardi, F. Hecht, "More pressure in the finite element discretization of the Stokes problem", *ESAIM, Math. Model. Numer. Anal.* **34** (2000), no. 5, p. 953-980.

[8] D. Boffi, F. Brezzi, M. Fortin, *Mixed and hybrid finite element methods and applications*, Springer-Verlag, 2013.

[9] A.-S. Bonnet-Ben Dhia, P. Ciarlet Jr., "Notes de cours sur les méthodes variationnelles pour l'analyse et la résolution de problèmes non coercifs", M.Sc. AMS lecture notes (ENSTA-IPP), 2024, https://hal.science/hal-03764649/.

[10] C. Brennecke, A. Linke, C. Merdon, J. Schöberl, "Optimal and pressure independent $L^2$ velocity error estimates for a modified Crouzeix-Raviart Stokes element with BDM reconstructions", *J. Comput. Math.* **33** (2015), no. 2, p. 191-208.

[11] S. C. Brenner, "Poincaré-Friedrichs inequalities for piecewise $H^1$ functions", *SIAM J. Numer. Anal.* **41** (2003), no. 1, p. 306-324.

[12] F. Brezzi, J. Douglas, L. D. Marini, "Two families of mixed finite elements for second order elliptic problems", *Numer. Math.* **47** (1985), no. 2, p. 217-235.

[13] E. Burman, P. Hansbo, "Stabilized Crouzeix-Raviart element for the Darcy-Stokes problem", *Numer. Methods Partial Differ. Equations* **21** (2005), no. 5, p. 986-997.

[14] C. Carstensen, S. Sauter, "Critical functions and inf-sup stability of Crouzeix-Raviart elements", *Comput. Math. Appl.* **108** (2022), p. 12-23.

[15] —, "Crouzeix-Raviart triangular elements are inf-sup stable", *Math. Comput.* **91** (2022), no. 337, p. 2041-2057.

[16] L. Chesnel, P. Ciarlet Jr., "*T*-coercivity and continuous Galerkin methods: application to transmission problems with sign changing coefficients", *Numer. Math.* **124** (2013), no. 1-29.

[17] P. Ciarlet Jr., "T-coercivity: Application to the discretization of Helmhotz-like problems", *Comput. Math. Appl.* **64** (2012), no. 1, p. 22-24.

[18] —, "Analysis of the Scott-Zhang interpolation in the fractional order Sobolev spaces", *J. Numer. Math.* **21** (2013), no. 3, p. 173-180.

[19] P. Ciarlet Jr., E. Jamelot, "Explicit T-coercivity for the Stokes problem: a coercive finite element discretization", to appear in *Comput. Math. Appl.*, 2025, https://inria.hal.science/hal-04414789.

[20] P. Ciarlet Jr., E. Jamelot, F. D. Kpadonou, "Domain decomposition methods for the diffusion equation with low-regularity solution", *Comput. Math. Appl.* **74** (2017), no. 10, p. 2369-2384.

[21] M. Costabel, M. Dauge, "On the inequalities of Babuška-Aziz, Friedrichs and Horgan-Payne", *Arch. Ration. Mech. Anal.* **217** (2015), p. 873-898.

[22] M. Crouzeix, P.-A. Raviart, "Conforming and nonconforming finite element methods for solving the stationary Stokes equations", *Rev. Franç. Autom. Inform. Rech. Opérat.* **7** (1973), no. 3, p. 33-75.

[23] E. A. Dari, R. G. Durán, C. Padra, "Error estimators for nonconforming finite element approximations of the Stokes problem", *Math. Comput.* **64** (1995), no. 211, p. 1017-1033.

[24] L. Diening, J. Storn, T. Tscherpel, "Fortin operator for Taylor-Hood element", *Numer. Math.* **150** (2022), no. 2, p. 671-689.

[25] W. Dörfler, M. Ainsworth, "Reliable a posteriori error control for nonconforming finite element approximation of Stokes flow", *Math. Comput.* **74** (2005), no. 252, p. 1599-1619.

[26] A. Ern, J.-L. Guermond, *Finite elements I*, Texts in Applied Mathematics, vol. 72, Springer, 2021.

[27] —, *Finite elements II*, Texts in Applied Mathematics, vol. 73, Springer, 2021.
We report in Table 2 (resp. 3) the values of $\epsilon_0^\nu(\mathbf{u})$ in the case $\nu = 1.00 e - 3$ (resp. $\nu = 1.00 e - 4$) for mesh sizes. We observe that when there is no projection, $\epsilon_0^\nu(\mathbf{u})$ is independent of $\nu$, whereas using the projection, $\epsilon_0^\nu(\mathbf{u})$ is proportional to $\nu$.

<table><caption>TABLE 2. Values of $\epsilon_0^\nu(\mathbf{u})$ for $\nu = 1.00 e - 3$</caption><thead><tr><th>$h$</th><th>$\mathbf{X}_{CR}$</th><th>$\mathbf{X}_{CR} + \Pi_{RT_0}$</th><th>$\mathbf{X}_{FS}$</th><th>$\mathbf{X}_{FS} + \Pi_{RT_1}$</th></tr></thead><tbody><tr><td>5.00 e - 2</td><td>1.32 e - 3</td><td>2.74 e - 5</td><td>4.73 e - 6</td><td>5.05 e - 7</td></tr><tr><td>2.50 e - 2</td><td>3.30 e - 4</td><td>6.93 e - 6</td><td>5.06 e - 7</td><td>6.42 e - 8</td></tr><tr><td>1.25 e - 2</td><td>8.25 e - 5</td><td>1.74 e - 6</td><td>6.31 e - 8</td><td>8.10 e - 9</td></tr><tr><td>6.25 e - 3</td><td>2.04 e - 5</td><td>4.35 e - 7</td><td>7.44 e - 9</td><td>1.03 e - 9</td></tr><tr><td>Rate</td><td>$h^{2.00}$</td><td>$h^{1.99}$</td><td>$h^{3.08}$</td><td>$h^{2.97}$</td></tr></tbody></table>

<table><caption>TABLE 3. Values of $\epsilon_0^\nu(\mathbf{u})$ for $\nu = 1.00 e - 4$</caption><thead><tr><th>$h$</th><th>$\mathbf{X}_{CR}$</th><th>$\mathbf{X}_{CR} + \Pi_{RT_0}$</th><th>$\mathbf{X}_{FS}$</th><th>$\mathbf{X}_{FS} + \Pi_{RT_1}$</th></tr></thead><tbody><tr><td>5.00 e - 2</td><td>1.32 e - 3</td><td>2.74 e - 6</td><td>4.70 e - 6</td><td>5.05 e - 8</td></tr><tr><td>2.50 e - 2</td><td>3.30 e - 4</td><td>6.93 e - 7</td><td>5.10 e - 7</td><td>6.43 e - 9</td></tr><tr><td>1.25 e - 2</td><td>8.25 e - 5</td><td>1.74 e - 7</td><td>6.37 e - 8</td><td>8.11 e - 10</td></tr><tr><td>6.25 e - 3</td><td>2.04 e - 5</td><td>4.36 e - 8</td><td>7.51 e - 9</td><td>9.77 e - 11</td></tr><tr><td>Rate</td><td>$h^{2.00}$</td><td>$h^{1.99}$</td><td>$h^{3.08}$</td><td>$h^{2.99}$</td></tr></tbody></table>

Let us consider Stokes problem (2.3) with a low-regular velocity. Let $\Omega = (0, 1)^2$, $S_0 = (0.5, 0.5)$, and $(r, \theta)$ be the polar coordinates centred on $S_0$. We set:

$$ \mathbf{u} = r^{\alpha}\mathbf{e}_{\theta}, \quad p = r \quad \text{so that} \quad \mathbf{f} := -\nu \Delta \mathbf{u} + \text{grad} \, p = \nu(1 - \alpha^2) r^{\alpha-2} \mathbf{e}_{\theta} + \mathbf{e}_r. $$

We report in Table 4 the values of $\epsilon_0^\nu(\mathbf{u})$ for $\nu = 1.00 e - 4$, and for different for mesh sizes, with $\alpha = 1$ and $\alpha = 0.49$. For $\alpha = 1$, $\mathbf{u} = (-y, x)^T \in \mathbf{H}^2(\Omega)$. For $\alpha = 0.49$, $\mathbf{u} \in \bigcap_{0<s<\alpha} \mathbf{H}^{1+s}(\Omega)$, hence $\mathbf{u} \notin \mathbf{H}^2(\Omega)$. It seems that the Raviart-Thomas projection is less efficient in that last case.

In order to enhance the numerical results, one could also use a posteriori error esti-

<table><caption>TABLE 4. Values of $\epsilon_0^\nu(\mathbf{u})$, regular and low-regular velocity, $\nu = 1.00 e - 4$.</caption><thead><tr><th></th><th colspan="2">α = 1</th><th colspan="2">α = 0.45</th></tr><tr><th>$h$</th><th>$\mathbf{X}_{FS}$</th><th>$\mathbf{X}_{FS} + \Pi_{RT_1}$</th><th>$\mathbf{X}_{FS}$</th><th>$\mathbf{X}_{FS} + \Pi_{RT_1}$</th></tr></thead><tbody><tr><td>1.00 e - 1</td><td>3.03 e - 5</td><td>2.81 e - 6</td><td>3.05 e - 5</td><td>3.94 e - 6</td></tr><tr><td>5.00 e - 2</td><td>4.34 e - 6</td><td>1.54 e - 6</td><td>4.57 e - 6</td><td>2.15 e - 6</td></tr><tr><td>2.50 e - 2</td><td>4.72 e - 7</td><td>2.41 e - 8</td><td>9.70 e - 7</td><td>8.52 e - 7</td></tr><tr><td>Rate</td><td>$h^{3.00}$</td><td>$h^{3.43}$</td><td>$h^{2.48}$</td><td>$h^{1.11}$</td></tr></tbody></table>

mators to adapt the mesh near point $S_0$ (see [40, 41] for $k = 1$ and [42] for $k = 2$). Alternatively, using the nonconforming Crouzeix-Raviart mixed finite element method, one can build a divergence-free basis, as described in [43] for $k = 1$. When $k = 1$, following the initial work of [44], one can also add $P^1$-Lagrange basis functions to the space of the discrete pressures as explained in [45]. The discrete velocity is
[28] M. S. Fabien, J. Guzmán, M. Neilan, A. Zytoon, "Low-order divergence-free approximations for the Stokes problem on Worsey-Farin and Powell-Sabin splits", *Comput. Methods Appl. Mech. Eng.* **390** (2022), p. 114444.

[29] M. Fortin, M. Soulie, "A non-conforming piecewise quadratic finite element on triangles", *Int. J. Numer. Methods Eng.* **19** (1983), no. 4, p. 505-520.

[30] D. Gallistl, "Rayleigh-Ritz approximation of the inf-sup constant for the divergence", *Math. Comput.* **88** (2019), no. 315, p. 73-89.

[31] G. N. Gatica, *A Simple Introduction to the Mixed Finite Element Method: Theory and Applications*, SpringerBriefs in Mathematics, Springer, 2014.

[32] V. Girault, P.-A. Raviart, *Finite element methods for Navier-Stokes equations*, Springer-Verlag, 1986.

[33] L. Giret, "Non-Conforming Domain Decomposition for the Multigroup Neutron SPN Equation", PhD Thesis, Université Paris-Saclay, 2018, https://hal.science/tel-01936967/.

[34] F. Hecht, "Construction d'une base de fonctions $P_1$ non conforme à divergence nulle dans $\mathbb{R}^3$", *RAIRO, Anal. Numér.* **15** (1981), no. 2, p. 119-150.

[35] E. Jamelot, "Nonconforming mixed finite elements code to solve Stokes Problem, 2D, $k = 1, 2$", https://github.com/cea-trust-platform/stokes_ncfem, 2022.

[36] E. Jamelot, P. Ciarlet Jr., "Fast non-overlapping Schwarz domain decomposition methods for solving the neutron diffusion equation", *J. Comput. Phys.* **241** (2013), p. 445-463.

[37] E. Jamelot, P. Ciarlet Jr., S. Sauter, "Stability of the $P_1$$_{nc}$ – ($P^0$ + $P^1$) element", to appear in the ENUMATH 2023 proceedings, 2025, https://hal.science/hal-04414894.

[38] A. Linke, "On the Role of the Helmholtz-Decomposition in Mixed Methods for Incompressible Flows and a New Variational Crime", *Comput. Methods Appl. Mech. Engrg.* **268** (2014), p. 782-800.

[39] P.-A. Raviart, J.-M. Thomas, "A mixed finite element method for second order elliptic problems", in *Mathematical aspects of finite element methods*, Lecture Notes in Mathematics, vol. 606, Springer, 1977, p. 292-315.

[40] S. Sauter, "The inf-sup constant for hp-Crouzeix-Raviart triangular elements", preprint, 2022, https://arxiv.org/abs/2204.01270v2.

[41] ——, "The inf-sup constant for hp-Crouzeix-Raviart triangular elements", *Comput. Math. Appl.* **149** (2023), p. 49-70.

[42] S. Sauter, C. Torres, "On the Inf-Sup Stability of Crouzeix-Raviart Stokes Elements in 3D", *Math. Comput.* **92** (2023), no. 341, p. 1033-1059.

[43] L. R. Scott, M. S. Vogelius, "Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials", *RAIRO, Modélisation Math. Anal. Numér.* **19** (1985), no. 1, p. 111-143.

[44] L. R. Scott, S. Zhang, "Finite element interpolation of nonsmooth functions satisfying boundary conditions", *Math. Comp.* **54** (1990), p. 483-493.

[45] R. Stenberg, "Error analysis of some finite element methods for the Stokes problem", *Math. Comp.* **54** (1990), p. 495-508.

[46] L. Tartar, *An introduction to Sobolev spaces and interpolation spaces*, Lecture Notes of the Unione Matematica Italiana, vol. 3, Springer, 2007.

[47] C. Taylor, T. Hood, "A numerical solution of the Navier-Stokes equations using the finite element technique", *Comput. Fluids* **1** (1973), p. 73-100.

[48] M. Vohralík, "On the discrete Poincaré-Friedrichs inequities for nonconforming approximations of the Sobolev space $H^1$", *Numer. Funct. Anal. Optim.* **26** (2005), p. 925-952.

[49] S. Zhang, "A new family of stable mixed finite elements for the 3D Stokes equations", *Math. Comput.* **74** (2005), p. 543-554.
then such that $u_h \in H(\text{div}; 0; \Omega)$. Notice that using conforming finite elements, the Scott-Vogelius finite elements [46, 47, 48] produce velocity approximations that are exactly divergence free.

The code used to get the numerical results can be downloaded on GitHub [49].

## 5. CONCLUSION

We analysed the discretization of Stokes problem with nonconforming finite elements in light of the T-coercivity theory. Furthermore, we obtained local stability estimates for order 1 in 2 or 3 dimension without mesh regularity assumption; and for order 2 in 2 dimension in the case of a shape-regular triangulation sequence. This local approach, splitting the normal and the tangential components could help to generalize our results to order $k \ge 3$ (using maybe also other internal moment conservation). This is ongoing work. We then provided numerical results to illustrate the importance of using $H(\text{div})$-conforming projection. Further, we intend to extend the study to other mixed finite element methods.

## ACKNOWLEDGEMENTS

The author acknowledges Mahran Rihani and Albéric Lefort.

## REFERENCES

[1] V. Girault and P.-A. Raviart. *Finite element methods for Navier-Stokes equations*. Springer-Verlag, 1986.

[2] M. Crouzeix and P.-A. Raviart. Conforming and nonconforming finite element methods for solving the stationary Stokes equations. *RAIRO, Sér. Anal. Numér.*, 7(3):33-75, 1973.

[3] M. Fortin and M. Soulie. A non-conforming piecewise quadratic finite element on triangles. *International Journal for Numerical Methods in Engineering*, 19(4):505-520, 1983.

[4] P. Ciarlet Jr. T-coercivity: Application to the discretization of Helmholtz-like problems. *Computers & Mathematics with Applications*, 64(1):22-24, 2012.

[5] E. Jamelot and P. Ciarlet, Jr. Fast non-overlapping Schwarz domain decomposition methods for solving the neutron diffusion equation. *Journal of Computational Physics*, 241:445-463, 2013.

[6] P. Ciarlet Jr., E. Jamelot, and F. D. Kpadonou. Domain decomposition methods for the diffusion equation with low-regularity solution. *Computers & Mathematics with Applications*, 74(10):2369-2384, 2017.

[7] L. Giret. *Non-Conforming Domain Decomposition for the Multigroup Neutron SPN Equation*. PhD thesis, Université Paris-Saclay, 2018.

[8] A. Ern and J.-L. Guermond. *Finite elements II*, volume 73 of *Texts in Applied Mathematics*. Springer, 2021.

[9] L. Chesnel and P. Ciarlet Jr. T-coercivity and continuous Galerkin methods: application to transmission problems with sign changing coefficients. *Numerische Mathematik*, 124:1-29, 2013.

[10] C. Bernardi, M. Costabel, M. Dauge, and V. Girault. Continuity properties of the inf-sup constant for the divergence. *SIAM Journal on Mathematical Analysis*, 48(2):1250-1271, 2016.

[11] D. Gallistl. Rayleigh-Ritz approximation of the inf-sup constant for the divergence. *Mathematics of Computation*, 88(315):73-89, 2019.

[12] M. Costabel and M. Dauge. On the inequalities of Babuška-Aziz, Friedrichs and Horgan-Payne. *Arch. Rational Mech. Anal.*, 217:873-898, 2015.

[13] A.-S. Bonnet-Ben Dhia and P. Ciarlet Jr. Méthodes variationnelles pour l'analyse de problèmes non coercifs, 2023. M.Sc. AMS lecture notes (ENSTA-IPP).

[14] M. Barré and P. Ciarlet Jr. The T-coercivity approach for mixed problems. to appear in *C. R. Acad. Sci. Paris, Ser. I.*, 2024.
[15] P. Ciarlet Jr. and E. Jamelot. The *T*-Coercivity Approach For Solving Stokes Problem: Stabilization of Finite Element Pairs, 2024.

[16] C. Taylor and T. Hood. Numerical solution of the Navier-Stokes equations using the finite element technique. *Computers & Fluids*, 1:73-100, 1973.

[17] D. Boffi, F. Brezzi, and M. Fortin. *Mixed and hybrid finite element methods and applications*. Springer-Verlag, 2013.

[18] A. Ern and J.-L. Guermond. *Finite elements I*, volume 72 of *Texts in Applied Mathematics*. Springer, 2021.

[19] S. Sauter. The inf-sup constant for hp-Crouzeix-Raviart triangular elements, 2022.

[20] S. C. Brenner. Poincaré-Friedrichs inequalities for piecewise $H^1$ functions. *SIAM J. Numer. Anal.*, 41(1):306-324, 2003.

[21] M. Vohralík. On the discrete Poincaré-Friedrichs inequities for nonconforming approximations of the Sobolev space $H^1$. *Numer. Funct. Anal. Optim.*, 26:925-952, 2005.

[22] E. Burman and P. Hansbo. Stabilized Crouzeix-Raviart element for the Darcy-Stokes problem. *Numerical Methods for Partial Differential Equations*, 21(5):986-997, 2005.

[23] S. Sauter. The inf-sup constant for hp-Crouzeix-Raviart triangular elements. *Computers and Mathematics with Applications*, 149:49-70, 2023.

[24] T. Apel, S. Nicaise, and J. Schöberl. Crouzeix-Raviart type finite elements on anisotropic meshes. *Numerische Mathematik*, 89(2):193-223, 2001.

[25] C. Carstensen and S. Sauter. Critical functions and inf-sup stability of Crouzeix-Raviart elements. *Computers and Mathematics with Applications*, 108:12-23, 2022.

[26] L. R. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions. *Math. Comp.*, 54:483-493, 1990.

[27] P. Ciarlet Jr. Analysis of the Scott-Zhang interpolation in the fractional order Sobolev spaces. *Journal of Numerical Mathematics*, 21(3):173-180, 2013.

[28] L. Tartar. *An introduction to Sobolev spaces and interpolation spaces*, volume 3 of *Lecture Notes of the Unione Matematica Italiana*. Springer, 2007.

[29] S. Sauter and C. Torres. On the Inf-Sup Stability of Crouzeix-Raviart Stokes Elements in 3D, 2023.

[30] R. Stenberg. Error analysis of some finite element methods for the Stokes problem. *Math. Comp.*, 54:495-508, 1990.

[31] Á. Baran and G. Stoyan. Gauss-Legendre elements: a stable, higher order non-conforming finite element family. *Computing*, 79(1):1-21, 2007.

[32] C. Carstensen and S. Sauter. Crouzeix-Raviart triangular elements are inf-sup stable, 2022.

[33] L. Diening, J. Storn, and T. Tscherpel. Fortin operator for Taylor-Hood element. *Numerische Mathematik*, 150(2):671-689, 2022.

[34] A. Linke. On the Role of the Helmholtz-Decomposition in Mixed Methods for Incompressible Flows and a New Variational Crime. *Comput. Methods Appl. Mech. Engrg.*, 268:782-800, 2014.

[35] P.-A. Raviart and J.-M. Thomas. A mixed finite element method for second order elliptic problems. In *Mathematical aspects of finite element methods*, volume 606 of *Lecture Notes in Mathematics*, pages 292-315. Springer, 1977.

[36] F. Brezzi, J. Douglas, and L. D. Marini. Two families of mixed finite elements for second order elliptic problems. *Numerische Mathematik*, 47(2):217-235, 1985.

[37] C. Brennecke, A. Linke, C. Merdon, and J. Schöberl. Optimal and pressure independent $L^2$ velocity error estimates for a modified Crouzeix-Raviart Stokes element with BDM reconstructions. *Journal of Computational Mathematics*, 33(2):191-208, 2015.

[38] T. Apel, V. Kempf, A. Linke, and C. Merdon. A nonconforming pressure-robust finite element method for the Stokes equations on anisotropic meshes. *IMA Journal of Numerical Analysis*, 42(1):392-416, 2022.

[39] G. N. Gatica. *A Simple Introduction to the Mixed Finite Element Method: Theory and Applications*. SpringerBriefs in Mathematics. Springer, 2014.

[40] E. Dari, R. Durán, and C. Padra. Error estimators for nonconforming finite element approximations of the Stokes problem. *Mathematics of Computation*, 64(211):1017-1033, 1995.

[41] W. Dörfler and M. Ainsworth. Reliable a posteriori error control for nonconforming finite element approximation of Stokes flow. *Mathematics of Computation*, 74(252):1599-1619, 2005.
[42] M. Ainsworth, A. Allendes, G. R. Barrenechea, and R. Rankin. Computable error bounds for nonconforming Fortin–Soulie finite element approximation of the Stokes problem. *IMA Journal of Numerical Analysis*, 32(2):417–447, 2011.

[43] F. Hecht. Construction d'une base de fonctions $P_1$ non conforme à divergence nulle dans $\mathbb{R}^3$. *RAIRO, Sér. Anal. Numér.*, 15(2):119–150, 1981.

[44] C. Bernardi and F. Hecht. More pressure in the finite element discretization of the Stokes problem. *ESAIM: M2AN*, 34(5):953–980, 2000.

[45] E. Jamelot, P. Ciarlet, and S. Sauter. Stability of the $P_{nc}^1 - (P^0+P^1)$ element. In *ENUMATH 2023 - The European Conference on Numerical Mathematics and Advanced Applications*, Lisbonne, Portugal, 2023.

[46] L. R. Scott and M. Vogelius. Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials. *RAIRO, Sér. Anal. Numér.*, 19(1):111–143, 1985.

[47] S. Zhang. A new family of stable mixed finite elements for the 3D Stokes equations. *Mathe-matics of Computation*, 74:543–554, 2005.

[48] M. Fabien, J. Guzmán, M. Neilan, and A. Zytoon. Low-order divergence-free approxima-tions for the Stokes problem on Worsey–Farin and Powell–Sabin splits. *Computer Methods in Applied Mechanics and Engineering*, 390:114444, 2022.

[49] E. Jamelot. Nonconforming mixed finite elements code to solve Stokes Problem, 2D, $k = 1, 2$, 2022.