A Quasi-Optimal Crouzeix-Raviart Discretization of the Stokes Equations

TL;DR

This paper introduces a modified Crouzeix-Raviart discretization for the Stokes equations that achieves quasi-optimal error bounds, independent of pressure errors, with computational costs comparable to standard methods.

Contribution

It proposes a quasi-optimal discretization of the Stokes equations that maintains error bounds independent of pressure and is computationally efficient.

Findings

Velocity error is proportional to the best approximation error.

Velocity and pressure can be computed separately in 2D domains.

Numerical experiments confirm theoretical error bounds.

Abstract

We present a modification of the Crouzeix-Raviart discretization of the Stokes equations in arbitrary dimension which is quasi-optimal, in the sense that the error of the discrete velocity field in a broken $H^1$-norm is proportional to the error of the best approximation to the analytical velocity field. In particular, the velocity error is independent of the pressure error and the discrete velocity field is element-wise solenoidal. Moreover, the sum of the velocity error times the viscosity plus the pressure $L^2$-error is proportional to the sum of the respective best errors. All proportionality constants are bounded in terms of shape regularity and do not depend on the viscosity. For simply connected two-dimensional domains, the velocity and pressure can be computed separately. The modification only affects the right-hand side aka load vector. The cost for building the modified load vector is proportional to the cost for building the standard load vector. Some numerical experiments illustrate our theoretical results.

Tables

Table 1. Table 6.1. Ratios γ n ​ ( 𝐮 ) subscript 𝛾 𝑛 𝐮 \gamma_{n}(\mathbf{u}) from § 6.1 and γ n m ​ ( 𝐮 ) superscript subscript 𝛾 𝑛 𝑚 𝐮 \gamma_{n}^{m}(\mathbf{u}) , m ∈ { 10 , 20 , 40 } 𝑚 10 20 40 m\in\{10,20,40\} , from § 6.2 for the standard ( std ) and the modified ( mod ) Crouzeix-Raviart discretizations.

		m = 10	m = 20	m = 40
n	std    mod	std    mod	std    mod	std    mod
2	1.37    2.07	1.39    2.03	1.39    2.03	1.39    2.03
3	1.48    2.06	1.50    2.04	1.50    2.04	1.50    2.04
4	1.54    2.05	1.55    2.05	1.55    2.05	1.55    2.05
5	1.57    2.05	1.57    2.05	1.57    2.05	1.57    2.05
6	1.58    2.05	1.58    2.05	1.58    2.05	1.58    2.05

Table 2. Table 6.3. Values δ n ​ ( 𝐮 ) subscript 𝛿 𝑛 𝐮 \delta_{n}(\mathbf{u}) and δ n ​ ( p ) subscript 𝛿 𝑛 𝑝 \delta_{n}(p) from § 6.4 and corresponding EOCs.

n	$\#{\breve{𝒯}}_n$	${\delta }_n\left(𝐮\right)$	${\mathrm{EOC}}_n\left(𝐮\right)$	${\delta }_n\left(p\right)$	${\mathrm{EOC}}_n\left(p\right)$
3	512	6.092e-02		4.339e-02
4	2048	3.673e-02	0.37	2.571e-02	0.38
5	8192	2.135e-02	0.39	1.455e-02	0.41
6	32768	1.206e-02	0.41	8.021e-03	0.43
7	131072	6.670e-03	0.43	4.349e-03	0.44

Equations

$$a(u,v)=\ell (v)$$

$$a_{\mathcal{T}}(u_{\mathcal{T}},v_{\mathcal{T}})={\ell }_{\mathcal{T}}(v_{\mathcal{T}})$$

$$a_{\mathcal{T}}(u_{\mathcal{T}},v_{\mathcal{T}})=\ell (E{v}_{\mathcal{T}})$$

$$a_{\mathcal{T}}(w_{\mathcal{T}},E{v}_{\mathcal{T}})=a_{\mathcal{T}}(w_{\mathcal{T}},v_{\mathcal{T}})$$

$${\Vert u-u_{\mathcal{T}}\Vert }_{\mathcal{T}}\le c_{\mathrm{qopt}}\underset{v_{\mathcal{T}}\in V_{\mathcal{T}}}{\inf }{\Vert u-v_{\mathcal{T}}\Vert }_{\mathcal{T}}.$$

$$\nu {\int }_{\Omega }\nabla \mathbf{u}:\nabla \mathbf{v}-{\int }_{\Omega }p\mathrm{div}\mathbf{v}$$

$${\int }_{\Omega }q\mathrm{div}\mathbf{u}$$

$$\begin{array}{rl}\nu {\int }_{\Omega }{\nabla }_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}:{\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}-{\int }_{\Omega }{p}_{\mathcal{T}}{\mathrm{div}}_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}& ={\int }_{\Omega }\mathbf{f}\cdot {\mathbf{v}}_{\mathcal{T}}\forall {\mathbf{v}}_{\mathcal{T}}\in CR(\mathcal{T})\\ {\int }_{\Omega }{q}_{\mathcal{T}}{\mathrm{div}}_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}& =0\forall q_{\mathcal{T}}\in S^{0,-1}(\mathcal{T})\cap L_0^2(\Omega ).\end{array}$$

$$\begin{array}{c}V=\{\mathbf{u}\in H_0^1(\Omega {)}^d:\mathrm{div}\mathbf{u}=0\}\\ a(\mathbf{u},\mathbf{v})=\nu {\int }_{\Omega }\nabla \mathbf{u}:\nabla \mathbf{v}\text{and}\ell (\mathbf{v})={\int }_{\Omega }\mathbf{f}\cdot \mathbf{v}.\end{array}$$

$$\begin{array}{c}{V}_{\mathcal{T}}=\{{\mathbf{u}}_{\mathcal{T}}\in CR(\mathcal{T}):{\mathrm{div}}_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}=0\}\\ a_{\mathcal{T}}({\mathbf{u}}_{\mathcal{T}},{\mathbf{v}}_{\mathcal{T}})=\nu {\int }_{\Omega }{\nabla }_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}:{\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}\text{and}{\ell }_{\mathcal{T}}({\mathbf{v}}_{\mathcal{T}})={\int }_{\Omega }\mathbf{f}\cdot {\mathbf{v}}_{\mathcal{T}}.\end{array}$$

$$\begin{array}{rl}\nu {\int }_{\Omega }{\nabla }_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}:{\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}-{\int }_{\Omega }{p}_{\mathcal{T}}{\mathrm{div}}_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}& ={\int }_{\Omega }\mathbf{f}\cdot E{\mathbf{v}}_{\mathcal{T}}\forall {\mathbf{v}}_{\mathcal{T}}\in CR(\mathcal{T})\\ {\int }_{\Omega }{q}_{\mathcal{T}}{\mathrm{div}}_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}& =0\forall q_{\mathcal{T}}\in S^{0,-1}(\mathcal{T})\cap L_0^2(\Omega )\end{array}$$

$$\begin{array}{c}{V}_{\mathcal{T}}=\{{\mathbf{u}}_{\mathcal{T}}\in CR(\mathcal{T}):{\mathrm{div}}_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}=0\}\\ a_{\mathcal{T}}({\mathbf{u}}_{\mathcal{T}},{\mathbf{v}}_{\mathcal{T}})=\nu {\int }_{\Omega }{\nabla }_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}:{\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}\text{and}{\ell }_{\mathcal{T}}({\mathbf{v}}_{\mathcal{T}})={\int }_{\Omega }\mathbf{f}\cdot E{\mathbf{v}}_{\mathcal{T}}.\end{array}$$

$${\int }_F{\mathbf{v}}_{\mathcal{T}|K_{F1}}=\left|F\right|{\mathbf{v}}_{\mathcal{T}|K_{F1}}(m_F)=\left|F\right|{\mathbf{v}}_{\mathcal{T}|K_{F2}}(m_F)={\int }_F{\mathbf{v}}_{\mathcal{T}|K_{F2}},$$

$${\int }_{\Omega }{\nabla }_{\mathcal{T}}{\mathbf{w}}_{\mathcal{T}}:{\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}=\sum _{K\in \mathcal{T}}\sum _{F\in {\mathcal{F}}_K}{\int }_F\nabla ({\mathbf{w}}_{\mathcal{T}}{)}_{|K}{\mathbf{n}}_K\cdot {\mathbf{v}}_{\mathcal{T}}$$

$${\int }_{\Omega }{\nabla }_{\mathcal{T}}{\mathbf{w}}_{\mathcal{T}}:{\nabla }_{\mathcal{T}}E{\mathbf{v}}_{\mathcal{T}}=\sum _{K\in \mathcal{T}}\sum _{F\in {\mathcal{F}}_K}{\int }_F\nabla ({\mathbf{w}}_{\mathcal{T}}{)}_{|K}{\mathbf{n}}_K\cdot E{\mathbf{v}}_{\mathcal{T}}$$

$${\int }_{F}E{\mathbf{v}}_{\mathcal{T}}={\int }_F{\mathbf{v}}_{\mathcal{T}}$$

$$A{\mathbf{v}}_{\mathcal{T}}=\sum _{z\in {\mathcal{V}}_{\Omega }}{\lambda }_z{\mathbf{v}}_{\mathcal{T}|K_z}(z).$$

$${\psi }_F=\frac{(2d-1)!}{(d-1)!\left|F\right|}\prod _{z\in {\mathcal{V}}_F}{\lambda }_z.$$

$$B{\mathbf{v}}_{\mathcal{T}}=\sum _{F\in {\mathcal{F}}_{\Omega }}{\psi }_F{\int }_F{\mathbf{v}}_{\mathcal{T}}.$$

$$C{\mathbf{v}}_{\mathcal{T}}=A{\mathbf{v}}_{\mathcal{T}}+B({\mathbf{v}}_{\mathcal{T}}-A{\mathbf{v}}_{\mathcal{T}}).$$

$${\int }_{F}C{\mathbf{v}}_{\mathcal{T}}={\int }_{F}B{\mathbf{v}}_{\mathcal{T}}+\underset{=0}{\underbrace{{\int }_F(A{\mathbf{v}}_{\mathcal{T}}-BA{\mathbf{v}}_{\mathcal{T}})}}={\int }_F{\mathbf{v}}_{\mathcal{T}}$$

$${\int }_K\mathrm{div}(C{\mathbf{v}}_{\mathcal{T}})={\int }_K{\mathrm{div}}_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}$$

$$\mathrm{div}(E{\mathbf{v}}_{\mathcal{T}})={\mathrm{div}}_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}$$

$${\int }_K\nabla {\mathbf{u}}_K:\nabla \mathbf{v}-{\int }_K{p}_K\mathrm{div}\mathbf{v}$$

$${\int }_{K}q\mathrm{div}{\mathbf{u}}_K$$

$${\Vert \nabla S_{K}r\Vert }_K\le c{\Vert r\Vert }_K.$$

$$(S_{K}r{)}_{|F}=0\text{and}\mathrm{div}(S_{K}r)=r$$

$$E{\mathbf{v}}_{\mathcal{T}}=C{\mathbf{v}}_{\mathcal{T}}-\sum _{K\in \mathcal{T}}{S}_K\left(\mathrm{div}(C{\mathbf{v}}_{\mathcal{T}})-{\mathrm{div}}_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}\right)$$

$$\Vert \nabla \mathbf{u}-{\nabla }_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}\Vert \le c_1\underset{{\mathbf{v}}_{\mathcal{T}}\in V_{\mathcal{T}}}{\inf }\Vert \nabla \mathbf{u}-{\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}\Vert =c_1\underset{{\mathbf{v}}_{\mathcal{T}}\in CR(\mathcal{T})}{\inf }\Vert \nabla \mathbf{u}-{\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}\Vert .$$

$${\Vert \mathbf{f}\Vert }_{C{R}^{\ast }}:=\underset{{\mathbf{v}}_{\mathcal{T}}\in CR(\mathcal{T})\setminus \{0\}}{\sup }\frac{{\int }_{\Omega }\mathbf{f}\cdot {\mathbf{v}}_{\mathcal{T}}}{\Vert {\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}\Vert }<\infty .$$

$$\begin{array}{rl}\nu {\int }_{\Omega }{\nabla }_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}:\nabla (E{\mathbf{v}}_{\mathcal{T}})-{\int }_{\Omega }{p}_{\mathcal{T}}\mathrm{div}(E{\mathbf{v}}_{\mathcal{T}})& =\nu {\int }_{\Omega }{\nabla }_{\mathcal{T}}{\mathbf{u}}_{\mathcal{T}}:{\nabla }_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}-{\int }_{\Omega }{p}_{\mathcal{T}}{\mathrm{div}}_{\mathcal{T}}{\mathbf{v}}_{\mathcal{T}}\\ & ={\int }_{\Omega }\mathbf{f}\cdot E{\mathbf{v}}_{\mathcal{T}}\\ & =\nu {\int }_{\Omega }\nabla \mathbf{u}:\nabla (E{\mathbf{v}}_{\mathcal{T}})-{\int }_{\Omega }p\mathrm{div}(E{\mathbf{v}}_{\mathcal{T}})\end{array}$$

Full text

1. Introduction

The Crouzeix-Raviart discretization [8] is a well-established nonconforming finite element method for the Stokes equations. It consists of piecewise constant pressures and piecewise affine velocity fields that are continuous at the barycentres of inter-element faces and that vanish at the barycentres of element faces on the boundary. This discretization is inf-sup stable on general simplicial meshes and provides a first-order approximation of both the analytical velocity and the analytical pressure in a rather straight-forward way. Other remarkable properties are that the discrete velocity field is element-wise solenoidal aka locally conservative and that, on simply connected two-dimensional domains, the velocity and pressure can be computed separately (cf. [3], [5, §VI.8] and §5). Still, one issue is that the nonconformity leads to a consistency error that can be bounded only under regularity assumptions on the analytical solution or the load functional, cf. [20]. Moreover, the velocity error depends on the pressure error, as pointed out by Linke [14].

Combining the approaches proposed in [14] and [21], we construct a modified Crouzeix-Raviart discretization, which is quasi-optimal and pressure-robust. More precisely, the error of the discrete velocity field in a broken $H^{1}$-norm is proportional to the error of the best approximation to the analytical velocity field (cf. Theorem 4.2). Moreover, the error of the discrete pressure in the $L^{2}$-norm is bounded, up to a constant, by the error of the best approximation to the analytical pressure plus the velocity error times the viscosity (cf. Theorem 4.6). All involved constants are independent of the viscosity and bounded in terms of the shape parameter of the underlying mesh. Thus, the error of the discrete velocity field is independent of the pressure, similarly to [14] and unlike the standard discretization of Crouzeix and Raviart. Moreover, our estimates do not involve additional regularity, in contrast to the ones proved by Linke et al. [14, 15].

The standard and the new discretization only differ in the right-hand side aka load vector. Consequently, the discrete velocity field is element-wise solenoidal also in the modified discretization and can be computed separately from the discrete pressure on simply connected two-dimensional domains. The additional cost for computing the modified load vector is proportional to the cost for computing the standard load vector.

The main idea to construct the modified load vector can be described as follows. We employ the smoothing operator of [21, 22] to map Crouzeix-Raviart vector fields into piecewise polynomial, globally continuous vector fields before applying the analytical load functional. Since this operator does not guarantee pressure-robustness, we correct it by a locally computable stable right inverse of the divergence. More precisely, we solve a discrete Stokes problem with Scott-Vogelius elements [11, 16, 18, 19, 23, 24] on a barycentric refinement of each mesh element (cf. (4.7) and (4.10)). Using the controvariant Piola transformation, it is actually sufficient to solve a fixed number of such local problems on a reference configuration (cf. §5). The resulting operator has the additional property that element-wise solenoidal vector fields are mapped into exactly solenoidal vector fields. As already observed in [14], this is decisive to achieve pressure-robustness.

The remainder of this article is organized as follows. In §2 we briefly recall the abstract setting and the main result of [21]. In §§3 and 4 we then present the standard Crouzeix-Raviart discretization and our modification. In §5 we discuss the realization and the additional costs of our modification, as well as the possibility to decouple the computation of velocity and pressure. Finally, in §6 we illustrate and complement our abstract results by means of some numerical experiments in dimension two.

2. Quasi-optimal nonconforming methods

For completeness and a better understanding, we outline the strategy of [21] to design quasi-optimal nonconforming methods for symmetric elliptic problems, in a form adapted to our needs.

Consider a Hilbert space $V$ with scalar product $a(\cdot,\cdot)$. Given any continuous linear functional $\ell$ on $V$, we are looking for the unique function $u\in V$ such that

[TABLE]

holds for all $v\in V$.

For the discretization of problem (2.1), we consider a finite-dimensional space $V_{\mathcal{T}}$. We assume that $a$ extends to a scalar product $a_{\mathcal{T}}(\cdot,\cdot)$ on $V+V_{\mathcal{T}}$ and denote by $\left\lVert\cdot\right\rVert_{\mathcal{T}}=a_{\mathcal{T}}(\cdot,\cdot)^{\frac{1}{2}}$ the induced energy norm. We also replace the load $\ell$ of (2.1) by a linear functional $\ell_{\mathcal{T}}$ on $V_{\mathcal{T}}$. Thus, we look for the unique function $u_{\mathcal{T}}\in V_{\mathcal{T}}$ such that

[TABLE]

holds for all $v_{\mathcal{T}}\in V_{\mathcal{T}}$. We say that this is a (possibly) nonconforming discretization of (2.1), because $V_{\mathcal{T}}$ is not required to be a subspace of $V$.

In standard nonconforming discretizations, like the one in §3 below, the definition of $\ell_{\mathcal{T}}$ often requires that $\ell$ in (2.1) has some extra-regularity. Hence, we cannot extend it to all continuous functionals on $V$. This generates a consistency error that cannot be bounded only in terms of the dual norm of $\ell$, cf. [20, Remark 4.9]. Consequently, the error $\left\lVert u-u_{\mathcal{T}}\right\rVert_{\mathcal{T}}$ is not proportional to the best approximation error $\inf_{v_{\mathcal{T}}\in V_{\mathcal{T}}}\left\lVert u-v_{\mathcal{T}}\right\rVert_{\mathcal{T}}$.

To overcome this drawback, Veeser and Zanotti suggest to consider a linear operator $E:V_{\mathcal{T}}\longrightarrow V$ and look for the function $u_{\mathcal{T}}\in V_{\mathcal{T}}$ such that

[TABLE]

holds for all $v_{\mathcal{T}}\in V_{\mathcal{T}}$. In this context, $E$ is also called smoothing operator, as the nonconformity of $V_{\mathcal{T}}$ often arises from a lack of smoothness. Notice that the discrete load $\ell_{\mathcal{T}}$ in (2.3) is well-defined for all continuous linear functionals on $V$. If $E$ also satisfies

[TABLE]

for all $w_{\mathcal{T}},v_{\mathcal{T}}\in V_{\mathcal{T}}$, then the solution of (2.3) is quasi-optimal (in the norm $\left\lVert\cdot\right\rVert_{\mathcal{T}}$), in the sense that

[TABLE]

Furthermore, the best value of the constant $c_{\mathrm{qopt}}$ is the operator norm of $E$, see [21, Corollary 2.7].

Remark 2.1** (Quasi-optimality and related notions).**

Quasi-optimality extends the well-known Céa’s lemma and implies that the error $\left\lVert u-u_{\mathcal{T}}\right\rVert_{\mathcal{T}}$ of (2.3) is equivalent to the best error $\inf_{v_{\mathcal{T}}\in V_{\mathcal{T}}}\left\lVert u-v_{\mathcal{T}}\right\rVert_{\mathcal{T}}$, in view of the inclusion $u_{\mathcal{T}}\in V_{\mathcal{T}}$, see e.g. [4, Section 2.8] or [20]. We mention that other authors call quasi-optimal estimates in the form $\left\lVert u-u_{\mathcal{T}}\right\rVert_{\mathcal{T}}\leq c_{\mathrm{qopt}}\inf_{v_{\mathcal{T}}\in V_{\mathcal{T}}}\left(\left\lVert u-v_{\mathcal{T}}\right\rVert_{\mathcal{T}}+\mathrm{AG}_{\mathcal{T}}(u-v_{\mathcal{T}})\right)$, where the norm $\left\lVert\cdot\right\rVert_{\mathcal{T}}$ is augmented with $\mathrm{AG}_{\mathcal{T}}(\cdot)$. The augmentation typically involves additional regularity of the solution $u$ or the load $\ell$ and is of higher-order, see Carstensen et al. [6, 7] and Linke et al. [15]. Such estimates are weaker than (2.5), in that the additional regularity required by $\mathrm{AG}_{\mathcal{T}}(\cdot)$ obstructs a further bound of the right-hand side solely in terms of the best error $\inf_{v_{\mathcal{T}}\in V_{\mathcal{T}}}\left\lVert u-v_{\mathcal{T}}\right\rVert_{\mathcal{T}}$.

3. The standard Crouzeix-Raviart discretization

In what follows, $\Omega$ is a connected bounded polyhedron in $\mathbb{R}^{d}$, $d\geq 2$, with Lipschitz-continuous boundary. We denote by $\left\lVert\cdot\right\rVert$ the $L^{2}$-norm on $\Omega$. A subscript to $\left\lVert\cdot\right\rVert$ indicates that we consider the $L^{2}$-norm only on the set specified by the subscript. Let $\mathcal{T}$ be a shape-regular, face-to-face simplicial partition of $\Omega$. We denote by $\mathcal{F}$ and $\mathcal{V}$ the faces and vertices, respectively, of the elements in $\mathcal{T}$. A subscript to $\mathcal{F}$ or $\mathcal{V}$ indicates that we consider only those faces or vertices that are contained in the set specified by the subscript. We denote by $\left\lvert K\right\rvert$ and $\left\lvert F\right\rvert$ the $d$-dimensional Lebesgue measure of an element $K$ and the $(d-1)$-dimensional Hausdorff measure of a face $F$, respectively. We write $h_{K}$ for the diameter of an element $K$. We associate with $\mathcal{T}$ a so-called barycentric refinement $\mathcal{M}$, which is obtained by connecting the vertices and the barycentre of every element $K\in\mathcal{T}$, see [11]. $\mathcal{M}_{K}$ stands for the restriction of $\mathcal{M}$ to an element $K\in\mathcal{T}$ and consists of $d+1$ simplices. We denote by $c$ a generic nondecreasing function of the shape parameter of $\mathcal{T}$. Such function does not need to be the same at different occurrences.

For any integer $k\geq 0$, we denote by $\mathbb{P}_{k}$ the space of polynomials of degree at most $k$ and set $S^{k,-1}(\mathcal{T})=\{\varphi\in L^{2}(\Omega):\forall K\in\mathcal{T}\;\;\varphi_{|K}\in\mathbb{P}_{k}\}$, where $\varphi_{|K}$ denotes the restriction of $\varphi$ to $K$. If $k\geq 1$ we set $S^{k,0}(\mathcal{T})=S^{k,-1}(\mathcal{T})\cap C(\Omega)\subset H^{1}(\Omega)$. The spaces $S^{k,-1}(\mathcal{M}_{K})$ and $S^{k,0}(\mathcal{M}_{K})$ are defined similarly with $\mathcal{T}$ replaced by $\mathcal{M}_{K}$. The (vector-valued) Crouzeix-Raviart space $CR(\mathcal{T})$ consists of all vector fields in $S^{1,-1}(\mathcal{T})^{d}$ that are continuous at the barycentres of interior faces and that vanish at the barycentres of boundary faces. Note that $CR(\mathcal{T})\not\subset H^{1}_{0}(\Omega)^{d}$ due to the missing global continuity and the violation of the boundary condition.

Denoting by $L^{2}_{0}(\Omega)$ the space of all $L^{2}$-functions with mean value zero on $\Omega$ and by $:$ and $\cdot$ the inner products of matrices and vectors respectively, the standard variational formulation of the Stokes equations with viscosity $\nu$ and load $\mathbf{f}\in L^{2}(\Omega)^{d}$ consists in finding $\mathbf{u}\in H^{1}_{0}(\Omega)^{d}$ and $p\in L^{2}_{0}(\Omega)$ such that

[TABLE]

Similarly, denoting by $\nabla_{\mathcal{T}}$ and $\operatorname{div}_{\mathcal{T}}$ the element-wise gradient and divergence, respectively, the standard Crouzeix-Raviart discretization of problem (3.1) consists in finding $\mathbf{u}_{\mathcal{T}}\in CR(\mathcal{T})$ and $p_{\mathcal{T}}\in S^{0,-1}(\mathcal{T})\cap L^{2}_{0}(\Omega)$ such that

[TABLE]

Problems (3.1) and (3.2) each admit a unique solution, see [5, Example II.1.1, Example VI.3.10] or [8]. The discrete velocity field $\mathbf{u}_{\mathcal{T}}$ is element-wise solenoidal, i.e. $\operatorname{div}_{\mathcal{T}}\mathbf{u}_{\mathcal{T}}=0$. If the analytical velocity field $\mathbf{u}$ is in $H^{2}(\Omega)^{d}$ and the analytical pressure $p$ is in $H^{1}(\Omega)$, then the error of the velocity measured in the broken $H^{1}$-norm $\left\lVert\nabla_{\mathcal{T}}\cdot\right\rVert$ and the error of the pressure measured in the $L^{2}$-norm decay linearly with respect to the mesh-size [3, 8]. In proving this estimate, one has to cope with the consistency error of the discretization due to the missing global continuity of the discrete velocity fields. Since, contrary to problem (3.1), problem (3.2) is not defined for general $\mathbf{f}\in H^{-1}(\Omega)^{d}$, it is not fully stable, meaning that its consistency error cannot be bounded in terms of the $H^{-1}$-norm of $\mathbf{f}$, cf. [20] and Remark 4.3 below. Moreover, the fact that element-wise solenoidal discrete velocities are in general not exactly solenoidal entails that (3.2) is not pressure-robust, i.e. the velocity $H^{1}$-error depends on the pressure $L^{2}$-error, see [14]. Both issues are addressed in the next section.

Problems (3.1) and (3.2) do not fit into the framework of §2, since they are in saddle-point form. Yet, testing the first equation of (3.1) with divergence-free functions, we obtain a reduced problem for the analytical velocity field, which fits into (2.1) with space, scalar product and load functional

[TABLE]

Similarly, one can reduce (3.2) to a problem for the discrete velocity field alone, which fits into (2.2) with

[TABLE]

Notice that $V_{\mathcal{T}}\not\subset V$ and that the functional $\ell$ can be extended to $\mathbf{f}\in H^{-1}(\Omega)^{d}$ but $\ell_{\mathcal{T}}$ not. The stiffness matrix of problem (2.2) with (3.4) is symmetric positive definite but its condition number grows like $\operatorname{O}(h^{-4})$ if $\mathcal{T}$ is quasi-uniform. If $\Omega$ is a simply connected two-dimensional domain, there is a basis of $V_{\mathcal{T}}$ consisting of locally supported vector fields and the discrete pressure can be computed from the discrete velocity by sweeping through elements (cf. §5 and Algorithm 1).

4. The modified Crouzeix-Raviart discretization

We now propose a modified version of (3.2) and prove our main results.

4.1. Construction of the discretization

Motivated by the abstract framework of §2 and by the discussion in §3, we modify the standard Crouzeix-Raviart discretization (3.2) as follows: Find $\mathbf{u}_{\mathcal{T}}\in CR(\mathcal{T})$ and $p_{\mathcal{T}}\in S^{0,-1}(\mathcal{T})\cap L^{2}_{0}(\Omega)$ such that

[TABLE]

where $E:CR(\mathcal{T})\longrightarrow H^{1}_{0}(\Omega)^{d}$ is a linear operator (i.e. a smoothing operator, in the terminology of §2).

Since (3.2) and (4.1) only differ in the right-hand side, [5, Example VI.3.10] or [8] again imply that problem (4.1) admits a unique solution. In particular, the discrete velocity field $\mathbf{u}_{\mathcal{T}}$ is element-wise solenoidal, showing that $\mathbf{u}_{\mathcal{T}}\in V_{\mathcal{T}}$. Testing the first equation with functions from $V_{\mathcal{T}}$, we derive a reduced problem for $\mathbf{u}_{\mathcal{T}}$, that fits into (2.2) and (2.3) with

[TABLE]

In contrast to (3.4), here $\ell_{\mathcal{T}}$ can be extended to $\mathbf{f}\in H^{-1}(\Omega)^{d}$.

We aim at constructing the operator $E$ so that the error of the discrete velocity field in the broken $H^{1}$-norm is proportional to the best approximation error to the analytical velocity. For this purpose, we preliminarily observe that (4.2) is a nonconforming discretization of (3.3). Hence, the results recalled in §2 indicate that $E$ should map $V_{\mathcal{T}}$ into $V$ and satisfy (2.4). Such conditions are sufficient to achieve (2.5), with the best constant $c_{\mathrm{qopt}}$ given by the operator norm of $E$. For this reason, we also require that $E$ is $H^{1}$-stable, i.e. $\left\lVert\nabla E\mathbf{v}_{\mathcal{T}}\right\rVert\leq c\left\lVert\nabla_{\mathcal{T}}\mathbf{v}_{\mathcal{T}}\right\rVert$ for all $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$.

Since the midpoint-rule is exact for affine functions, we have, for all $F\in\mathcal{F}_{\Omega}$ and $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$,

[TABLE]

where $K_{F1}$ and $K_{F2}$ are the two elements sharing the face $F$ and $m_{F}$ is its barycentre. The same observation reveals $\int_{F}\mathbf{v}_{\mathcal{T}}=0$ for $F\in\mathcal{F}_{\partial\Omega}$. Therefore, the integral $\int_{F}\mathbf{v}_{\mathcal{T}}$ is defined without ambiguity for all faces and vanishes for boundary faces. Integrating by parts element-wise, one can easily check that, for all $\mathbf{w}_{\mathcal{T}},\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$,

[TABLE]

and

[TABLE]

where $\mathbf{n}_{K}$ is the outward unit normal vector of $K$. Thus, since $\nabla_{\mathcal{T}}\mathbf{w}_{\mathcal{T}}$ is element-wise constant, a sufficient condition for (2.4) is

[TABLE]

for all $F\in\mathcal{F}$ and $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$. This identity always holds on faces $F\in\mathcal{F}_{\partial\Omega}$, because of the homogeneous boundary conditions in $H^{1}_{0}(\Omega)^{d}$ and $CR(\mathcal{T})$.

We first construct a vector version of the smoothing operator in [21, §3.2]. For every interior vertex $z\in\mathcal{V}_{\Omega}$, we denote by $\lambda_{z}\in S^{1,0}(\mathcal{T})$ the conforming first-order nodal basis function associated with the evaluation at $z$. We also choose an element $K_{z}\in\mathcal{T}$ such that $z\in K_{z}$ and keep it fixed in what follows. With this notation, we define a simplified averaging operator $A:CR(\mathcal{T})\longrightarrow S^{1,0}(\mathcal{T})^{d}\cap H^{1}_{0}(\Omega)^{d}$ by

[TABLE]

Note that $A\mathbf{v}_{\mathcal{T}}$ is set to zero at the vertices on $\partial\Omega$.

Next, we associate a face-bubble $\psi_{F}\in S^{d,0}(\mathcal{T})$ with every interior face $F\in\mathcal{F}_{\Omega}$

[TABLE]

The function $\psi_{F}$ is normalized so that $\int_{F^{\prime}}\psi_{F}=\delta_{FF^{\prime}}$ for all $F^{\prime}\in\mathcal{F}$. Then, we define a bubble operator $B:CR(\mathcal{T})\longrightarrow S^{d,0}(\mathcal{T})^{d}\cap H^{1}_{0}(\Omega)^{d}$ by

[TABLE]

The normalization of the bubble functions $\psi_{F}$ guarantees that $B$ satisfies (4.3). Yet, this operator is not $H^{1}$-stable, cf. [21, Remark 3.5]. Thus, we combine $A$ and $B$ and define a smoothing operator $C:CR(\mathcal{T})\longrightarrow S^{d,0}(\mathcal{T})^{d}\cap H^{1}_{0}(\Omega)^{d}$ by

[TABLE]

According to [21, Proposition 3.3], this operator is $H^{1}$-stable. Moreover, rearranging terms and exploiting the definition of $B$, we see that

[TABLE]

for all $F\in\mathcal{F}$ and $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$. This confirms that $C$ satisfies (4.3).

Due to the Gauss theorem, identity (4.4) entails, in particular,

[TABLE]

for all $K\in\mathcal{T}$ and $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$. Unfortunately, this does not imply that $C$ maps $V_{\mathcal{T}}$ into $V$, because $\operatorname{div}(C\mathbf{v}_{\mathcal{T}})\in S^{d-1,-1}(\mathcal{T})$, cf. [22, Remark 3.14]. Recall the barycentric refinement $\mathcal{M}$ of $\mathcal{T}$. The crucial point of our modification of the Crouzeix-Raviart discretization is the construction of another $H^{1}$-stable smoothing operator $E:CR(\mathcal{T})\longrightarrow S^{d,0}(\mathcal{M})^{d}\cap H^{1}_{0}(\Omega)^{d}$, correcting $C$, which preserves the validity of (4.4) and additionally satisfies

[TABLE]

for all $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$. Indeed, this condition entails that $E$ maps $V_{\mathcal{T}}$ into $V$. The desired correction is achieved by solving on every element $K\in\mathcal{T}$ a discrete Stokes problem, based on stable Scott-Vogelius elements on $\mathcal{M}_{K}$, with homogeneous momentum equation and suitable inhomogeneous continuity equation involving $\operatorname{div}(C\mathbf{v}_{\mathcal{T}})$.

To make things precise, we consider, for every element $K\in\mathcal{T}$ and every function $r\in L^{2}_{0}(K)$, the following discrete Stokes problem with unit viscosity, which consists in finding $\mathbf{u}_{K}\in S^{d,0}(\mathcal{M}_{K})^{d}\cap H^{1}_{0}(K)^{d}$ and $p_{K}\in S^{d-1,-1}(\mathcal{M}_{K})\cap L^{2}_{0}(K)$ such that

[TABLE]

Denote by $S_{K}$ the mapping $r\mapsto\mathbf{u}_{K}$. As a consequence of [11, Theorem 3.1], problem (4.7) is uniquely solvable (entailing that $S_{K}$ is well-defined) and we have

[TABLE]

Assuming $r\in S^{d-1,-1}(\mathcal{M}_{K})\cap L^{2}_{0}(K)$ and extending both $r$ and $S_{K}r$ to zero outside $K$, we infer also

[TABLE]

for all $F\in\mathcal{F}$. In particular, the latter property entails that the sum $\sum_{K\in\mathcal{T}}S_{K}$ is a stable global right inverse of the divergence, which is defined on a subspace of $S^{d-1,-1}(\mathcal{M})\cap L^{2}_{0}(\Omega)$ and can be computed locally.

In view of (4.5), the restriction of $\operatorname{div}(C\mathbf{v}_{\mathcal{T}})-\operatorname{div}_{\mathcal{T}}(\mathbf{v}_{\mathcal{T}})$ to each element $K\in\mathcal{T}$ is in $S^{d-1,-1}(\mathcal{M}_{K})\cap L^{2}_{0}(K)$. Therefore, we define the announced smoothing operator $E:CR(\mathcal{T})\longrightarrow S^{d,0}(\mathcal{M})^{d}\cap H^{1}_{0}(\Omega)^{d}$ by

[TABLE]

for every $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$. The previously observed $H^{1}$-stability of $C$ and (4.8) ensure that $E$ is $H^{1}$-stable and its norm is bounded only in terms of the shape parameter of $\mathcal{T}$. Furthermore, condition (4.3) is fulfilled thanks to (4.4) and the first part of (4.9), while (4.6) follows from the second part of (4.9).

Remark 4.1** (Divergence-free pairs).**

The Scott-Vogelius pair employed in (4.7) is divergence-free, in the sense that the divergence maps the discrete velocity space onto the discrete pressure space. This property is needed to compute the local right inverses $S_{K}$ of the divergence, which are then used to correct the operator $C$. An alternative construction, based on a different divergence-free pair, can be found in [15]. The analysis of the Scott-Vogelius pair was initiated in [18, 19], where it was pointed out that stability can be obtained only for certain combinations of the mesh geometry and the polynomial degree $k$. The stability on the barycentric refinement of a simplicial mesh was proved by Qin [16] in 2D for $k\geq 2$ and by Zhang [24] in 3D for $k\geq 3$. The recent paper [11] by Guzmán and Neilan generalizes these results in $\mathbb{R}^{d}$ for $k\geq d$.

4.2. Error estimates

Consider the modified Crouzeix-Raviart discretization (4.1) with the operator $E$ from (4.10). The above-mentioned properties of $E$ imply our first main result.

Theorem 4.2** (Quasi-optimal velocity error).**

Denote by $(\mathbf{u},p)$ and $(\mathbf{u}_{\mathcal{T}},p_{\mathcal{T}})$ the unique solutions of problems (3.1) and (4.1), respectively. There is a constant $\operatorname{c}_{1}$, which only depends on the shape parameter of $\mathcal{T}$ and not on the viscosity $\nu$, such that

[TABLE]

Proof.

The first inequality follows from identities (4.3) and (4.6) and the $H^{1}$-stability of $E$, in view of the abstract discussion in section 2 or [21, Corollary 2.7]. The identity of the best errors can be checked with the help of the so-called Crouzeix-Raviart quasi-interpolation operator $I_{\mathcal{T}}:H^{1}_{0}(\Omega)^{d}\to CR(\mathcal{T})$, which is defined by $\int_{F}I_{\mathcal{T}}\mathbf{u}=\int_{F}\mathbf{u}$ for all faces $F\in\mathcal{F}$. Integrating by parts element-wise, we see that $\int_{K}\nabla_{\mathcal{T}}(I_{\mathcal{T}}\mathbf{u})=\int_{K}\nabla\mathbf{u}$ and $\int_{K}\operatorname{div}_{\mathcal{T}}(I_{\mathcal{T}}\mathbf{u})=\int_{K}\operatorname{div}\mathbf{u}$ for all $K\in\mathcal{T}$. The latter identity entails $I_{\mathcal{T}}\mathbf{u}\in V_{\mathcal{T}}$, because $\mathbf{u}\in V$ and $\operatorname{div}_{\mathcal{T}}(I_{\mathcal{T}}\mathbf{u})$ is element-wise constant. Combining this fact with the first identity implies $\inf_{\mathbf{v}_{\mathcal{T}}\in V_{\mathcal{T}}}\left\lVert\nabla\mathbf{u}-\nabla_{\mathcal{T}}\mathbf{v}_{\mathcal{T}}\right\rVert=\left\lVert\nabla\mathbf{u}-\nabla_{\mathcal{T}}(I_{\mathcal{T}}\mathbf{u})\right\rVert=\inf_{\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})}\left\lVert\nabla\mathbf{u}-\nabla_{\mathcal{T}}\mathbf{v}_{\mathcal{T}}\right\rVert$. ∎

Theorem 4.2 combines the advantages of [14, 15] and [21], because the velocity $H^{1}$-error is independent of the pressure and proportional to the error of the best approximation to the analytical velocity. This implies that the velocity error of the standard discretization (3.2) can be smaller, but not arbitrarily smaller, than the error of the modified one. The constant $c_{1}$ quantifies the maximum gap, which is uniformly bounded for shape regular meshes. Conversely, the next remark shows that (4.1) can significantly outperform (3.2) for certain loads.

Remark 4.3** (Nonconforming discretizations without smoothing).**

The error estimate of Theorem 4.2, combined with the triangle inequality, reveals that the solution $\mathbf{u}_{\mathcal{T}}$ of (4.1) depends continuously on the analytical velocity $\mathbf{u}$ in the broken $H^{1}$-norm. This property hinges on the use of the smoothing operator $E$, which maps the Crouzeix-Raviart test functions into $H^{1}_{0}(\Omega)^{d}$, before the application of the load functional. In particular, $\mathbf{u}_{\mathcal{T}}$ is well-defined for all loads $\mathbf{f}\in H^{-1}(\Omega)^{d}$. In contrast, the standard Crouzeix-Raviart discretization (3.2) is only defined under the regularity assumption

[TABLE]

Since $\left\lVert\cdot\right\rVert_{CR^{*}}$ is not equivalent to the $H^{-1}$-norm, the discrete velocity $\mathbf{u}_{\mathcal{T}}$ of (3.2) does not depend continuously on $\mathbf{u}$ in the broken $H^{1}$-norm. Therefore, the velocity $H^{1}$-error of the standard discretization can be arbitrarily larger than the best error $\inf_{\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})}\left\lVert\nabla_{\mathcal{T}}(\mathbf{u}-\mathbf{v}_{\mathcal{T}})\right\rVert$, provided $\left\lVert\mathbf{f}\right\rVert_{CR^{*}}$ is sufficiently larger than the $H^{-1}$-norm of $\mathbf{f}$. This observation is not specific to (3.2) and concerns any other nonconforming discretization without smoothing, including the ones in [14, 15].

The following remarks shed additional light on the error estimate of Theorem 4.2, in connection with the results available for other existing discretizations of the Stokes equations.

Remark 4.4** (Quasi-optimality and divergence-free pairs).**

The quasi-optimal estimate of Theorem 4.2 distinguishes (4.1) from nonconforming methods without smoothing but also from several other conforming finite element methods for (3.1). Indeed, to our best knowledge, estimates in this form have been previously obtained only with conforming and divergence-free pairs, like the one of Scott and Vogelius, see Remark 4.1. Moreover, if we restrict the attention to first-order discretizations, only few such pairs are known to be stable, either requiring special mesh geometries or the use of non-polynomial basis functions, cf. [10, 11, 25]. Theorem 4.2 shows that quasi-optimality is actually not restricted to conforming and divergence-free pairs but can be achieved by means of a proper discretization of the load functional.

Remark 4.5** (Shape regularity and anisotropic meshes).**

Acosta and Durán considered in [1] possibly anisotropic meshes, fulfilling maximum angle conditions, in dimension two and three. For certain smooth solutions of (3.1), they proved that the velocity error of (3.2) in the broken $H^{1}$-norm converges with maximum decay rate. A counterpart of this result cannot be derived for (4.1) by Theorem 4.2. In fact, the constant $\operatorname{c}_{1}$ in (4.11) depends on the shape parameter of $\mathcal{T}$ and it does not seem possible to avoid such dependence if the operator $E$ is defined as in (4.10), cf. §6.2. Still, even in case $\operatorname{c}_{1}$ is indeed large, it is not obvious that (3.2) always outperforms (4.1). This is illustrated by a numerical experiment in §6.

The abstract framework of §2 cannot be used to bound the pressure $L^{2}$-error. For this purpose, we combine the previous bound of the velocity error, the properties of $E$ and standard techniques for saddle point problems. We first observe that integration by parts element-wise, together with (4.3), the first equations of (3.1) and (4.1) and the fact that $E$ maps into $H^{1}_{0}(\Omega)^{d}$ imply

[TABLE]

for all $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$. Thus, we derive the identity

[TABLE]

Next, we compare $p_{\mathcal{T}}$ with the $L^{2}$-projection $\pi_{\mathcal{T}}p\in S^{0,-1}(\mathcal{T})$ of $p$

[TABLE]

and use (4.3) and the Gauss theorem to obtain

[TABLE]

for all $\mathbf{v}_{\mathcal{T}}\in CR(\mathcal{T})$. Finally, the uniform stability of the Crouzeix-Raviart element [5, Example VI.3.10], [8] implies

[TABLE]

where $\beta=\inf_{q\in L^{2}_{0}(\Omega)\setminus\{0\}}\sup_{\mathbf{v}\in H^{1}_{0}(\Omega)^{d}\setminus\{0\}}\frac{\int_{\Omega}q\operatorname{div}\mathbf{v}}{\left\lVert q\right\rVert\left\lVert\nabla\mathbf{v}\right\rVert}$ is the inf-sup constant of the divergence operator. Recall that $\beta$ does not depend on the viscosity $\nu$ and $\beta^{-1}$ is bounded in terms of the ratio $\mathrm{diam}(\Omega)/R$, provided $\Omega$ is star-shaped with respect to a ball of radius $R$, cf. [2]. Consequently, the estimate in Theorem 4.6 below depends on this ratio, while the one in Theorem 4.2 is independent of it.

Combining (4.13), (4.15), (4.14), (4.12), the $H^{1}$-stability of $E$ and Theorem 4.2 proves our second main result.

Theorem 4.6** (Pressure error).**

Denote by $(\mathbf{u},p)$ and $(\mathbf{u}_{\mathcal{T}},p_{\mathcal{T}})$ the unique solutions of problems (3.1) and (4.1), respectively, and let $c_{1}$ be as in Theorem 4.2. There are two constants $\operatorname{c}_{2}$ and $\operatorname{c}_{3}$, which only depend on the shape parameter of $\mathcal{T}$ and the inf-sup constant $\beta$ of the divergence operator but not on the viscosity $\nu$, such that

[TABLE]

Estimates of the pressure error in this form are well-established for all conforming stable pairs and not only for the divergence-free ones, see [5, Theorem II.2.1]. Correspondingly, our proof exploits (4.3) and the $H^{1}$-stability of $E$, but not (4.6). In contrast to similar results in [8] and [14, 15], Theorem 4.6 does not assume additional regularity of the solution.

5. Practical aspects

5.1. Assembly of the modified discretization

Assume that $\{\mathbf{v}_{\mathcal{T}}^{1},\dots,\mathbf{v}_{\mathcal{T}}^{N}\}$ and $\{\mathbf{w}_{\mathcal{M}}^{1},\dots,\mathbf{w}_{\mathcal{M}}^{M}\}$ are bases of $CR(\mathcal{T})$ and $S^{d,0}(\mathcal{M})^{d}\cap H^{1}_{0}(\Omega)^{d}$, respectively. Problem (4.1) requires the computation of $\int_{\Omega}\mathbf{f}\cdot E\mathbf{v}_{\mathcal{T}}^{i}$ for $i=1,\dots,N$. Let $\mathbb{E}\in\mathbb{R}^{N\times M}$ and $\overrightarrow{\mathbf{f}}\in\mathbb{R}^{M}$ be defined as follows

[TABLE]

for $i=1,\dots,N$ and $j=1,\dots,M$. Few algebraic manipulations reveal

[TABLE]

We aim at showing that both $\mathbb{E}$ and $\overrightarrow{\mathbf{f}}$ can be computed with $\operatorname{O}(N)$ operations, if one uses standard nodal bases of $CR(\mathcal{T})$ and $S^{d,0}(\mathcal{M})^{d}\cap H^{1}_{0}(\Omega)^{d}$.

The nodal basis $\{\varphi_{F}\mathbf{e}_{1},\ldots,\varphi_{F}\mathbf{e}_{d}:F\in\mathcal{F}_{\Omega}\}$ of $CR(\mathcal{T})$ is known to be a convenient choice for assembling the left-hand side of (3.2) and (4.1). Here $\{\mathbf{e}_{1},\ldots,\mathbf{e}_{d}\}$ stands for the standard basis of $\mathbb{R}^{d}$ and each $\varphi_{F}\in S^{1,-1}(\mathcal{T})$ is uniquely defined by the following properties: $\varphi_{F}$ is continuous at barycentres of interior faces and is normalized so that $\varphi_{F}(m_{F^{\prime}})=\left\lvert F\right\rvert^{-1}\delta_{FF^{\prime}}$ for all $F^{\prime}\in\mathcal{F}$. Note that the midpoint quadrature rule implies $\int_{F^{\prime}}\varphi_{F}=\delta_{FF^{\prime}}$ for all $F^{\prime}\in\mathcal{F}$. Similarly, we consider the basis $\{\varphi_{z}\mathbf{e}_{1},\ldots,\varphi_{z}\mathbf{e}_{d}:z\in\mathcal{V}_{\Omega}^{d}(\mathcal{M})\}$ of $S^{d,0}(\mathcal{M})^{d}\cap H^{1}_{0}(\Omega)^{d}$, where $\mathcal{V}_{\Omega}^{d}(\mathcal{M})$ denotes the interior Lagrange nodes of degree $d$ of $\mathcal{M}$ and $\varphi_{z}$ is the nodal basis function of $S^{d,0}(\mathcal{M})$ associated with the evaluation at $z$.

With the latter basis, the computation of $\overrightarrow{\mathbf{f}}$ requires $\operatorname{O}(M)=\operatorname{O}(N)$ operations. Moreover, each entry of $\mathbb{E}$ can be computed in $\operatorname{O}(1)$ operations and the maximum number of nonzero entries in each row of $\mathbb{E}$ is bounded by a constant depending on the shape parameter of $\mathcal{T}$. This entails that the cost for computing $\mathbb{E}$ is $\operatorname{O}(N)$ operations. Since the left-hand side of (3.2) is the same as in (4.1), we infer that the cost for assembling the modified discretization is the same as the one for assembling the standard one.

To check our claim concerning the matrix $\mathbb{E}$, we first observe that each basis function $\varphi_{F}\mathbf{e}_{k}$, with $F\in\mathcal{F}_{\Omega}$ and $k\in\{1,\ldots,d\}$, is supported in the union $\omega_{F}$ of the two elements sharing $F$. According to [21, §3.2], the function $C(\varphi_{F}\mathbf{e}_{k})$ is supported in the union of all elements touching $\omega_{F}$ and can be computed with $\operatorname{O}(1)$ operations. Hence, to obtain $E(\varphi_{F}\mathbf{e}_{k})$ from $C(\varphi_{F}\mathbf{e}_{k})$, one has to solve the local Stokes problem (4.7) with load $r=\operatorname{div}(C(\varphi_{F}\mathbf{e}_{k}))-\operatorname{div}_{\mathcal{T}}(\varphi_{F}\mathbf{e}_{k})$ on each element $K\in\mathcal{T}$ such that $K\cap\omega_{F}\neq\emptyset$. For this purpose, an efficient strategy is to precompute the solution of the local problem on a reference element $\widehat{K}$ with load $\widehat{r}$, where $\widehat{r}$ varies in a basis $\{\widehat{r}_{1},\ldots,\widehat{r}_{m}\}$ of $\mathbb{P}_{d-1}\cap L^{2}_{0}(\widehat{K})$. Then, the solution of (4.7) can be simply obtained by means of the controvariant Piola’s transformation, see [5, §III.1.3]. This confirms that $E(\varphi_{F}\mathbf{e}_{k})$ is supported in $\{K\in\mathcal{T}:K\cap\omega_{F}\neq\emptyset\}$ and can be computed from $\varphi_{F}\mathbf{e}_{k}$ with $\operatorname{O}(1)$ operations. We conclude recalling the definition of $\mathbb{E}$ and our choice of the nodal basis of $S^{d,0}(\mathcal{M})^{d}\cap H^{1}_{0}(\Omega)^{d}$.

5.2. Solution of the modified discretization

As for the standard Crouzeix-Raviart discretization (3.2), the computation of the velocity and pressure solving the modified problem (4.1) can be decoupled, if $\Omega$ is a simply connected two-dimensional domain. This procedure is essentially known in the literature. We show how to adapt it to the modified problem, for the sake of completeness. An extension to three dimensions is given in [12].

First recall that the discrete velocity $\mathbf{u}_{\mathcal{T}}$ solving (4.1) is in the subspace $V_{\mathcal{T}}$, consisting of element-wise solenoidal Crouzeix-Raviart functions, cf. (4.2). For a simply connected two-dimensional domain, a basis of $V_{\mathcal{T}}$ is given by the union of the sets $\{\varphi_{F}\mathbf{t}_{F}:F\in\mathcal{F}_{\Omega}\}$ and $\{\mathbf{w}_{z}:z\in\mathcal{V}_{\Omega}\}$, see [5, Example VI.8.1, Figure VI.35]. Here, $\mathbf{t}_{F}$ is a unit tangent vector to $F$ and $\mathbf{w}_{z}$ is a vortex around $z$, defined by

[TABLE]

where $\mathcal{F}_{z}$ consists of all edges meeting at $z$ and $\mathbf{n}_{F_{z}}$ is a unit normal vector to $F_{z}$, oriented counterclockwise with respect to $z$. This basis can be used to compute $\mathbf{u}_{\mathcal{T}}$ from problem (2.3) with (4.2).

Next, let $\mathbf{n}_{F}$, $F\in\mathcal{F}_{\Omega}$, be a unit normal vector to $F$ and denote by ${\mathbb{J}}_{F}(\cdot)$ the jump across $F$ in direction $\mathbf{n}_{F}$. Assuming that $\mathbf{u}_{\mathcal{T}}$ is known, we test the first equation of (4.1) with $\mathbf{v}_{\mathcal{T}}=\varphi_{F}\mathbf{n}_{F}$. Rearranging terms, the Gauss theorem yields

[TABLE]

for all $F\in\mathcal{F}_{\Omega}$. This set of conditions determines the discrete pressure $p_{\mathcal{T}}$ uniquely, because the first equation of (4.1) is automatically fulfilled for all test functions $\mathbf{v}_{\mathcal{T}}=\varphi_{F}\mathbf{t}_{F}$ with $F\in\mathcal{F}_{\Omega}$.

Let $K_{0},K\in\mathcal{T}$ be connected by a path of triangles $(K_{i})_{i=1}^{n}$ such that $K_{n}=K$ and $F_{i}:=K_{i}\cap K_{i-1}$, $i=1,\dots,n$, is an interior edge of $\mathcal{T}$. Since $p_{\mathcal{T}}$ is element-wise constant, we have

[TABLE]

where $\mathbf{n}_{K_{i}}$ is the outward unit normal vector of $K_{i}$. This identity has two interesting consequences. First, the sum $\sum_{i=1}^{n}\mathbf{n}_{F_{i}}\cdot\mathbf{n}_{K_{i}}{\mathbb{J}}_{F_{i}}(p_{\mathcal{T}})$ vanishes for $K=K_{0}$, meaning that the choice of the path connecting two triangles is irrelevant. Second, comparing with (5.1), we see that the value of $p_{\mathcal{T}}$ in any triangle $K\neq K_{0}$ depends only on the load $\mathbf{f}$, the discrete velocity $\mathbf{u}_{\mathcal{T}}$ and the value of $p_{\mathcal{T}}$ in $K_{0}$. Therefore, we can compute $p_{\mathcal{T}}$ as follows. Defining $\overline{p}_{\mathcal{T}}:=p_{\mathcal{T}}-p_{\mathcal{T}|K_{0}}$, we have $\overline{p}_{\mathcal{T}|K_{0}}=0$ and, for all $K\neq K_{0}$,

[TABLE]

where $(K_{i})_{i=1}^{n}$ is any path connecting $K_{0}$ and $K$. The constraint $\int_{\Omega}p_{\mathcal{T}}=0$ further implies $\sum_{K\in\mathcal{T}}\left\lvert K\right\rvert\overline{p}_{\mathcal{T}|K}=-\left\lvert\Omega\right\rvert p_{\mathcal{T}|K_{0}}$, showing that $p_{\mathcal{T}}$ can be easily recovered from $\overline{p}_{\mathcal{T}}$. The following algorithm implements this procedure in a way that requires a number of operations proportional to $\dim(S^{0,-1}(\mathcal{T}))$.

6. Numerical experiments

In this section we report and discuss the results of four numerical experiments, that are intended to illustrate and partially complement Theorems 4.2 and 4.6. In particular, we compare the modified Crouzeix-Raviart discretization (4.1), with $E$ as in (4.10), and the standard one (3.2), whenever the latter is defined. If an exact solution is available, we compute also the best approximation $H^{1}$-error to the analytical velocity and the best approximation $L^{2}$-error to the analytical pressure,

[TABLE]

As mentioned in the proof of Theorem 4.2, the former is given by

[TABLE]

where $I_{\mathcal{T}}$ is the Crouzeix-Raviart quasi-interpolation operator.

All experiments have been implemented in ALBERTA 3.0 [13, 17] and concern the two-dimensional Stokes equations (3.1), posed in the unit square, with unit viscosity, i.e.

[TABLE]

6.1. Smooth solution

We first consider a test case with smooth solution, namely

[TABLE]

where $\operatorname{curl}(v):=(\partial_{x_{2}}v,-\partial_{x_{1}}v)$. We solve (3.2) and (4.1) on the following sequence $(\mathcal{T}_{n})_{n\geq 0}$ of uniform meshes. We divide $\Omega$ into $2^{n}\times 2^{n}$ squares, with edges parallel to the lines $x_{2}=0$ and $x_{1}=0$ and edge length $2^{-n}$. Then, we obtain $\mathcal{T}_{n}$ by drawing, for each square, the diagonal parallel to the line $x_{1}=x_{2}$, see Figure 6.1. Since $(\mathbf{u},p)\in H^{2}(\Omega)^{2}\times H^{1}(\Omega)$, both $e_{\mathcal{T}_{n}}(\mathbf{u})$ and $e_{\mathcal{T}_{n}}(p)$ converge to zero with maximum decay rate $2^{-n}\approx(\#\mathcal{T}_{n})^{-0.5}$.

To assess the quality of the standard and the modified Crouzeix-Raviart discretizations, we compute the ratios

[TABLE]

where $(\mathbf{u}_{\mathcal{T}_{n}},p_{\mathcal{T}_{n}})\in CR(\mathcal{T}_{n})\times S^{0,-1}(\mathcal{T}_{n})$ denotes either the solution of (3.2) or the one of (4.1). Some values of $\gamma_{n}(\mathbf{u})$ and $\gamma_{n}(p)$ are displayed in the first column of Tables 6.2 and 6.2, respectively. They indicate that the velocity $H^{1}$-error of the modified discretization is larger than the one of the standard discretization by a factor between $1.3$ and $1.5$. Instead, the corresponding pressure $L^{2}$-errors are nearly of the same size for sufficiently large $n$ and the error of the modified discretization is smaller for the first values of $n$. More generally, one can expect that both discretizations perform similarly for smooth solutions of (3.1), on shape regular sequences of meshes.

6.2. Smooth solution and anisotropic meshes

The purpose of this experiment is to compare the standard Crouzeix-Raviart discretization (3.2) and the modified one (4.1) on sequences of meshes with increasing shape parameter. To this end, we consider the same exact solution as in §6.1 and the following sequence $(\mathcal{T}_{n}^{m})_{n\geq 0}$ of meshes with prescribed anisotropy $m\in\mathbb{N}$. For any fixed $m$, we divide $\Omega$ into $(m2^{n})\times 2^{n}$ rectangles, with edges parallel to the lines $x_{2}=0$ and $x_{1}=0$. The length of the edges is $2^{-n}/m$ and $2^{-n}$, respectively. Then, we obtain $\mathcal{T}_{n}^{m}$ by drawing, for each rectangle, the diagonal parallel to the line $x_{2}=mx_{1}$, see Figure 6.1. For $m=1$, this is nothing else than the mesh $\mathcal{T}_{n}$ considered in the previous experiment. The diameter of all triangles in $\mathcal{T}_{n}^{m}$ is proportional to $2^{-n}$.

Since the maximum angle in all meshes is $\pi/2$, we expect and observe that the velocity $H^{1}$-error and the pressure $L^{2}$-error of the standard discretization (and, consequently, also $e_{\mathcal{T}_{n}^{m}}(\mathbf{u})$ and $e_{\mathcal{T}_{n}^{m}}(p)$) converge to zero with the maximum decay rate $2^{-n}\approx(\#\mathcal{T}_{n}^{m}/m)^{-0.5}$, irrespective of $m$, cf. [1]. The same result cannot be inferred from the quasi-optimal estimates in Theorems 4.2 and 4.6, because the shape parameter of $\mathcal{T}_{n}^{m}$ is proportional to $m$. Moreover, it does not appear possible to improve such estimates, because we numerically computed the best constant in (4.8) for this type of triangles and observed a linear dependence on $m$. (Recall that this constant enters into the bound of $\operatorname{c}_{1}$ in Theorem 4.2).

Proceeding as before, we compute the ratios

[TABLE]

both for the standard discretization and the modified one. Some values of $\gamma_{n}^{m}(\mathbf{u})$ and $\gamma_{n}^{m}(p)$ are displayed in Tables 6.2 and 6.2, respectively, for $m\in\{10,20,40\}$. They indicate that, in this specific example, the performance of both discretizations for large $m$ remains close to the best possible and is similar to the one observed in §6.1 for $m=1$.

6.3. Rough pressure

This experiment aims at illustrating the pressure-robustness of the modified discretization (4.1). Hence, we consider a test case with smooth analytical veclocity and rough analytical pressure, namely

[TABLE]

We construct the initial triangulation $\widehat{\mathcal{T}}_{0}$ by drawing the diagonals of $\Omega$, see Figure 6.2. The intersection of the diagonals is taken as newest vertex for all triangles in $\widehat{\mathcal{T}}_{0}$. Each one of the successive meshes $\widehat{\mathcal{T}}_{n}$, $n\geq 1$, is obtained from the previous one through two global refinements by newest vertex bisection.

Notice that the load $\mathbf{f}=-\Delta\mathbf{u}+\nabla p$ has a singular part concentrated on the line $l:=\{\pi^{-1}\}\times(0,1)$. Thus, $\mathbf{f}$ does not belong to $L^{2}(\Omega)^{2}$. Still, it is possible to extend the standard Crouzeix-Raviart discretization (3.2) to this case because, for all $n\geq 0$, each edge of $\widehat{\mathcal{T}}_{n}$ intersects $l$ in at most one point.

Since $\mathbf{u}\in H^{2}(\Omega)^{2}$, Theorem 4.2 and standard interpolation estimates entail that $e_{\widehat{T}_{n}}(\mathbf{u})$ and the velocity $H^{1}$-error of (4.1) converge to zero with the maximum decay rate $(\#\widehat{\mathcal{T}}_{n})^{-0.5}$. Moreover, our numerical results indicate that the ratio $\gamma_{n}(\mathbf{u})$, defined as in (6.1) with $\widehat{\mathcal{T}}_{n}$ in place of $\mathcal{T}_{n}$, is nearly 2. In contrast, the data displayed in Figure 6.3 show that the velocity $H^{1}$-error of (3.2) is impaired by the low regularity of the analytical pressure and converges approximately with rate $(\#\widehat{\mathcal{T}}_{n})^{-0.25}$. The pressure $L^{2}$-errors of both discretizations are quite close to the best $L^{2}$-error $e_{\widehat{\mathcal{T}}_{n}}(p)$ and converge approximately with decay rate $(\#\widehat{\mathcal{T}}_{n})^{-0.25}$.

6.4. Rough load

Similarly as in the previous experiment, we now consider the Stokes equations (3.1) with a load $\mathbf{f}\notin L^{2}(\Omega)^{2}$. Indeed, we assume that the action of $\mathbf{f}$ on $H^{1}_{0}(\Omega)^{2}$-functions is given by

[TABLE]

where $\mathbf{t}=(0,1)$. The main difference from §6.3 is that here $\mathbf{f}$ is concentrated on the line $l=\{0.5\}\times(0,1)$ and its density is not constant (nor piecewise constant) in the $1$-dimensional Hausdorff measure.

Let $\breve{\mathcal{T}}_{n}$, $n\geq 0$, be obtained from the mesh $\widehat{\mathcal{T}}_{n}$ in §6.3 through a global refinement by newest vertex bisection, see Figure 6.2. Since some edges of $\breve{\mathcal{T}}_{n}$ are contained in the support $l$ of $\mathbf{f}$, the standard Crouzeix-Raviart discretization $\eqref{E:CR}$ cannot be extended to this case. (The same observation applies also to the discretization proposed in [14, 15].) In contrast, the modified discretization of §4 is well-defined, because $\mathbf{f}\in H^{-1}(\Omega)^{2}$. Thus, we only solve (4.1).

As the analytical solution is not available, we only monitor the difference between the approximations $(\mathbf{u}_{\breve{\mathcal{T}}_{n-1}},p_{\breve{\mathcal{T}}_{n-1}})$ and $(\mathbf{u}_{\breve{\mathcal{T}}_{n}},p_{\breve{\mathcal{T}}_{n}})$, obtained on $\breve{\mathcal{T}}_{n-1}$ and $\breve{\mathcal{T}}_{n}$, respectively. Hence, we compute

[TABLE]

for $n\geq 1$. We estimate the decay rate of $\delta_{n}(\cdot)$ in terms of $\#\breve{\mathcal{T}}_{n}$ through the so-called experimental order of convergence

[TABLE]

for $n\geq 2$. Some values of $\delta_{n}(\mathbf{u})$ and $\delta_{n}(p)$, together with the corresponding EOCs, are displayed in Table 6.3.

The load $\mathbf{f}$ is in $H^{-1/2-\varepsilon}(\Omega)$ for all $\varepsilon>0$ and this implies that the corresponding analytical solution is in $H^{3/2-\varepsilon}(\Omega)\times H^{1/2-\varepsilon}(\Omega)$, owing to the shift theorem of [9]. Thus, one may expect that $\mathrm{EOC}_{n}(\mathbf{u})$ and $\mathrm{EOC}_{n}(p)$ are nearly $0.25$. The values in Table 6.3 actually indicate a higher decay rate, both of $\delta_{n}(\mathbf{u})$ and $\delta_{n}(p)$. A possible explanation is that, focusing for instance on the best $H^{1}$-error to the analytical velocity $\mathbf{u}$, one has, for all $n\geq 0$,

[TABLE]

where $\Omega_{-}=(0,0.5)\times(0,1)$ and $\Omega_{+}=(0.5,1)\times(0,1)$, cf. [21, eq. (1.2)]. Therefore, according to Theorem 4.2, the modified discretization (4.1) is potentially able to exploit additional regularity of $\mathbf{u}_{|\Omega_{-}}$ and $\mathbf{u}_{|\Omega_{+}}$ beyond the one of $\mathbf{u}$.

Acknowledgements

We thank Andreas Veeser for many inspiring discussions and the unknown reviewers for their comments which led to substantial improvements of the present article.