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

HAL
open science

# The pressure-robust virtual element methods with reconstruction for the Stokes equation

Jikun Zhao, Wenhao Zhu, Zhaonan Dong, Shipeng Mao

► To cite this version:

Jikun Zhao, Wenhao Zhu, Zhaonan Dong, Shipeng Mao. The pressure-robust virtual element methods with reconstruction for the Stokes equation. 2025. hal-05288123

HAL Id: hal-05288123

https://inria.hal.science/hal-05288123v1

Preprint submitted on 29 Sep 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
# The pressure-robust virtual element methods with reconstruction for the Stokes equation

Jikun Zhao<sup>a,*</sup>, Wenhao Zhu<sup>a</sup>, Zhaonan Dong<sup>b</sup>, Shipeng Mao<sup>d</sup>

<sup>a</sup>School of Mathematics and Statistics, Zhengzhou University, Zhengzhou 450001, PR China

<sup>b</sup>Inria, 48 Rue Barrault, 75647 Paris, France and CERMICS, Ecole nationale des ponts et chaussées, IP Paris, 6 & 8 avenue B. Pascal, 77455 Marne-la-Vallée, France.

<sup>c</sup>State Key Laboratory of Mathematical Sciences, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, PR China

<sup>d</sup>School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, PR China

## Abstract

Based on the idea of reconstruction, we develop a unified theoretical framework on the pressure-robust virtual element methods for Stokes equation. The reconstruction can be done on any shape-regular simplicial subdivision of polytopal elements and consists of two steps: the interpolation onto $H(\text{div})$-conforming finite element space and the solution of a local problem element by element. By using this reconstruction we only modify the right-hand side of original virtual element methods in order to obtain the pressure-robust schemes with optimal convergence. This framework is applied to the conforming and nonconforming virtual element methods of arbitrary polynomial degree in two and three dimensions. In particular, we present the pressure-robust scheme for the nonconforming virtual element method and give some explanations on developing the pressure-robust conforming virtual element method by using this reconstruction. Numerical experiments confirm the optimal convergence rate for this nonconforming pressure-robust scheme and show that the velocity error is independent of the dynamic viscosity parameter as well as of large pressure gradients.

**Keywords:** Stokes equation; Virtual element method; Pressure-robust; Reconstruction

## 1. Introduction

The incompressible Navier-Stokes equations are fundamental in modeling fluid flow across a wide range of scientific and engineering applications. In regimes where fluid motion is slow, the Stokes equation provides a simplified approximation by neglecting inertial effects. On a bounded polytopal (polygonal/polyhedral) domain $\Omega \subset \mathbb{R}^d$ with $d = 2, 3$, we consider the basic Stokes equation with Dirichlet boundary conditions:

$$ egin{cases} -\nu \Delta \mathbf{u} + \nabla p = \mathbf{f}, & \text{in } \Omega, \\ \text{div } \mathbf{u} = 0, & \text{in } \Omega, \\ \mathbf{u} = \mathbf{0}, & \text{on } \partial\Omega, \end{cases} \quad (1.1) $$

where the positive parameter $\nu$ is the dynamic viscosity, $\mathbf{f}$ represents a given external force, $\mathbf{u}$ is the unknown velocity field and $p$ the pressure. The Stokes equation forms a linear system but requires special attention due to the incompressibility condition $\text{div } \mathbf{u} = 0$, which physically corresponds to the conservation of mass. Moreover, adding a gradient field to the external force alters the pressure solution without affecting the velocity, in symbols, that is

$$ \mathbf{f} \rightarrow \mathbf{f} + \nabla q \quad \Rightarrow \quad (\mathbf{u}, p) \rightarrow (\mathbf{u}, p+q). \quad (1.2) $$

It is highly desirable for numerical approximations to preserve both key properties. Therefore, on the one hand, the numerical approximation $\mathbf{u}_h$ of the velocity should satisfy the divergence-free condition, either

<footnote>*Corresponding author.</footnote>

<footnote>Email addresses: jkzhao@zzu.edu.cn (Jikun Zhao), zhuwenhao@gs.zzu.edu.cn (Wenhao Zhu), zhaonan.dong@inria.fr (Zhaonan Dong), maosp@lsec.cc.ac.cn (Shipeng Mao)</footnote>
in a conforming form $\mathrm{div} \, u_h = 0$ (exactly divergence-free) or at least in a nonconforming form $\mathrm{div}_h u_h = 0$ (discretely divergence-free). On the other hand, the method should be pressure-robust, meaning that the velocity error is independent of the pressure $p$ and the parameter $\nu$, thereby ensuring numerical stability even in the presence of large pressure gradients or small parameter $\nu$. For more detailed information, see [16].

Over the past few decades, various stable finite element methods (FEMs) have been successfully developed and widely applied to the Stokes equation. Among them, finite elements (FEs) such as the Taylor-Hood, mini, and Crouzeix-Raviart (CR) elements are theoretically well-established, but generally do not yield exactly divergence-free velocity fields, which limits their effectiveness in simulating incompressible flows [16, 27]. In contrast, $P_k$-$P_{k-1}$ element pairs are particularly appealing due to their ability to preserve key structural features of the Stokes equation, notably exact divergence-free velocity and pressure-robust behavior. It is worth noting that the low-order $P_1$-$P_0$ element pair does not satisfy the essential inf-sup condition and therefore cannot provide a stable approximation for the Stokes equation [15]. In contrast, the lowest-order CR-$P_0$ element pair naturally satisfies the inf-sup condition and enforces a divergence-free velocity field in a nonconforming sense, making it particularly appealing for engineering applications [12]. However, this lowest-order nonconforming FEM is not pressure-robust [16, 20], meaning that the velocity approximation may deteriorate when the pressure gradient $\nabla p$ is large in scale compared to the velocity $u$ or the positive parameter $\nu \ll 1$. To achieve pressure robustness in these classical mixed FEMs, discretization schemes based on the reconstruction of the right-hand side have recently attracted significant attention. The central idea behind this approach is to construct a computable reconstruction operator that maps discretely divergence-free test functions onto an exactly divergence-free function space while preserving some essential properties. In [19], a reconstruction operator was introduced by mapping the lowest-order CR element onto the Raviart-Thomas (RT) element on general simplex meshes in two and three dimensions. Unlike the lowest-order CR mixed FEM with discontinuous pressure approximations, Taylor-Hood and mini mixed FEMs use continuous pressures space requiring different reconstruction strategies. In order to develop the pressure-robust schemes of Taylor-Hood and mini mixed FEMs, in [17] Lederer et al. designed novel divergence-free reconstruction operators for both elements by solving specific local problems defined on various vertex patches. An important advantage of these strategies mentioned above is that only the right-hand side discretization is modified, leaving the stiffness matrix and the discrete inf-sup condition of the original method unchanged.

In recent years, numerical methods that relax the restrictions on mesh partitioning in FEMs and are applicable to general polytopal meshes have attracted considerable attention. As an extension of traditional FEMs, the virtual element method (VEM) is particularly convenient for handling meshes with hanging nodes by treating them as polytopes, making it highly promising for the development of adaptive algorithms. Meanwhile, the relaxation of meshing requirements opens up new possibilities for the numerical solution of partial differential equations. Unlike conventional FEMs that rely on polynomial shape functions, the VEM allows for non-polynomial shape functions, providing greater flexibility in constructing approximation spaces and enabling the design of divergence-free virtual elements (VEs). However, since VEM shape functions typically lack explicit expressions despite being well-defined in theory, numerical implementations rely on discrete bilinear forms that are carefully designed to be computable, ensuring the consistency and stability of the numerical scheme. A variety of VEs have been developed for the Stokes equation, including $H^1$-conforming Stokes elements [2, 6, 7, 22], $H(\mathrm{div})$-conforming but $H^1$-nonconforming elements [28], as well as CR-type $H^1$-nonconforming elements [26, 27, 29], etc. These VEs have also been successfully extended to related problems such as the Navier-Stokes, Brinkman, and Stokes-Darcy equations.

Although the aforementioned VEMs yield exactly or discretely divergence-free velocities, standard VE discretizations typically lack pressure robustness that is essential for incompressible flow simulations. This is because the body force term ($f, w_h$) is usually replaced by its discrete approximation ($f, \Pi_k^{0,h} w_h$) in the VE discretization, where $\Pi_k^{0,h}$ denotes the standard piecewise $L^2$-projection of degree $k$ applied to the velocity test function $w_h$. While this approach guarantees asymptotic pressure robustness, it nevertheless breaks the inherent pressure robustness (1.2) presented in the continuous Stokes equation (1.1). There have been some works on pressure-robust VEMs in the literature, which can be divided into two classes. The first class is carried out based on the reformulation of the variational problem. For example, Beirão da Veiga et al. [5] have developed a pressure-robust VEM for the Stokes equation with a homogeneous Newtonian fluid and unit viscosity, supporting arbitrary polynomial orders and polyhedral meshes. Their approach is based on reformulating the viscous term as the curl of the vorticity in the momentum balance equation and projecting the body force onto a discrete $H(\mathrm{curl})$ space. Mora et al. [21] proposed a high-order pressure-robust $C^1$-VEM
for solving the Brinkman problem formulated using the stream function. The second class is based on a reconstruction strategy inspired by the pressure-robust FE framework developed in [16, 17, 19, 20]. This strategy reconstructs the right-hand side term ($f, w_h$) as ($f, R_h w_h$), where $R_h$ is a reconstruction operator computable from the degrees of freedom (DoFs) of VE function $w_h$ and designed to preserve its essential properties (e.g. exactly divergence-free). Unlike the pressure-robust VEMs based on unconventional weak formulations proposed in [5, 21], the strategy of introducing operator $R_h$ offers the key advantage that only the discretization of the right-hand side is modified, while the stiffness matrix and the discrete inf-sup condition of the original method remain unchanged. Frerichs et al. [14] designed the reconstruction operator $R_h$ based on a $H(\text{div})$-conforming approximation obtained by solving local minimization problems defined on a RT FE space on local subtriangulations of polygons, where a pressure-robust conforming VEM of arbitrary order was developed. We mention that in [14] the reconstruction proposed by Frerichs et al. is presented in two dimensions and it is stated that the extension to three dimensions is possible. Wang et al. [23] directly constructed a lowest-order operator $R_h$ using a specialized $H(\text{div})$-conforming space defined on convex polygons and presented a lowest-order pressure-robust conforming VEM. Recently, by a similar reconstruction introduced in [23], Liang et al. [18] developed a nonconforming VEM for the incompressible poroelasticity model, which ensures parameter robustness and effectively overcomes Poisson locking on general polygonal meshes. For more related works, see [13, 24, 25].

The design of reconstruction operators can be naturally applied to other related problems to achieve pressure-robust VEMs. However, except for [14], almost all of the existing studies only focus on two-dimensional cases and the construction of low-order reconstruction operators. The construction of computable reconstruction operators typically relies on $H(\text{div})$-conforming FEs defined on simplex-based subtriangulations of polytopes in meshes. Meanwhile, the DoFs in VE spaces often do not directly correspond to those in standard $H(\text{div})$-conforming FE spaces, which significantly increases the complexity of constructing reconstruction operators in practice. This challenge is especially pronounced in three dimensions, where polyhedral boundaries consist of polygons, while simplex boundaries are triangles. Furthermore, the mismatch between the internal DoFs in these spaces further challenges the construction of high-order reconstruction operators. Besides, we mention that in [18, Section 3.2.3] only lowest-order reconstruction operator is designed to obtain the pressure-robust nonconforming VEM, so it is valid only in the lowest-order case.

In this paper, we introduce a unified framework on pressure-robust VEMs for the Stokes equation (1.1) based on the idea of reconstruction. This theoretical framework overcomes the difficulties mentioned above and can be applied to the conforming and nonconforming VEMs of arbitrary degree in two and three dimensions. Like the pressure-robust FEMs [16], we replace the right-hand side of original VEMs by ($f, R_h w_h$) with an elaborate reconstruction operator $R_h$, where the reconstruction can be done on any shaped-regular simplicial triangulation of polytope in meshes. In the abstract framework presented in Section 3, we give some concrete requirements on $R_h$, which contains computability, commutativity with piecewise $L^2$-projection, error estimates in the $L^2$-norm and polynomial orthogonality. The computability is necessary for VEMs. The commutativity ensures that the reconstruction maps a discretely divergence-free function to a exactly divergence-free function. The last two requirements ensure the optimal convergence of VE approximation. The process of reconstruction is divided into two steps presented in Section 4. The first step is to interpolate the VE function onto a $H(\text{div})$-conforming macro FE space element by element, which yields the commutativity and error estimates in the $L^2$-norm but fails to acquire the polynomial orthogonality. The second step is to solve a local problem on each element, of which the solution is used to amend the interpolation in order to achieve the polynomial orthogonality. We remark that the way of reconstruction presented in this paper is variously different from that in [14], although two types of reconstructions have the very similar properties. By contrast, our reconstruction is applicable to a wider range of cases and seems to be easier to implement in two and three dimensions. Further, on general simplicial subtriangulations of polygons such as barycentric subtriangulation, the lowest-order reconstruction presented in [14] requires solving a constrained minimization problem on each element, whereas the lowest-order reconstruction presented in this paper can be achieved only by the interpolation on the $H(\text{div})$-conforming macro FE space.

Compared with the existing works (e.g. [14, 18, 23, 24]), here we make some important progress on the pressure-robust VEMs, which are summed up as follows.

(1) The only requirement on the simplicial subdivision of polytope in meshes is to be shape-regular, where the reconstruction is established with some good properties mentioned above. This mesh condition is general and allows for the most popular polytopal meshes, especially for the three-dimensional meshes,
where the common barycenteric subdivision is feasible.

(2) The reconstruction is exploited in a general framework, such that the resulting pressure-robust VEMs can be of arbitrary degree, containing conforming and nonconforming VEMs in two and three dimensions.

(3) For the first time, the pressure-robust nonconforming VEM of arbitrary degree is developed in two and three dimensions.

Finally, we emphasize that the starting point of our reconstruction is a specialized $H(\text{div})$-conforming macro FE with the divergence on polytope and normal components on $(d-1)$-faces being polynomials, which is constructed via the FE de Rham complex on a shape-regular simplicial triangulation of polytope, see [11] for the details on this macro FE. For VE functions, we define a computable interpolation operator $I_h^{\text{div}}$ by this $H(\text{div})$-conforming macro FE, for which the commutativity is naturally satisfied. We also derive the interpolation error estimates for the operator $I_h^{\text{div}}$ in the $L^2$-norm. However, the interpolation $I_h^{\text{div}}$ does not satisfy the required polynomial $L^2$-orthogonality for higher-order cases ($k \ge 3$). Inspired by the reconstruction idea from [17], we resolve this issue by introducing a specially designed local problem defined on this $H(\text{div})$-conforming macro FE. The well-posedness of this local problem is analyzed by proving an inf-sup condition. Based on the abstract framework on reconstruction, we develop a concrete pressure-robust scheme for the CR-type nonconforming VE in Section 5. In the numerical tests presented in Section 6, we accomplish this reconstruction on a simple simplicial submesh obtained by decomposing each polygon into triangles by connecting its barycenter to the vertices. The numerical results show that the velocity error in the $H^1$-seminorm converges at the optimal order $O(h^k)$ independently of the pressure, while the pressure error in the $L^2$-norm also converges at the optimal order $O(h^k)$. Therein, for a high-order method ($k = 3$) we also verify the necessity of solving the local problem in the reconstruction by a numerical test where only the $H(\text{div})$-conforming FE interpolation $I_h^{\text{div}}$ is used but it fails to achieve the optimal convergence, specially for the velocity approximation.

## 2. Notation and the weak formulation

### 2.1. Notation

For any given bounded domain $D \subset \mathbb{R}^d$ with $d = 2, 3$ and a nonnegative integer $s$, the usual notations are used to stand for both the scalar Sobolev space $H^s(D)$ and the vector Sobolev space $H^s(D; \mathbb{R}^d) := H^s(D) \otimes \mathbb{R}^d$. These spaces are endowed with the norm $\| \cdot \|_{s,D}$ and the seminorm $| \cdot |_{s,D}$, while the notation $(\cdot, \cdot)_D$ is utilized to denote the usual $L^2(D)$ (or $L^2(D; \mathbb{R}^d)$) inner product. Let $H_0^s(D)$ be the closure of $C_0^\infty(D)$ with respect to the norm $\| \cdot \|_{s,D}$, and $L_0^2(D)$ consists of all functions in $L^2(D)$ with zero mean value. If $D$ is $\Omega$, we will get rid of the subscript in norms and inner product.

Let $P_s(D)$ stand for the set of all polynomials over $D$ with the total degree no more than $s$, and let $P_s(D; \mathbb{R}^d) := P_s(D) \otimes \mathbb{R}^d$ denote the vector version. The set of homogeneous polynomials of degree $s$ on $D$ is denoted by $P_s(D)$. We also use the symbol $C$ (or $C_1, C_2, \alpha, \beta, \dots$) to stand for a generic positive constant, of which the values may be changed at different places, but independent of the mesh size and the parameter $\nu$. For convenience, let the short symbol $a \lesssim b$ stand for the inequality $a \le Cb$ with the constant $C$ independent of $a, b$ and the mesh size. Let $a \asymp b$ mean $a \lesssim b$ and $b \lesssim a$.

For a column vector function $w = (w_i)_{1 \le i \le d}$, we define the divergence and gradient by

$$ \operatorname{div} w = \sum_{i=1}^{d} \frac{\partial w_i}{\partial x_i}, \qquad \nabla w = \left( \frac{\partial w_i}{\partial x_j} \right)_{1 \le i,j \le d}. $$

We define the space

$$ \mathbf{H}(\text{div}; D) = \{ w \in L^2(D; \mathbb{R}^d); \operatorname{div} w \in L^2(D) \}, $$

and

$$ \mathbf{H}_0(\text{div}; D) = \{ w \in \mathbf{H}(\text{div}; D); \text{tr} w = 0 \text{ on } \partial D \}, $$

where tr denotes the trace operator for $H(\text{div}; D)$. Moreover, for a scalar function $v$ and a vector function $w$, depending on the space dimension $d$, we will set

$$ \operatorname{curl} v = \left( \frac{\partial v}{\partial x_2}, -\frac{\partial v}{\partial x_1} \right)^{\top}, \quad d = 2, $$
and

$$ \mathrm{curl} \mathbf{w} = \nabla \times \mathbf{w} = \left( \frac{\partial w_3}{\partial x_2} - \frac{\partial w_2}{\partial x_3}, \frac{\partial w_1}{\partial x_3} - \frac{\partial w_3}{\partial x_1}, \frac{\partial w_2}{\partial x_1} - \frac{\partial w_1}{\partial x_2} \right)^{\top}, \quad d=3. $$

We present a decomposition of polynomial vector space. Let $x_* \in \mathbb{R}^d$ be an arbitrary reference point. In two dimensions we define

$$ \mathbf{G}_k^{x*}(D) = (\mathbf{x} - \mathbf{x}_*)^{\perp} \mathbb{P}_{k-1}(D), \quad \mathbf{x}^{\perp} = (-x_2, x_1)^{\top}, $$

and in three dimensions

$$ \mathbf{G}_k^{x*}(D) = (\mathbf{x} - \mathbf{x}_*) \times \mathbb{P}_{k-1}(D; \mathbb{R}^3), \quad \mathbf{x} = (x_1, x_2, x_3)^{\top}. $$

The following space decomposition

$$ \mathbb{P}_k(D; \mathbb{R}^d) = \nabla \mathbb{P}_{k+1}(D) \oplus \mathbf{G}_k^{x*}(D) $$

will be useful in what follows.

## 2.2. The continuous problem

Let us consider the bilinear forms

$$ \begin{aligned} a(\mathbf{v}, \mathbf{w}) &= \nu(\nabla \mathbf{v}, \nabla \mathbf{w}), && \forall \mathbf{v}, \mathbf{w} \in H^1(\Omega; \mathbb{R}^d), \\ b(\mathbf{w}, q) &= -(\mathrm{div} \mathbf{w}, q), && \forall \mathbf{w} \in H^1(\Omega; \mathbb{R}^d), \forall q \in L_0^2(\Omega). \end{aligned} $$

Then the corresponding weak formulation of problem (1.1) reads: find $(\mathbf{u}, p) \in H_0^1(\Omega; \mathbb{R}^d) \times L_0^2(\Omega)$ such that

$$ \left\{ \begin{array}{ll} a(\mathbf{u}, \mathbf{w}) + b(\mathbf{w}, p) = (\mathbf{f}, \mathbf{w}), & \forall \mathbf{w} \in H_0^1(\Omega; \mathbb{R}^d), \\ b(\mathbf{u}, q) = 0, & \forall q \in L_0^2(\Omega). \end{array} \right. \qquad (2.1) $$

It is well-known that the existence and uniqueness of the solution to problem (2.1) is ensured by the subsequent two stability conditions

$$ \begin{gathered} a(\mathbf{w}, \mathbf{w}) \geq \nu\alpha \|\mathbf{w}\|_1^2, \quad \forall \mathbf{w} \in \mathcal{Z}, \\ \sup_{0 \neq \mathbf{w} \in H_0^1(\Omega; \mathbb{R}^d)} \frac{b(\mathbf{w}, q)}{\|\mathbf{w}\|_1} \geq \beta \|q\|_0, \quad \forall q \in L_0^2(\Omega), \end{gathered} $$

where the subspace

$$ \mathcal{Z} = \{ \mathbf{w} \in H_0^1(\Omega; \mathbb{R}^d); b(\mathbf{w}, q) = 0, \forall q \in L_0^2(\Omega) \} = \{ \mathbf{w} \in H_0^1(\Omega; \mathbb{R}^d); \mathrm{div} \, \mathbf{w} = 0 \}. $$

# 3. The strategy on developing pressure-robust VEM

## 3.1. Polytopal meshes and simplicial subdivisions of polytopes

Let $\mathcal{T}_h$ be a decomposition of $\Omega$ into non-overlapping polytopes and $h$ stand for the maximum diameter of elements in $\mathcal{T}_h$. The symbol $\mathcal{V}_h$ denotes the set of all vertices in $\mathcal{T}_h$, $\mathcal{E}_h$ the set of all $(d-2)$-faces in $\mathcal{T}_h$ and $\mathcal{F}_h$ the set of all $(d-1)$-faces in $\mathcal{T}_h$. For an element $K$ in $\mathcal{T}_h$, the diameter of $K$ is denoted by $h_K$. $\mathcal{V}(K)$ denotes the set of vertices of $K$, $\mathcal{E}(K)$ the set of $(d-2)$-faces of $K$ and $\mathcal{F}(K)$ the set of $(d-1)$-faces of $K$. Thus in two dimensions, we have $\mathcal{V}_h = \mathcal{E}_h$ and $\mathcal{V}(K) = \mathcal{E}(K)$. Moreover, let $\mathbf{n}_K$ denote the exterior unit normal vector of the polytope $K \in \mathcal{T}_h$.

For every $(d-1)$-face $F \in \mathcal{F}_h$, the diameter of $F$ is denoted by $h_F$. Let $\mathbf{n}_F$ denote a unit normal vector of $F$, which is given arbitrarily but fixed for interior $(d-1)$-faces and the exterior normal of $\Omega$ for boundary $(d-1)$-faces. For an internal $(d-1)$-face $F$, we define the jump of function $w$ over $F$ by

$$ [w]|_F = (w|_K)|_F - (w|_L)|_F, $$
where $F$ is shared by two polytopes $K$ and $L$ and $n_F$ points from $K$ to $L$. When $F$ is a boundary $(d-1)$-face, we simply set $[w]_F = w|_F$. The following assumption on the mesh $\mathcal{T}_h$ is standard:

**H0.** There exists a positive constant $\gamma$, such that

• for every element $K \in \mathcal{T}_h$ and every face $F \in \mathcal{F}(K)$, the ratio between the diameters $h_F$ and $h_K$ is bigger than $\gamma$,

• each element $K \in \mathcal{T}_h$ and each face $F \in \mathcal{F}_h$ is star-shaped with the uniformly bounded constant $\gamma$.

To develop a pressure-robust VEM with a computable reconstruction operator, we assume that there exists a shape-regular simplicial mesh $\mathcal{T}_h^*$ such that each $K \in \mathcal{T}_h$ is a union of some simplexes in $\mathcal{T}_h^*$. For every polytope $K \in \mathcal{T}_h$, let $\mathcal{T}_K$ be the quasi-uniform simplicial partition of $K$ induced from $\mathcal{T}_h^*$. The symbol $\mathcal{E}(\mathcal{T}_K)$ denotes the set of all $(d-2)$-faces in the simplicial partition $\mathcal{T}_K$ and $\mathcal{F}(\mathcal{T}_K)$ the set of all $(d-1)$-faces in $\mathcal{T}_K$. Hereafter we use $T$ to represent a simplex, and $K$ to denote a general polytope.

Set

$$ \mathcal{F}^{\partial}(\mathcal{T}_K) = \{F \in \mathcal{F}(\mathcal{T}_K); F \subset \partial K\}, \quad \mathcal{F}^{0}(\mathcal{T}_K) = \{F \in \mathcal{F}(\mathcal{T}_K); F \not\subset \partial K\}. $$

There is generally a significant difference between $\mathcal{F}^{\partial}(\mathcal{T}_K)$ and $\mathcal{F}(K)$. The simpler geometry in two dimensions allows one to preserve the original edges of polygonal element as edges in the subtriangulation $\mathcal{T}_K$, that is, we easily find a special triangulation such that

$$ \mathcal{F}^{\partial}(\mathcal{T}_K) = \mathcal{F}(K), $$

whereas such a construction is generally impossible in three dimensions. A significant advantage of this special subdivision is that it effectively avoids disturbing the DoFs associated with the original boundary structure, which greatly facilitates the construction of the reconstruction operator. One possible approach is to employ constrained Delaunay triangulation on the polygon $K$. However, this method is sensitive to the geometric complexity of the polygon, especially in the presence of hanging nodes. Thanks to the assumption of mesh regularity in the VEM, another possible approach is to construct a simplicial subdivision of the polytope $K \in \mathcal{T}_h$ by introducing its barycenter $x_K$ as an interior vertex (see Figure 1 for the examples).

In three dimensions, it is noted that the elements in $\mathcal{F}^{\partial}(\mathcal{T}_K)$ are triangles, whereas those in $\mathcal{F}(K)$ are general polygons. This means that

$$ \mathcal{F}^{\partial}(\mathcal{T}_K) \neq \mathcal{F}(K), $$

which brings significant challenges in designing the reconstruction operator. Consequently, to construct a suitable simplicial subdivision, one has to modify the original geometric structure of the element boundary $\partial K$. At this point, the constrained Delaunay triangulation cannot be directly used to generate simplicial subdivisions of polyhedron. In contrast, the barycentric subdivision additionally involves the barycenters of the faces $F \in \mathcal{F}(K)$ as new vertices in $\mathcal{T}_K$ (see Figure 2), thereby enabling a valid simplicial decomposition based on the assumption of mesh regularity in the VEM. Obviously, the barycentric subdivision is feasible in practical computations.

<figure><img src="image_3.png" /><figcaption>Figure 1: Examples of subtriangulations on polygon $K$.</figcaption></figure>

## 3.2. The common virtual element method

To continue our work, we first present an abstract framework on the common VEMs. Let $W_h$ denote the VE space for the velocity, with $W_h|_K \subset H^1(K; \mathbb{R}^d)$ as the restriction of $W_h$ onto a polytope $K \in \mathcal{T}_h$. The discrete
<figure><img src="image_4.png" /><figcaption>Figure 2: Barycentric subdivision on polyhedron K with barycenter $x_K$ and one of its faces F with barycenter $x_F$.</figcaption></figure>

setting employs the broken gradient $\nabla_h : W_h \to L^2(\Omega; \mathbb{R}^{d \times d})$ and the broken divergence $\text{div}_h : W_h \to L^2(\Omega)$ in the sense that

$$ (\nabla_h w_h)|_K := \nabla(w_h)|_K, \quad (\text{div}_h w_h)|_K := \text{div}(w_h)|_K, \quad \forall w_h \in W_h, \forall K \in \mathcal{T}_h. $$

We assume that the functions in $W_h$ are at least weakly continuous, i.e., for any function $w_h \in W_h$, it holds

$$ \int_F [w_h] \cdot q_h \, ds = 0, \quad q_h \in \mathbb{P}_{k-1}(F; \mathbb{R}^d), \quad \forall F \in \mathcal{F}_h. \quad (3.1) $$

Let us assume that the discrete seminorm $|\cdot|_{1,h}$ defines a norm on $W_h$, where

$$ |w_h|_{1,h}^2 := (\nabla_h w_h, \nabla_h w_h) = \frac{1}{\nu} \sum_{K \in \mathcal{T}_h} a^K(w_h, w_h), \quad \forall w_h \in W_h, $$

and the local bilinear form is defined by

$$ a^K(w_h, w_h) = \nu(\nabla w_h, \nabla w_h)_K. $$

The pressure is approximated using the standard discontinuous piecewise polynomial space

$$ Q_h = \{q_h \in L_0^2(\Omega); q_h|_K \in \mathbb{P}_{k-1}(K), \forall K \in \mathcal{T}_h\}. $$

For any $q \in L_0^2(\Omega)$, the piecewise $L^2$-projection $\Pi_{k-1}^{0,h}q \in Q_h$ is defined by

$$ (q - \Pi_{k-1}^{0,h}q, w)_K = 0, \quad \forall w \in \mathbb{P}_{k-1}(K), \forall K \in \mathcal{T}_h. $$

Since VE functions are not explicitly known, VEMs require carefully designed computable discrete bilinear forms to ensure the accuracy and reliability of the discrete scheme. To this end, for any polytope $K \in \mathcal{T}_h$, we assume that the computable local discrete bilinear form $a_h^K(\cdot, \cdot) : W_h|_K \times W_h|_K \to \mathbb{R}$ has the following properties [4]

* k-consistency:

$$ a_h^K(\mathbf{q}_k, \mathbf{w}_h) = a^K(\mathbf{q}_k, \mathbf{w}_h), \quad \forall \mathbf{q}_k \in \mathbb{P}_k(K; \mathbb{R}^d), \forall \mathbf{w}_h \in W_h; \quad (3.2) $$

* stability:

$$ \alpha_* a^K(\mathbf{w}_h, \mathbf{w}_h) \le a_h^K(\mathbf{w}_h, \mathbf{w}_h) \le \alpha^* a^K(\mathbf{w}_h, \mathbf{w}_h), \quad \forall \mathbf{w}_h \in W_h. \quad (3.3) $$

Then, the global discrete bilinear form $a_h(\cdot, \cdot) : W_h \times W_h \to \mathbb{R}$ is defined by

$$ a_h(v_h, w_h) = \sum_{K \in \mathcal{T}_h} a_h^K(v_h, w_h), \quad \forall v_h, w_h \in W_h. $$

Another computable bilinear form $b_h(\cdot, \cdot)$ can be simply defined as follows

$$ b_h(\mathbf{w}_h, q_h) = -(\Pi_{k-1}^{0,h} \text{div}_h \mathbf{w}_h, q_h) = -(\text{div}_h \mathbf{w}_h, q_h), \quad \forall \mathbf{w}_h \in W_h, \forall q_h \in Q_h. $$
For the right-hand side, a common and natural choice is to approximate $(f, w_h)$ by $(f, \Pi_k^{0,h} w_h)$. We write the common VE discretization of problem (1.1): find $(u_h, p_h) \in W_h \times Q_h$ such that

$$ \left\{ \begin{array}{ll} a_h(u_h, w_h) + b_h(w_h, p_h) = (f, \Pi_k^{0,h} w_h), & \forall w_h \in W_h, \\ b_h(u_h, q_h) = 0, & \forall q_h \in Q_h. \end{array} \right. \qquad (3.4) $$

However, this choice breaks the pressure robustness of the VEM, even when using an $H^1$-conforming VE with exactly divergence-free property [14]. The standard error estimates for the discrete scheme (3.4) have been well established for the conforming and nonconforming VEMs in [7, 14, 22, 26, 27].

**Theorem 1.** Let $f \in H^{k+1}(\Omega; \mathbb{R}^d)$, $(u, p) \in (H_0^1(\Omega; \mathbb{R}^d) \cap H^{k+1}(\Omega; \mathbb{R}^d)) \times (L_0^2(\Omega) \cap H^k(\Omega))$ and $(u_h, p_h) \in W_h \times Q_h$ be the solutions of problem (1.1) and problem (3.4), respectively. Then we have

$$ \begin{gathered} |u - u_h|_{1,h} \lesssim h^k |u|_{k+1} + \frac{1}{\nu} (h^{k+2} |f|_{k+1} + E_h), \\ \|p - p_h\|_0 \lesssim h^k |p|_k + \nu h^k |u|_{k+1} + h^{k+2} |f|_{k+1} + E_h, \end{gathered} $$

where the consistency error is defined by

$$ E_h := \sup_{0 \neq w_h \in W_h} \frac{\nu(\nabla_h u, \nabla_h w_h) - (\operatorname{div}_h w_h, p) - (f, w_h)}{|w_h|_{1,h}}. $$

We remark that, it holds $E_h = 0$ for the conforming VE ($W_h \subset H_0^1(\Omega; \mathbb{R}^d)$). For the nonconforming VE, the weak continuity (3.1) implies that $E_h$ can be bounded by the optimal rate $h^k$ but the resulting upper bound depends on the pressure $p$.

### 3.3. The pressure-robust virtual element method based on reconstruction

Motivated by the reconstruction idea from the pressure-robust FEMs [16, 20], a promising approach is to employ a reconstruction operator to restore pressure robustness within the VE framework. Unfortunately, as pointed out in [14], the reconstruction operators developed for FEMs are not directly applicable to general VEMs. A brief explanation is provided below to clarify this fact.

The Raviart-Thomas (RT) FE space on a subtriangulation $\mathcal{T}_K$ of polytope $K$ is defined by

$$ W_{k-1}^{\mathrm{RT}}(K) = \{w \in \mathcal{H}(\mathrm{div}; K); w|_T \in \mathbb{P}_{k-1}(T; \mathbb{R}^d) + x\bar{\mathbb{P}}_{k-1}(T), \forall T \in \mathcal{T}_K\}. $$

For any function $w_h \in W_h|_K$, we define the standard RT interpolation $\Pi_K^{\mathrm{RT}} w_h \in W_{k-1}^{\mathrm{RT}}(K)$ such that

$$ \bullet \ (\Pi_K^{\mathrm{RT}} w_h \cdot n_F, q_{k-1})_F = (w_h \cdot n_F, q_{k-1})_F, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(F), \forall F \in \mathcal{F}(\mathcal{T}_K), $$

$$ \bullet \ (\Pi_K^{\mathrm{RT}} w_h, q_{k-2})_T = (w_h, q_{k-2})_T, \quad \forall q_{k-2} \in \mathbb{P}_{k-2}(T; \mathbb{R}^d), \forall T \in \mathcal{T}_K. $$

The following properties of the standard RT interpolation operator are well known:

$$ \bullet \ \mathrm{div} \Pi_K^{\mathrm{RT}} w_h = \Pi_{k-1}^{0, \mathcal{T}_K} \mathrm{div} \ w_h, \quad \text{in } T \in \mathcal{T}_K, \qquad (3.5) $$

$$ \bullet \ \|w_h - \Pi_K^{\mathrm{RT}} w_h\|_{0,T} \lesssim h_T |w_h|_{1,T}, \quad \forall T \in \mathcal{T}_K, \qquad (3.6) $$

$$ \bullet \ (w_h - \Pi_K^{\mathrm{RT}} w_h, q_{k-2})_T = 0, \quad \forall q_{k-2} \in \mathbb{P}_{k-2}(T; \mathbb{R}^d), \forall T \in \mathcal{T}_K, \qquad (3.7) $$

where $\Pi_{k-1}^{0, \mathcal{T}_K}|_T = \Pi_{k-1}^{0, T}$ is the $L^2$-projection on $\mathbb{P}_{k-1}(T)$ for $T \in \mathcal{T}_K$.

For a Stokes VE function $w_h$, we do not know its closed form, although some limited information

$$ (w_h \cdot n_F, q_{k-1})_F, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(F), \forall F \in \mathcal{F}(K), $$

and

$$ (w_h, q_{k-2})_K, \quad \forall q_{k-2} \in \mathbb{P}_{k-2}(K; \mathbb{R}^d), $$

can be obtained by the DoFs of $w_h$. Generally, it will appear

$$ \mathcal{F}(\mathcal{T}_K) \neq \mathcal{F}(K), $$

which leads to the incompatibility of $\Pi_K^{\mathrm{RT}} w_h$. Therefore, it is necessary to construct a computable interpolation operator with the essential properties of the RT interpolation.
**Remark 1.** The RT interpolation $\Pi^{\text{RT}}$ is usually used to obtain the pressure-robust FEM by simply replacing the usual right-hand side term $(f, w_h)$ by $(f, \Pi^{\text{RT}} w_h)$ in the discontinuous pressure approximation case, and $(f, R_h^{\text{FEM}} v_h)$ in the continuous pressure approximation case, where $R_h^{\text{FEM}}$ is some reconstruction operator based on the RT interpolation and some local problem on vertex patches (refer to [17] for the details). However, this technique cannot be directly applied to the design of pressure-robust VEM. The reason is that the standard RT interpolation is not computable for VEM, as mentioned above.

In this paper, we aim to propose a computable reconstruction operator $R_h$ with the following abstract properties

* commutativity with $\Pi_{k-1}^{0,h}$:

$$ \mathrm{div} R_h w_h = \Pi_{k-1}^{0,h} \mathrm{div}_h w_h, \quad \forall w_h \in W_h, \qquad (3.8) $$

* approximation error estimate:

$$ \|w_h - R_h w_h\|_{0,K} \lesssim h_K |w_h|_{1,K}, \quad \forall w_h \in W_h, \forall K \in \mathcal{T}_h, \qquad (3.9) $$

* polynomial $L^2$-orthogonality:

$$ (w_h - R_h w_h, q_{k-2})_K = 0, \quad \forall w_h \in W_h, \forall q_{k-2} \in \mathbb{P}_{k-2}(K; \mathbb{R}^d), \forall K \in \mathcal{T}_h. \qquad (3.10) $$

Here, $W_h$ is a VE space for velocity.

**Remark 2.** We note that the properties (3.8)-(3.10) of $R_h$ hold on the whole polytope $K$, not on each simplex in the simplicial subdivision $\mathcal{T}_K$, which enables us to design a computable reconstruction operator $R_h$.

Then, with $R_h$ we give the pressure-robust scheme for VEMs of problem (1.1): find $(u_h, p_h) \in W_h \times Q_h$ such that

$$ \left\{ \begin{array}{ll} a_h(u_h, w_h) + b_h(w_h, p_h) = (f, R_h w_h), & \forall w_h \in W_h, \\ b_h(u_h, q_h) = 0, & \forall q_h \in Q_h. \end{array} \right. \qquad (3.11) $$

Observing the stability (3.3), it is not difficult to see that $a_h(\cdot, \cdot)$ is continuous with respect to the norm $|\cdot|_{1,h}$ and coercive

$$ a_h(w_h, w_h) \ge \alpha_* \nu |w_h|_{1,h}^2, \quad \forall w_h \in W_h. \qquad (3.12) $$

We further suppose that the discrete space pair $W_h \times Q_h$ satisfies a discrete inf-sup condition

$$ \sup_{0 \neq w_h \in W_h} \frac{b_h(w_h, q_h)}{|w_h|_{1,h}} \ge \beta_* \|q_h\|_0, \quad \forall q_h \in Q_h. \qquad (3.13) $$

Under these conditions, we obtain the existence and uniqueness of the solution to the pressure-robust VEM. Moreover, the discrete inf-sup condition (3.13), together with the classical analysis in [15], yields the following lemma.

**Lemma 1.** Let the VE spaces $W_h$ and $Q_h$ fulfill the discrete inf-sup condition (3.13), then it holds that for all $w \in \mathbb{Z}$

$$ \inf_{w_h \in \mathbb{Z}_h} |w - w_h|_{1,h} \lesssim \inf_{v_h \in W_h} |w - v_h|_{1,h}, $$

where

$$ \mathbb{Z}_h := \{w_h \in W_h; b_h(w_h, q_h) = 0, \forall q_h \in Q_h\}. $$

### 3.4. Error estimates for the pressure-robust scheme

We first recall the following standard approximation result from piecewise polynomials, refer to [8].

**Lemma 2.** For every function $w \in H^s(\Omega; \mathbb{R}^d)$ with $0 \le m \le s \le k + 1$, there exists a piecewise k-degree polynomial $w_\pi$ such that

$$ |w - w_\pi|_{m,K} \lesssim h_K^{s-m} |w|_{s,K}, \quad \forall K \in \mathcal{T}_h. \qquad (3.14) $$

Next, we assume that the following interpolation estimates hold.
**Lemma 3.** Let $K \in \mathcal{T}_h$ and $\mathbf{w} \in H^s(K; \mathbb{R}^d)$ with $1 \le s \le k + 1$. Then there exists $\mathbf{w}_I \in \mathbf{W}_h$ such that

$$ \|\mathbf{w} - \mathbf{w}_I\|_{m,K} \lesssim h_K^{s-m} |\mathbf{w}|_{s,K}, \quad m = 0, 1. $$

Then, using the abstract properties (3.8)-(3.10) of $\mathbb{R}_h$, we establish the following error estimates for the VE approximation from the pressure-robust scheme (3.11).

**Theorem 2.** Let $(\mathbf{u}, p) \in (H_0^1(\Omega; \mathbb{R}^d) \cap H^{k+1}(\Omega; \mathbb{R}^d)) \times (L_0^2(\Omega) \cap H^k(\Omega))$ and $(\mathbf{u}_h, p_h) \in \mathbf{W}_h \times \mathcal{Q}_h$ be the solutions of problem (1.1) and problem (3.11), respectively. Then it holds that

$$ |\mathbf{u} - \mathbf{u}_h|_{1,h} \lesssim h^k |\mathbf{u}|_{k+1} + N_h, \quad (3.15) $$

$$ \|p - p_h\|_0 \lesssim h^k |p|_k + \nu (h^k |\mathbf{u}|_{k+1} + N_h), \quad (3.16) $$

where the nonconforming error is defined by

$$ N_h := \sup_{\mathbf{0} \neq \mathbf{w}_h \in \mathbf{W}_h} \frac{\sum_{F \in \mathcal{F}_h} \left( \frac{\partial \mathbf{u}}{\partial \mathbf{n}_F}, [\mathbf{w}_h] \right)_F}{|\mathbf{w}_h|_{1,h}}. $$

*Proof.* For an arbitrary $\mathbf{v}_h \in \mathcal{Z}_h$, let $\eta_h := \mathbf{v}_h - \mathbf{u}_h \in \mathcal{Z}_h$, then $\mathbb{R}_h \eta_h \in \mathcal{H}_0(\text{div}; \Omega)$ satisfies

$$ \text{div} \, \mathbb{R}_h \eta_h = \Pi_{k-1}^{0,h} \text{div}_h \eta_h = 0, \quad (3.17) $$

due to the commutativity (3.8). Combining (3.17) and the first equation in (3.11) leads to

$$ a_h(\mathbf{u}_h, \eta_h) = (\mathbf{f}, \mathbb{R}_h \eta_h) = (-\nu \Delta \mathbf{u} + \nabla p, \mathbb{R}_h \eta_h) = -(\nu \Delta \mathbf{u}, \mathbb{R}_h \eta_h). $$

Then, we can conclude from the coercivity (3.12) and k-consistency (3.2) of $a_h(\cdot, \cdot)$ that

$$ \begin{aligned} \alpha_* \nu |\eta_h|_{1,h}^2 &\le a_h(\eta_h, \eta_h) \\ &= a_h(\mathbf{v}_h, \eta_h) - a_h(\mathbf{u}_h, \eta_h) \\ &= \sum_{K \in \mathcal{T}_h} (a_h^K(\mathbf{v}_h - \mathbf{u}_\pi, \eta_h) + a_h^K(\mathbf{u}_\pi, \eta_h)) + (\nu \Delta \mathbf{u}, \mathbb{R}_h \eta_h) \\ &= \sum_{K \in \mathcal{T}_h} (a_h^K(\mathbf{v}_h - \mathbf{u}_\pi, \eta_h) + a_h^K(\mathbf{u}_\pi - \mathbf{u}, \eta_h)) + \left( \sum_{K \in \mathcal{T}_h} a_h^K(\mathbf{u}, \eta_h) + (\nu \Delta \mathbf{u}, \mathbb{R}_h \eta_h) \right), \end{aligned} \quad (3.18) $$

where $\mathbf{u}_\pi$ is a piecewise k-degree polynomial approximation to $\mathbf{u}$ with the estimates (3.14).

Using the stability condition (3.3) for $a_h(\cdot, \cdot)$, we have

$$ \left| \sum_{K \in \mathcal{T}_h} (a_h^K(\mathbf{v}_h - \mathbf{u}_\pi, \eta_h) + a^K(\mathbf{u}_\pi - \mathbf{u}, \eta_h)) \right| \lesssim \nu \left( |\mathbf{v}_h - \mathbf{u}_\pi|_{1,h} + |\mathbf{u}_\pi - \mathbf{u}|_{1,h} \right) |\eta_h|_{1,h}. \quad (3.19) $$

From the integration by parts, it follows that

$$ \begin{aligned} \frac{1}{\nu} \left( \sum_{K \in \mathcal{T}_h} a^K(\mathbf{u}, \eta_h) + (\nu \Delta \mathbf{u}, \mathbb{R}_h \eta_h) \right) &= \sum_{K \in \mathcal{T}_h} \left( \frac{\partial \mathbf{u}}{\partial \mathbf{n}_K}, \eta_h \right)_{\partial K} - (\Delta \mathbf{u}, \eta_h - \mathbb{R}_h \eta_h) \\ &= \sum_{F \in \mathcal{F}_h} \left( \frac{\partial \mathbf{u}}{\partial \mathbf{n}_F}, [\eta_h] \right)_F - (\Delta \mathbf{u}, \eta_h - \mathbb{R}_h \eta_h). \end{aligned} \quad (3.20) $$

Moreover, using the polynomial $\mathcal{L}^2$-orthogonality (3.10), the approximation error estimate (3.9) and the projection error estimates, we obtain

$$ (\Delta \mathbf{u}, \eta_h - \mathbb{R}_h \eta_h) = (\Delta \mathbf{u} - \Pi_{k-2}^{0,h} \Delta \mathbf{u}, \eta_h - \mathbb{R}_h \eta_h) \lesssim h^k |\mathbf{u}|_{k+1} |\eta_h|_{1,h}. \quad (3.21) $$
Combining (3.18)-(3.21), we derive

$$|\eta_h|_{1,h} \lesssim |v_h - u_\pi|_{1,h} + |u_\pi - u|_{1,h} + N_h + h^k|u|_{k+1}.$$

It follows from the triangle inequality and Lemma 2 that

$$|u - u_h|_{1,h} \le |u - v_h|_{1,h} + |\eta_h|_{1,h} \lesssim |u - v_h|_{1,h} + h^k|u|_{k+1} + N_h.$$

Therefore, we have

$$|u - u_h|_{1,h} \lesssim \inf_{v_h \in \mathcal{Z}_h} |u - v_h|_{1,h} + h^k|u|_{k+1} + N_h.$$

Applying Lemma 1 and Lemma 3 yields the velocity error estimate (3.15).

Next, we consider the pressure error. For any $w_h \in W_h$, the property (3.8) and the k-consistency (3.2) of $a_h(\cdot, \cdot)$ lead to

$$
\begin{align}
b_h(w_h, \Pi_{k-1}^{0,h} p) - b_h(w_h, p_h) & = -(\mathrm{div}_h w_h, \Pi_{k-1}^{0,h} p) + a_h(u_h, w_h) - (f, R_h w_h) \nonumber \\
& = (\mathrm{div}_h w_h, p) - (\mathrm{div}_h w_h, \Pi_{k-1}^{0,h} p) + a_h(u_h, w_h) + (\nu \Delta u, R_h w_h) \nonumber \\
& = \sum_{K \in \mathcal{T}_h} a_h^K (u_h - u_\pi, w_h) - \sum_{K \in \mathcal{T}_h} a^K (u - u_\pi, w_h) + \sum_{K \in \mathcal{T}_h} a^K (u, w_h) + (\nu \Delta u, R_h w_h). \tag{3.22}
\end{align}
$$

The stability condition (3.3) of $a_h(\cdot, \cdot)$ leads to

$$\left| \sum_{K \in \mathcal{T}_h} a_h^K (u_h - u_\pi, w_h) - \sum_{K \in \mathcal{T}_h} a^K (u - u_\pi, w_h) \right| \lesssim \nu \left( |u - u_h|_{1,h} + |u - u_\pi|_{1,h} \right) |w_h|_{1,h}. \quad (3.23)$$

From the proofs of (3.20) and (3.21), it follows that

$$\sum_{K \in \mathcal{T}_h} a^K(u, w_h) + (\nu \Delta u, R_h w_h) \lesssim \nu \sum_{F \in \mathcal{F}_h} \left( \frac{\partial u}{\partial n_F}, [w_h] \right)_F + \nu h^k |u|_{k+1} |w_h|_{1,h}. \quad (3.24)$$

Combining (3.22)-(3.24) and Lemma 2 leads to

$$b_h(w_h, \Pi_{k-1}^{0,h} p) - b_h(w_h, p_h) \lesssim \nu \left( |u - u_h|_{1,h} + h^k|u|_{k+1} + N_h \right) |w_h|_{1,h}.$$

By using the discrete inf-sup condition (3.13), we obtain

$$\|\Pi_{k-1}^{0,h} p - p_h\|_0 \le \frac{1}{\beta_*} \sup_{0 \ne w_h \in W_h} \frac{b_h(w_h, \Pi_{k-1}^{0,h} p) - b_h(w_h, p_h)}{|w_h|_{1,h}} \lesssim \nu \left( |u - u_h|_{1,h} + h^k |u|_{k+1} + N_h \right).$$

Thus, (3.16) is obtained by the triangle inequality and the estimate (3.15) of $|u - u_h|_{1,h}$. The proof is now complete. □

**Remark 3.** For the conforming VE, it is easy to see that $N_h = 0$ and the velocity error estimates is independent of the pressure $p$. For the nonconforming VE, the weak continuity (3.1) implies that $N_h$ can be bounded by the optimal rate $h^k$ and the resulting upper bound is independent of the pressure $p$. The VE scheme (3.11) is pressure-robust.

# 4. Reconstruction

## 4.1. Interpolation for VE to $H(\mathrm{div})$-conforming FE

In order to define a proper reconstruction operator, we first introduce a $H(\mathrm{div})$-conforming macro FE space. Given a $d$-dimensional polytope $K \in \mathcal{T}_h$ with a quasi-uniform simplicial partition $\mathcal{T}_K$, let $W_h^{d-2}(K)$ be the FE space with degree $k$ defined over $\mathcal{T}_K$, taking the form of the Lagrange element space for $d=2$ and the
first kind of Nédélec element space for $d = 3$. This space plays a key role in the construction of the discrete de Rham complex. Furthermore we define the space

$$ \dot{W}_k^{d-2}(K) := \{w \in W_k^{d-2}(K); \text{tr } w = 0 \text{ on } \partial K\}, $$

where tr denotes the trace operator for $H^1(K)$ when $d = 2$, and the trace operator for $H(\text{curl}; K)$ when $d = 3$. In addition, we define the space

$$ W_{k-1}^{L^2}(K) := \{w \in L^2(K); w|_T \in \mathbb{P}_{k-1}(T), \forall T \in \mathcal{T}_K\}. $$

Set $\dot{W}_{k-1}^{L^2}(K) = W_{k-1}^{L^2}(K) \cap L_0^2(K)$ and $\dot{W}_{k-1}^{\text{RT}}(K) = W_{k-1}^{\text{RT}}(K) \cap H_0(\text{div}; K)$. We recall the local finite element de Rham complexes presented in [3, 11].

**Lemma 4.** *Finite element complexes*

$$ W_k^{d-2}(K) \xrightarrow{\text{curl}} W_{k-1}^{\text{RT}}(K) \xrightarrow{\text{div}} W_{k-1}^{L^2}(K) \to 0, $$

$$ \dot{W}_k^{d-2}(K) \xrightarrow{\text{curl}} \dot{W}_{k-1}^{\text{RT}}(K) \xrightarrow{\text{div}} \dot{W}_{k-1}^{L^2}(K) \to 0, $$

are exact.

Based on the result of Lemma 4, we define the space

$$ \dot{D}_{k-1}^{\text{RT}}(K) := \{w \in \dot{W}_{k-1}^{\text{RT}}(K); \text{div } w = 0\} = \text{curl } \dot{W}_k^{d-2}(K). $$

In order to construct a reconstruction operator that satisfies property (3.8), we introduce the space $(k \ge 1)$

$$ \widetilde{W}_{k-1}^{\text{div}}(K) = \{w \in W_{k-1}^{\text{RT}}(K); \text{div } w \in \mathbb{P}_{k-1}(K)\}, $$

along with the following DoFs

$$ \bullet \text{ the moments } \frac{1}{|F|}(w \cdot n_F, q_{k-1})_F, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(F), \forall F \in \mathcal{F}^\partial(\mathcal{T}_K), \qquad (4.1) $$

$$ \bullet \text{ the moments } \frac{1}{|K|}(\text{div } w, q_{k-1})_K, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(K)/\mathbb{R}, \qquad (4.2) $$

$$ \bullet \text{ the moments } \frac{1}{|K|}(w, q)_K, \quad \forall q \in \dot{D}_{k-1}^{\text{RT}}(K). \qquad (4.3) $$

The following lemma was established in [11].

**Lemma 5.** *The space $\widetilde{W}_{k-1}^{\text{div}}(K)$ can be decomposed as a direct sum as follows*

$$ \widetilde{W}_{k-1}^{\text{div}}(K) = \text{curl } W_k^{d-2}(K) \oplus (x - x_K)\mathbb{P}_{k-1}(K). $$

Then the set of DoFs (4.1)-(4.3) is unisolvent for space $\widetilde{W}_{k-1}^{\text{div}}(K)$.

It should be noted that DoFs (4.1) depend on the face $F \in \mathcal{F}^\partial(\mathcal{T}_K)$, while VEs generally only provide information associated with $F \in \mathcal{F}(K)$. To ensure the computability of the reconstruction operator, it is necessary to impose additional constraints on the space $\widetilde{W}_{k-1}^{\text{div}}(K)$. Following [11], we now introduce a $H(\text{div})$-conforming macro space

$$ W_{k-1}^{\text{div}}(K) = \{w \in \widetilde{W}_{k-1}^{\text{div}}(K); w \cdot n_F|_F \in \mathbb{P}_{k-1}(F), \forall F \in \mathcal{F}(K)\}, $$

equipped with the following DoFs

$$ \bullet \text{ the moments } \frac{1}{|F|}(w \cdot n_F, q_{k-1})_F, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(F), \forall F \in \mathcal{F}(K), \qquad (4.4) $$

$$ \bullet \text{ the moments } \frac{1}{|K|}(\text{div } w, q_{k-1})_K, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(K)/\mathbb{R}, \qquad (4.5) $$

$$ \bullet \text{ the moments } \frac{1}{|K|}(w, q)_K, \quad \forall q \in \dot{D}_{k-1}^{\text{RT}}(K). \qquad (4.6) $$
It follows from Lemma 5 that DoFs (4.4)-(4.6) are unisolvent for space $\tilde{W}_{k-1}^{\text{div}}(K)$.

For any VE function $w_h \in W_h|_K \subset H^1(K; \mathbb{R}^d)$, we define its local interpolation function $I_K^{\text{div}} w_h \in W_{k-1}^{\text{div}}(K)$ such that

$$(I_K^{\text{div}} w_h \cdot n_F, q_{k-1})_F = (w_h \cdot n_F, q_{k-1})_F, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(F), \forall F \in \mathcal{F}(K), \quad (4.7)$$

$$(\text{div} I_K^{\text{div}} w_h, q_{k-1})_K = (\text{div} w_h, q_{k-1})_K, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(K)/\mathbb{R}, \quad (4.8)$$

$$(I_K^{\text{div}} w_h, q)_K = (\Pi_k^{0,K} w_h, q)_K, \quad \forall q \in \mathring{D}_{k-1}^{\text{RT}}(K). \quad (4.9)$$

Then, the global $H(\text{div})$-conforming FE space is given by

$$W_{k-1}^{\text{div}} = \{w \in H_0(\text{div}; \Omega); w|_K \in W_{k-1}^{\text{div}}(K), \forall K \in \mathcal{T}_h\},$$

and the global interpolation $I_h^{\text{div}}: W_h \to W_{k-1}^{\text{div}}$ by setting $I_h^{\text{div}}|_K = I_K^{\text{div}}$, $K \in \mathcal{T}_h$.

**Remark 4.** For VE function $w_h \in W_h$, we require that the interpolation $I_K^{\text{div}} w_h$ are computable from the VE DoFs.

We next present some properties of the space $W_{k-1}^{\text{div}}$. For any $w_h \in W_h$, noting that

$$(\text{div} w_h, q_0)_K = (w_h \cdot n_K, q_0)_{\partial K}, \quad \forall q_0 \in \mathbb{P}_0(K),$$

and using the definitions of DoFs (4.4)-(4.5), we obtain

$$(\text{div} w_h, q_{k-1})_K = (\text{div} I_K^{\text{div}} w_h, q_{k-1})_K, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(K).$$

Therefore, it holds that

$$(\Pi_{k-1}^{0,h} \text{div}_h w_h, q_h) = (\text{div}_h w_h, q_h) = \sum_{K \in \mathcal{T}_h} (\text{div} I_K^{\text{div}} w_h, q_h)_K = (\text{div} I_h^{\text{div}} w_h, q_h), \quad \forall q_h \in Q_h,$$

which implies the commutativity with $\Pi_{k-1}^{0,h}$

$$\text{div} I_h^{\text{div}} w_h = \Pi_{k-1}^{0,h} \text{div}_h w_h. \quad (4.10)$$

**Remark 5.** Note that the $H(\text{div})$ macro FE here is a little different from the origin one from the reference [11] where the divergence of macro FE is a polynomial of degree $(k-2)$. This is to ensure the commutativity (3.8) with $\Pi_{k-1}^{0,h}$.

We now turn our attention to the interpolation error estimate for $I_h^{\text{div}}$. As discussed in [11], we can derive the following lemma.

**Lemma 6.** For $w \in W_{k-1}^{\text{div}}(K)$, it holds that

$$\|w\|_{0,K} \lesssim h_K \|\text{div } w\|_{0,K} + \sup_{0 \neq q \in \mathring{D}_{k-1}^{\text{RT}}(K)} \frac{(w, q)_K}{\|q\|_{0,K}} + \sum_{F \in \mathcal{F}(K)} h_F^{1/2} \|w \cdot n_F\|_{0,F}.$$

With the help of the above lemma, we can derive the following interpolation error estimate.

**Lemma 7.** Let $w \in H^1(K; \mathbb{R}^d)$. Then it holds that

$$\|w - I_K^{\text{div}} w\|_{0,K} \lesssim h_K |w|_{1,K}. \quad (4.11)$$

*Proof.* Let $w_0 = \Pi_0^{0,K} w$, we have $w_0 - I_K^{\text{div}} w \in W_{k-1}^{\text{div}}(K)$. It follows from the commutativity (4.10) that

$$\| \text{div}(w_0 - I_K^{\text{div}} w) \|_{0,K} = \| \Pi_{k-1}^{0,K} \text{div} w \|_{0,K} \lesssim |w|_{1,K}. \quad (4.12)$$

By the definition of $I_K^{\text{div}}$, we can deduce that

$$(I_K^{\text{div}} w, q)_K = (\Pi_k^{0,K} w, q)_K, \quad \forall q \in \mathring{D}_{k-1}^{\text{RT}}(K).$$
Then, apply the Cauchy-Schwarz inequality and the boundedness of $\Pi_k^{0,K}$ to acquire

$$ \sup_{0 \neq q \in \bar{D}_{k-1}^{\text{RT}}(K)} \frac{(w_0 - \Pi_K^{\text{div}} w, q)_K}{\|q\|_{0,K}} = \sup_{0 \neq q \in \bar{D}_{k-1}^{\text{RT}}(K)} \frac{(w_0 - \Pi_k^{0,K} w, q)_K}{\|q\|_{0,K}} \le \|\Pi_k^{0,K}(w_0 - w)\|_{0,K} \lesssim h_K |w|_{1,K}. \quad (4.13) $$

Recalling $w_0 - \Pi_K^{\text{div}} w \in W_{k-1}^{\text{div}}(K)$ and the definition of $\Pi_K^{\text{div}}$, it holds that

$$ \|(w_0 - \Pi_K^{\text{div}} w) \cdot n_F\|_{0,F}^2 = ((w_0 - w) \cdot n_F, (w_0 - \Pi_K^{\text{div}} w) \cdot n_F)_F, $$

which, together with the Cauchy-Schwarz inequality, the trace inequality and standard projection error estimates, yields

$$ \begin{aligned} h_K^{1/2} \|(w_0 - \Pi_K^{\text{div}} w) \cdot n_F\|_{0,F} &\le h_K^{1/2} \|(w_0 - w) \cdot n_F\|_{0,F} \\ &\lesssim \|w_0 - w\|_{0,K} + h_K |w_0 - w|_{1,K} \\ &\lesssim h_K |w|_{1,K}. \end{aligned} \quad (4.14) $$

Therefore, recalling the triangle inequality

$$ \|w - \Pi_K^{\text{div}} w\|_{0,K} \le \|w - w_0\|_{0,K} + \|w_0 - \Pi_K^{\text{div}} w\|_{0,K}, $$

and the fact that $w_0 - \Pi_K^{\text{div}} w \in W_{k-1}^{\text{div}}(K)$, the estimate (4.11) follows from Lemma 6 and (4.12)-(4.14). $\square$

By applying the Green's formula

$$ (w_h, \nabla q_{k-1})_K = (w_h \cdot n_K, q_{k-1})_{\partial K} - (\text{div} w_h, q_{k-1})_K, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(K), $$

from the definition of $\Pi_K^{\text{div}}$ we obtain

$$ (w_h, \nabla q_{k-1})_K = (\Pi_K^{\text{div}} w_h, \nabla q_{k-1})_K, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(K). \quad (4.15) $$

For low-order cases ($k=1,2$), since $G_{k-2}^{x*}(K) = \{0\}$, the decomposition

$$ \mathbb{P}_{k-2}(K; \mathbb{R}^d) = \nabla \mathbb{P}_{k-1}(K) \oplus G_{k-2}^{x*}(K), $$

together with the property (4.15), implies that only for $k=1,2$,

$$ (w_h, q_{k-2})_K = (\Pi_K^{\text{div}} w_h, q_{k-2})_K, \quad \forall q_{k-2} \in \mathbb{P}_{k-2}(K; \mathbb{R}^d). $$

Therefore, for $k=1,2$, the interpolation operator $\Pi_K^{\text{div}}$ satisfies the polynomial $L^2$-orthogonality (3.10). However, for higher-order cases ($k \ge 3$), the interpolation operator $\Pi_h^{\text{div}}$ fails to satisfy the required polynomial $L^2$-orthogonality (3.10). Since the polynomial $L^2$-orthogonality (3.10) plays a critical role in achieving optimal convergence rate for the pressure-robust VEMs (3.11), next we will amend the interpolation by a solution of a specially designed local problem to ensure the property (3.10) holds.

## 4.2. A specifically designed local problem

Motivating by the idea proposed in [17], for $k \ge 3$ we define a local problem on each element $K \in \mathcal{T}_h$ with subtriangulation $\mathcal{T}_K$. In the quasi-uniform simplicial partition $\mathcal{T}_K$ of the polytope $K$, we arbitrarily select an interior $(d-1)$-face $F = \partial T_1 \cap \partial T_2$, where $T_1, T_2 \in \mathcal{T}_K$. We note that $F$ is a $(d-1)$-dimensional simplex. Let $x_F$ be the barycenter of $F$. Then we have the following decomposition of the polynomial space

$$ \mathbb{P}_{k-2}(K; \mathbb{R}^d) = \nabla \mathbb{P}_{k-1}(K) \oplus G_{k-2}^{x_F}(K), \quad (4.16) $$

where the subspace

$$ G_{k-2}^{x_F}(K) = \begin{cases} (x - x_F)^\perp \mathbb{P}_{k-3}(K), & d=2, \\ (x - x_F) \times \mathbb{P}_{k-3}(K; \mathbb{R}^3), & d=3. \end{cases} $$
We define the local problem: find $(\sigma_h^K, \lambda_h^K) \in \mathring{D}_{k-1}^{\text{RT}}(K) \times \mathbf{G}_{k-2}^{\mathbf{x}_F}(K)$ such that

$$ \left\{ \begin{array}{ll} a_1^K(\sigma_h^K, \tau) + b_1^K(\tau, \lambda_h^K) = 0, & \forall \tau \in \mathring{D}_{k-1}^{\text{RT}}(K), \\ b_1^K(\sigma_h^K, \mu) = (w_h - I_h^{\text{div}} w_h, \mu)_K, & \forall \mu \in \mathbf{G}_{k-2}^{\mathbf{x}_F}(K), \end{array} \right. \qquad (4.17) $$

where the local bilinear forms are

$$ a_1^K(\sigma_h^K, \tau) := (\sigma_h^K, \tau)_K, \quad \forall (\sigma_h^K, \tau) \in \mathring{D}_{k-1}^{\text{RT}}(K) \times \mathring{D}_{k-1}^{\text{RT}}(K), $$

$$ b_1^K(\tau, \lambda_h^K) := (\tau, \lambda_h^K)_K, \quad \forall (\tau, \lambda_h^K) \in \mathring{D}_{k-1}^{\text{RT}}(K) \times \mathbf{G}_{k-2}^{\mathbf{x}_F}(K), $$

Before discussing the well-posedness of local problem (4.17), we prove a basic result for polynomials (also see [9, Section 4.3] for a more general result).

**Lemma 8.** Let $T_1$ and $T_2$ be the shape-regular simplices sharing one $(d-1)$-face. For any polynomial $q \in \mathbb{P}_k(K)$, the following norm equivalence holds

$$ \|q\|_{0,T_1} \approx \|q\|_{0,T_2}. \qquad (4.18) $$

Thus on a quasi-uniform subdivision $\mathcal{T}_K$ of the polytope $K$, we have

$$ \|q\|_{0,T} \approx \|q\|_{0,K}, \quad \forall T \in \mathcal{T}_K. $$

<figure><img src="image_5.png" /><figcaption>Figure 3: Shape-regular triangles and concentric discs.</figcaption></figure>

*Proof.* Since $T_1$ and $T_2$ are shape-regular, there exist two concentric discs/balls $\mathfrak{B}_1$ and $\mathfrak{B}_2$ with $\mathfrak{B}_1 \subset \mathfrak{B}_2$ (see Figure 3 for the case $d=2$), such that the ratio of their radii is less than a constant $C_r$ and

$$ \mathfrak{B}_1 \subset T_1 \subset T_1 \cup T_2 \subset \mathfrak{B}_2. $$

Then for any polynomial $q \in \mathbb{P}_k(K)$, we have

$$ \|q\|_{0,T_2} \le \|q\|_{0,T_1 \cup T_2} \le \|q\|_{0,\mathfrak{B}_2} \lesssim \|q\|_{0,\mathfrak{B}_1} \lesssim \|q\|_{0,T_1}. $$

Similar argument yields

$$ \|q\|_{0,T_1} \lesssim \|q\|_{0,T_2}. $$

This concludes the proof of norm equivalence (4.18). □

Since the identity

$$ \|\tau\|_{\text{div},K} = \|\tau\|_{0,K}, \quad \forall \tau \in \mathring{D}_{k-1}^{\text{RT}}(K), $$

holds trivially, we proceed to derive the following inf-sup condition for the bilinear form $b_1^K(\cdot, \cdot)$ in $L^2$-norm.
**Lemma 9.** For $k \ge 3$, the bilinear form $b_1^K(\cdot, \cdot)$ satisfies the discrete inf-sup condition

$$ \sup_{0 \neq \tau \in \mathring{D}_{k-1}^{\text{RT}}(K)} \frac{b_1^K(\tau, \mu)}{\|\tau\|_{0,K}} \geq \beta_* \|\mu\|_{0,K}, \quad \forall \mu \in \mathcal{G}_{k-2}^{\mathbf{x}_F}(K), \qquad (4.19) $$

where the constant $\beta_*$ is independent of the mesh size $h_K$, but depends on quasi-uniform constants of the virtual triangulation $\mathcal{T}_K$.

**Proof.** For the interior $(d-1)$-face $F$ shared by two simplices $T_1$ and $T_2$ in $\mathcal{T}_K$, by $\mathbf{a}_0, \mathbf{a}_1, \dots, \mathbf{a}_{d-1}$ we denote the vertices of $F$, which are also vertices of both $T_1$ and $T_2$. Assume $\lambda_0, \lambda_1, \dots, \lambda_{d-1}$ are the barycenter coordinates with respect to $\mathbf{a}_0, \mathbf{a}_1, \dots, \mathbf{a}_{d-1}$. Define the face bubble function

$$ b_F := \lambda_0 \lambda_1 \dots \lambda_{d-1} \in H_0^1(\omega_F), $$

where $\omega_F = T_1 \cup T_2$. This face bubble function can be naturally extended by zero to the whole polytope $K$, denoted still by $b_F$ for convenience.

Let $\mathbf{r} = \mathbf{x} - \mathbf{x}_F$ for brevity. We start with the case $d=2$. For any $\mu = \mathbf{r}^{\perp} \in \mathcal{G}_{k-2}^{\mathbf{x}_F}(K)$ with $\xi \in \mathbb{P}_{k-3}(K)$, we define

$$ \tau = -\operatorname{curl}(b_F \xi) \in \mathring{D}_{k+d-4}^{\text{RT}}(K) \subset \mathring{D}_{k-1}^{\text{RT}}(K). $$

Using integration by parts we obtain

$$ (-\operatorname{curl}(b_F \xi), \mathbf{r}^{\perp} \xi)_K = (-\nabla(b_F \xi), \mathbf{r} \xi)_{\omega_F} = (b_F \xi, \operatorname{div}(\mathbf{r} \xi))_{\omega_F} = (b_F \xi, 2\xi)_{\omega_F} + (b_F \xi, \mathbf{r} \cdot \nabla \xi)_{\omega_F}, $$

and so

$$ \begin{aligned} b_1^K(\tau, \mu) &= 2(b_F, \xi^2)_{\omega_F} + \frac{1}{2}(\mathbf{r}b_F, \nabla\xi^2)_{\omega_F} \\ &= 2(b_F, \xi^2)_{\omega_F} - \frac{1}{2}(2b_F, \xi^2)_{\omega_F} - \frac{1}{2}(\mathbf{r} \cdot \nabla b_F, \xi^2)_{\omega_F} \\ &= (b_F, \xi^2)_{\omega_F} - \frac{1}{2}(\mathbf{r} \cdot \nabla b_F, \xi^2)_{\omega_F}. \end{aligned} \qquad (4.20) $$

Recalling the identity $(\mathbf{x} - \mathbf{a}_i) \cdot \nabla \lambda_j = \lambda_j - \delta_{ij}$ ($\delta_{ij}$ is the usual Kronecker symbol), we have

$$ \begin{aligned} \mathbf{r} \cdot \nabla b_F &= \frac{1}{2}((\mathbf{x} - \mathbf{a}_0) \cdot \nabla(\lambda_0 \lambda_1) + (\mathbf{x} - \mathbf{a}_1) \cdot \nabla(\lambda_0 \lambda_1)) \\ &= \frac{1}{2}((2\lambda_0 - 1)\lambda_1 + \lambda_0(2\lambda_1 - 1)) \\ &= 2\lambda_0 \lambda_1 - \frac{1}{2}(\lambda_0 + \lambda_1), \end{aligned} \qquad (4.21) $$

which, together with (4.20) and Lemma 8, yields

$$ \begin{aligned} b_1^K(\tau, \mu) &= (b_F, \xi^2)_{\omega_F} - \frac{1}{2}(2b_F, \xi^2)_{\omega_F} + \frac{1}{4}(\lambda_0 + \lambda_1, \xi^2)_{\omega_F} \\ &= \frac{1}{4}(\lambda_0 + \lambda_1, \xi^2)_{\omega_F} \ge \|\xi\|_{0,\omega_F}^2 \ge \|\xi\|_{0,K}^2. \end{aligned} \qquad (4.22) $$

The standard argument by the elementwise affine mapping leads to

$$ \|\mu\|_{0,K} \lesssim h_K \|\xi\|_{0,K}. \qquad (4.23) $$

We combine the inequalities (4.22)-(4.23) to get

$$ \frac{b_1^K(\tau, \mu)}{\|\tau\|_{0,K}} \ge \frac{b_1^K(\tau, \mu)}{h_K^{-1}\|\xi\|_{0,K}} \ge h_K \|\xi\|_{0,K} \ge \beta_* \|\mu\|_{0,K}. $$

Thus the discrete inf-sup condition (4.19) holds for $d=2$.
Next, we consider the case $d = 3$. For any $\boldsymbol{\mu} = \boldsymbol{r} \times \boldsymbol{\xi} \in G_{k-2}^{\boldsymbol{\mu}_F}(K)$ with $\boldsymbol{\xi} \in \mathbb{P}_{k-3}(K; \mathbb{R}^d)$, from [17, Theorem 20] we can assume $\mathrm{div} \boldsymbol{\xi} = 0$. Then we define

$$ \boldsymbol{\tau} = - \mathrm{curl} (b_F \boldsymbol{\xi}) \in \mathring{D}_{k+d-4}^{\mathrm{RT}}(K) = \mathring{D}_{k-1}^{\mathrm{RT}}(K). $$

Employing the integration by parts we get

$$ (-\mathrm{curl}(b_F \boldsymbol{\xi}), \boldsymbol{r} \times \boldsymbol{\xi})_K = (b_F \boldsymbol{\xi}, -\mathrm{curl}(\boldsymbol{r} \times \boldsymbol{\xi}))_K = (b_F \boldsymbol{\xi}, -\mathrm{curl}(\boldsymbol{r} \times \boldsymbol{\xi}))_{\omega_F}, $$

which, together with the standard vector calculus

$$ \mathrm{curl}(\boldsymbol{r} \times \boldsymbol{\xi}) = (\nabla \boldsymbol{r})\boldsymbol{\xi} - (\nabla \boldsymbol{\xi})\boldsymbol{r} + \boldsymbol{r} \operatorname{div} \boldsymbol{\xi} - \boldsymbol{\xi} \operatorname{div} \boldsymbol{r} = -2\boldsymbol{\xi} - (\nabla \boldsymbol{\xi})\boldsymbol{r}, $$

implies that

$$
\begin{align}
b_1^K (\boldsymbol{\tau}, \boldsymbol{\mu}) &= 2(b_F \boldsymbol{\xi}, \boldsymbol{\xi})_{\omega_F} + (b_F \boldsymbol{\xi}, (\nabla \boldsymbol{\xi})\boldsymbol{r})_{\omega_F} \nonumber \\
&= 2(b_F, \boldsymbol{\xi}^2)_{\omega_F} + \frac{1}{2}(b_F \boldsymbol{r}, \nabla \boldsymbol{\xi}^2)_{\omega_F} \nonumber \\
&= 2(b_F, \boldsymbol{\xi}^2)_{\omega_F} - \frac{1}{2}(\mathrm{div}(b_F \boldsymbol{r}), \boldsymbol{\xi}^2)_{\omega_F} \nonumber \\
&= 2(b_F, \boldsymbol{\xi}^2)_{\omega_F} - \frac{1}{2}(\boldsymbol{r} \cdot \nabla b_F + 3b_F, \boldsymbol{\xi}^2)_{\omega_F} \nonumber \\
&= \frac{1}{2} \left[ (b_F, \boldsymbol{\xi}^2)_{\omega_F} - (\boldsymbol{r} \cdot \nabla b_F, \boldsymbol{\xi}^2)_{\omega_F} \right]. \tag{4.24}
\end{align}
$$

It follows from the properties of the barycentric coordinate $\lambda_i$ that

$$
\begin{align}
\boldsymbol{r} \cdot \nabla b_F &= \frac{1}{3} ((\boldsymbol{x} - \boldsymbol{a}_0) \cdot \nabla (\lambda_0 \lambda_1 \lambda_2) + (\boldsymbol{x} - \boldsymbol{a}_1) \cdot \nabla (\lambda_0 \lambda_1 \lambda_2) + (\boldsymbol{x} - \boldsymbol{a}_2) \cdot \nabla (\lambda_0 \lambda_1 \lambda_2)) \nonumber \\
&= \frac{1}{3} ((3\lambda_0 - 1)\lambda_1\lambda_2 + (3\lambda_1 - 1)\lambda_2\lambda_0 + (3\lambda_2 - 1)\lambda_0\lambda_1) \nonumber \\
&= 3b_F - \frac{1}{3}(\lambda_1\lambda_2 + \lambda_2\lambda_0 + \lambda_0\lambda_1). \tag{4.25}
\end{align}
$$

Noting that $0 \le \lambda_0, \lambda_1, \lambda_2 \le 1$ and $\lambda_0 + \lambda_1 + \lambda_2 \le 1$ on $\omega_F = T_1 \cup T_2$, we get

$$ b_F = \lambda_0 \lambda_1 \lambda_2 \le \frac{1}{27}, $$

thus we have

$$ \frac{1}{3}(\lambda_1 \lambda_2 + \lambda_2 \lambda_0 + \lambda_0 \lambda_1) \ge (\lambda_0 \lambda_1 \lambda_2)^{\frac{2}{3}} \ge 3b_F. \tag{4.26} $$

Combining (4.25)-(4.26) leads to

$$ \boldsymbol{r} \cdot \nabla b_F \le 0, \quad \text{on } \omega_F, $$

which, together with (4.24) and Lemma 8, yields

$$ b_1^K (\boldsymbol{\tau}, \boldsymbol{\mu}) \ge \frac{1}{2} (b_F, \boldsymbol{\xi}^2)_{\omega_F} \ge \|\boldsymbol{\xi}\|_{0, \omega_F}^2 \ge \|\boldsymbol{\xi}\|_{0, K}^2. $$

Similar to the case $d = 2$, we have

$$ \|\boldsymbol{\mu}\|_{0,K} \lesssim h_K \|\boldsymbol{\xi}\|_{0,K}. $$

Thus we obtain

$$ \frac{b_1^K (\boldsymbol{\tau}, \boldsymbol{\mu})}{\|\boldsymbol{\tau}\|_{0,K}} \ge \frac{b_1^K (\boldsymbol{\tau}, \boldsymbol{\mu})}{h_K^{-1} \|\boldsymbol{\xi}\|_{0,K}} \ge h_K \|\boldsymbol{\xi}\|_{0,K} \ge \beta_* \|\boldsymbol{\mu}\|_{0,K}, $$

which implies that the discrete inf-sup condition (4.19) holds for $d = 3$. $\square$
**Remark 6.** If the barycentric subdivision is employed for the element K (see Figures 1c and 2), then the second polynomial space can be chosen as $G_{k-2}^{x_K}(K)$ with $x_K$ instead of $x_F$. The corresponding inf-sup condition can be proved by using the technique proposed in [17]. However, the proof therein is difficult to generalize to more general subdivisions.

Furthermore, we have the ellipticity of $a_1^K(\cdot, \cdot)$

$$ \|\tau\|_{0,K}^2 = a_1^K(\tau, \tau), \quad \forall \tau \in \dot{D}_{k-1}^{\text{RT}}(K). $$

Therefore, in view of the classical Brezzi theory, the saddle-point problem (4.17) has a unique solution $(\sigma_h^K, \lambda_h^K) \in \dot{D}_{k-1}^{\text{RT}}(K) \times G_{k-2}^{x_F}(K)$ with the estimates

$$ \|\sigma_h^K\|_{0,K} + \|\lambda_h^K\|_{0,K} \lesssim \|w_h - I_h^{\text{div}} w_h\|_{0,K}. \quad (4.27) $$

We proceed to analyze the properties of the solution $\sigma_h^K$ to the local problem (4.17). Since $\sigma_h^K \in \dot{D}_{k-1}^{\text{RT}}(K)$, we have

$$ \operatorname{div} \sigma_h^K = 0, \quad (4.28) $$

and

$$ (\sigma_h^K, \nabla q_{k-1})_K = (\sigma_h^K \cdot n_K, q_{k-1})_{\partial K} - (\operatorname{div} \sigma_h^K, q_{k-1})_K = 0, \quad \forall q_{k-1} \in \mathbb{P}_{k-1}(K). $$

Noting that

$$ (\sigma_h^K, \mu)_K = (w_h - I_h^{\text{div}} w_h, \mu)_K, \quad \forall \mu \in G_{k-2}^{x_F}(K), $$

and the decomposition (4.16), from (4.15) we obtain

$$ (\sigma_h^K, q_{k-2})_K = (w_h - I_h^{\text{div}} w_h, q_{k-2})_K, \quad \forall q_{k-2} \in \mathbb{P}_{k-2}(K; \mathbb{R}^d). \quad (4.29) $$

Finally, combining (4.27) and (4.11), we get a upper bound of $\sigma_h^K$

$$ \|\sigma_h^K\|_{0,K} \lesssim h_K |w_h|_{1,K}. \quad (4.30) $$

### 4.3. Definition of the reconstruction operator

For a given $w_h \in W_h$ and each element $K \in \mathcal{T}_h$, the local interpolation $I_h^{\text{div}} w_h \in W_{k-1}^{\text{div}}(K)$ is defined by (4.7)–(4.9), and $\sigma_h^K$ denotes the solution of the local problem (4.17). Then the local reconstruction operator $R_h^K$ is defined by

$$ R_h^K w_h := I_h^{\text{div}} w_h + \sigma_h^K. $$

The corresponding global operators are denoted by $R_h$, $I_h^{\text{div}}$, and $\sigma_h$, respectively. Based on the preceding construction, $R_h w_h$ can be explicitly computed and satisfies $R_h w_h \in H_0(\text{div}; \Omega)$.

We proceed to verify the three assumptions (3.8)–(3.10) for the operator $R_h$. From (4.10) and (4.28), it is clear to obtain the commutativity

$$ \operatorname{div} R_h w_h = \operatorname{div} I_h^{\text{div}} w_h + \operatorname{div} \sigma_h = \Pi_{k-1}^{0,h} \operatorname{div}_h w_h, \quad \forall w_h \in W_h. $$

Moreover, using (4.11) and (4.30), we obtain the error estimate

$$ \|w_h - R_h w_h\|_{0,K} \le \|w_h - I_h^{\text{div}} w_h\|_{0,K} + \|\sigma_h\|_{0,K} \lesssim h_K |w_h|_{1,K}, \quad \forall K \in \mathcal{T}_h. $$

Finally, the polynomial $L^2$-orthogonality (3.10) follows from (4.29).

## 5. Pressure-robust nonconforming VEM

To present the enhanced version of nonconforming VE [10], we introduce the usual $H^1$-projection $\Pi_k^{\nabla,K}$: $H^1(K) \to \mathbb{P}_k(K)$ by finding a unique polynomial $\Pi_k^{\nabla,K} w$ such that

$$ \left\{ \begin{array}{l} (\nabla \Pi_k^{\nabla,K} w, \nabla q_k)_K = (\nabla w, \nabla q_k)_K, \quad \forall q_k \in \mathbb{P}_k(K), \\ \\ \int_{\partial K} \Pi_k^{\nabla,K} w \,ds = \int_{\partial K} w \,ds. \end{array} \right. $$
The projection $\Pi_k^{\nabla,K}$ can act on the vector-valued functions similarly.

The scalar nonconforming VE space of order $k \ge 1$ on the element $K$ is defined as

$$W_k(K) = \{w \in H^1(K); \Delta w \in \mathbb{P}_k(K), (w - \Pi_k^{\nabla,K} w, q)_K = 0, \forall q \in \mathbb{P}_{k-2}(K), \frac{\partial w}{\partial \mathbf{n}_F} \in \mathbb{P}_{k-1}(F), \forall F \in \mathcal{F}(K)\},$$

where $\mathbb{P}_{k-2}(K)$ denotes the subspace of $\mathbb{P}_k(K)$ consisting of polynomials that are $L^2$-orthogonal to all polynomials in $\mathbb{P}_{k-2}(K)$. Then, we can define the vector version of nonconforming VE space as

$$W_k(K; \mathbb{R}^d) = W_k(K) \otimes \mathbb{R}^d$$

with the following DoFs

$$\bullet \text{ the moments } \frac{1}{|F|}(\mathbf{w}, \mathbf{q}_{k-1})_F, \quad \forall \mathbf{q}_{k-1} \in \mathbb{P}_{k-1}(F; \mathbb{R}^d), \forall F \in \mathcal{F}(K), \qquad (5.1)$$

$$\bullet \text{ the moments } \frac{1}{|K|}(\mathbf{w}, \mathbf{q}_{k-2})_K, \quad \forall \mathbf{q}_{k-2} \in \mathbb{P}_{k-2}(K; \mathbb{R}^d). \qquad (5.2)$$

Note that the $H^1$-projection $\Pi_k^{\nabla,K}$ and the $L^2$-projection $\Pi_k^{0,K}$ are computable using DoFs (5.1)–(5.2). The global VE space $\mathbf{W}_h$ is naturally defined by

$$\mathbf{W}_h = \{\mathbf{w} \in L^2(\Omega; \mathbb{R}^d); \mathbf{w}|_K \in W_k(K; \mathbb{R}^d), \forall K \in \mathcal{T}_h, ([\mathbf{w}], \mathbf{q})_F = 0, \forall \mathbf{q} \in \mathbb{P}_{k-1}(F; \mathbb{R}^d), \forall F \in \mathcal{F}_h\}.$$

In alignment with the conventional approach in the VEM framework, we give the following local bilinear forms by

$$a_h^K(\mathbf{v}_h, \mathbf{w}_h) = (\nabla \Pi_k^{\nabla,K} \mathbf{v}_h, \nabla \Pi_k^{\nabla,K} \mathbf{w}_h)_K + S^K(\mathbf{v}_h - \Pi_k^{\nabla,K} \mathbf{v}_h, \mathbf{w}_h - \Pi_k^{\nabla,K} \mathbf{w}_h).$$

The stabilization term $S^K(\cdot, \cdot)$ is a symmetric and positive definite bilinear form defined by

$$S^K(\mathbf{v}_h, \mathbf{w}_h) = h_K^{d-2} \sum_{i=1}^{N^K} \chi_i(\mathbf{v}_h) \chi_i(\mathbf{w}_h),$$

where $\chi_i$ denotes the $i$-th DoF of nonconforming VE. The bilinear form $S^K(\cdot, \cdot)$ satisfies the following norm equivalence

$$c_1 |\mathbf{w}_h|_{1,K}^2 \le S^K(\mathbf{w}_h, \mathbf{w}_h) \le c_2 |\mathbf{w}_h|_{1,K}^2, \quad \forall \mathbf{w}_h \in \text{Ker}(\Pi_k^{\nabla,K}) \cap W_k(K; \mathbb{R}^d).$$

Then we obtain the $k$-consistency (3.2) and stability (3.3) of the local discrete bilinear form $a_h^K(\cdot, \cdot)$. In the usual way, the computable global bilinear forms $a_h(\cdot, \cdot)$ and $b_h(\cdot, \cdot)$ are given by

$$a_h(\mathbf{v}_h, \mathbf{w}_h) = \sum_{K \in \mathcal{T}_h} a_h^K(\mathbf{v}_h, \mathbf{w}_h), \quad \forall \mathbf{v}_h, \mathbf{w}_h \in \mathbf{W}_h,$$

$$b_h(\mathbf{w}_h, q_h) = \sum_{K \in \mathcal{T}_h} b_h^K(\mathbf{w}_h, q_h), \quad \forall \mathbf{w}_h \in \mathbf{W}_h, \forall q_h \in Q_h,$$

where

$$b_h^K(\mathbf{w}_h, q_h) = -(\text{div}_h \mathbf{w}_h, q_h)_K = -(\Pi_{k-1}^{0,K} \text{div}_h \mathbf{w}_h, q_h)_K.$$

By the way, for any $\mathbf{w}_h \in W_k(K; \mathbb{R}^d)$ and $q_{k-1} \in \mathbb{P}_{k-1}(K)$, the DoFs (5.1)–(5.2) can be used to compute the integrals $(\text{div} \mathbf{w}_h, q_{k-1})_K$, since

$$(\text{div} \mathbf{w}_h, q_{k-1})_K = (\mathbf{w}_h \cdot \mathbf{n}_K, q_{k-1})_{\partial K} - (\mathbf{w}_h, \nabla q_{k-1})_K.$$

Therefore, the interpolation $\Pi_K^{\text{div}} \mathbf{w}_h$ and the right-hand side of the local problem (4.17) can be computed by the DoFs (5.1)–(5.2). Consequently, the reconstruction operator $R_h$ is computable by these VE DoFs (5.1)–(5.2) and satisfies properties (3.8)–(3.10).

With the help of $R_h$, we design the pressure-robust nonconforming VEM by the scheme (3.11). It is well-known that the interpolation result stated in Lemma 3 holds for the space $\mathbf{W}_h$. As discussed in [10, 27, 29], one can prove the discrete inf-sup condition (3.13) and derive the following estimate for the nonconforming error term $N_h$.
**Lemma 10.** Let $u \in H_0^1(\Omega; \mathbb{R}^d) \cap H^{k+1}(\Omega; \mathbb{R}^d)$ be the solution of (1.1) and $w_h \in W_h$. Then, it holds

$$N_h \lesssim h^k |u|_{k+1}. \qquad (5.3)$$

Substituting the nonconforming error estimate (5.3) into the abstract error estimate in Theorem 2, we immediately obtain the following result.

**Theorem 3.** Let $(u, p) \in (H_0^1(\Omega; \mathbb{R}^d) \cap H^{k+1}(\Omega; \mathbb{R}^d)) \times (L_0^2(\Omega) \cap H^k(\Omega))$ and $(u_h, p_h) \in W_h \times Q_h$ be the solutions of problem (1.1) and problem (3.11), respectively. Then we have the pressure-robust error estimates

$$|u - u_h|_{1,h} \lesssim h^k |u|_{k+1}, \qquad (5.4)$$

$$\|p - p_h\|_0 \lesssim h^k (|p|_k + \nu |u|_{k+1}). \qquad (5.5)$$

**Remark 7.** Similarly, for the conforming VEM [6, 7] applied to the Stokes equation, the reconstruction operator $R_h$ developed in Section 4 remains valid for designing the pressure-robust conforming VEM in two and three dimensions. In particular, for the conforming VE function $w_h$, the computability of the reconstruction $R_h w_h$ depends on the computability of the $L^2$-projection $\Pi_k^{0,F} w_h$ on $(d-1)$-faces $F$ and the $L^2$-projection $\Pi_k^{0,K} w_h$ in polytope $K \in \mathcal{T}_h$, which are ensured by using the enhancement technique [1, 22]. Besides, the properties (3.8)-(3.10) of $R_h$ can be obtained only by the definition of $R_h$, which is not related to the DoFs of the VE function. According to Theorem 2, the pressure-robust error estimates with $N_h = 0$ can be directly obtained for the resulting conforming VEM based on the reconstruction.

**Remark 8.** From the discussion in Remark 7, we see that the reconstruction technique developed in Section 4 can be directly used to design the pressure-robust scheme for many other VEMs with discontinuous pressure approximation presented in [2, 26, 27].

## 6. Numerical results

In this section, we provide some numerical tests to verify our theoretical results. The Stokes equation is solved by the pressure-robust nonconforming VEM (3.11) with $k = 1, 2, 3$ for different values of $\nu$ on the unit square $\Omega = (0, 1) \times (0, 1)$, while the right-hand side $f$ is chosen to yield the exact solution

$$u(x_1, x_2) = \begin{pmatrix} \sin(2\pi x_2) - \cos(2\pi x_1) \sin(2\pi x_2) \\ -\sin(2\pi x_1) + \sin(2\pi x_1) \cos(2\pi x_2) \end{pmatrix}, \quad p(x_1, x_2) = x_1 x_2^2 - \frac{1}{6}.$$

Since the VE function is not computable, we compute the absolute errors for the velocity in the discrete energy norm

$$\|u - u_h\|_{1,h} := \sqrt{\sum_{K \in \mathcal{T}_h} |u - \Pi_k^{\nabla,K} u_h|_{1,K}^2}$$

and for the pressure in the usual $L^2$-norm. The unit square $\Omega$ is partitioned into the Voronoi mesh $\mathcal{T}_h^1$ and the nonconvex polygonal mesh $\mathcal{T}_h^2$ as shown in Figure 4.

In order to demonstrate that the reconstruction operator $R_h$ enables the classical nonconforming VEM to achieve pressure robustness, we consider a range of dynamic viscosity values $\nu = 1, 10^{-3}, 10^{-6}, 10^{-9}$. Physically, smaller dynamic viscosity value $\nu$ leads to larger pressure gradient $\nabla p$, so this choice of $\nu$ effectively represents the case of large pressure gradient. The reconstruction operator $R_h$ is defined on a simplicial submesh, which is constructed by connecting the barycenter of each polygonal element to its vertices (see Figure 1c). Two discretizations of the right-hand side are tested: the pressure-robust CR-type nonconforming VEM (3.11) (with reconstruction) and the classical CR-type nonconforming VEM (3.4) (without reconstruction).

Figures 5–16 report the velocity errors in the discrete energy norm $\|\cdot\|_{1,h}$ and the pressure errors in the $L^2$-norm on two types of meshes $\mathcal{T}_h^1$ and $\mathcal{T}_h^2$ for both methods, with polynomial degrees $k = 1, 2, 3$ and various values of the dynamic viscosity $\nu$. For ease of verification, we plot the absolute error curves against the mesh size $h$ in log-log scale, with numbers near the triangles indicating the theoretical convergence rates. From the numerical results, we observe that both the velocity and pressure approximations achieve the expected $k$-order convergence for each method with $k = 1, 2, 3$, which is consistent with the theoretical predictions. As expected, the accuracy of the velocity approximation from the classical CR-type nonconforming VEM
<figure><img src="image_6.png" /><figcaption>Figure 4: The Voronoi (left) and nonconvex (right) meshes.</figcaption></figure>

deteriorates significantly as $\nu$ decreases, while the new method (3.11) with the reconstruction operator remains pressure-robust.

We next verify the necessity of solving the local problem (4.17) in the reconstruction. Theoretical analysis from Section 4 indicates that for low-order methods ($k = 1, 2$), the reconstruction without solving local problems can lead to the optimal convergence rate. However, for higher-order methods ($k \ge 3$), the solution of local problems ensure the $L^2$-orthogonality of the reconstruction operator and become essential to achieve the optimal convergence. To illustrate this, Tables 1 and 2 report numerical results for $k = 3$ using the nonconforming VEM (3.11), where the reconstruction is done only by using the $H(\text{div})$-conforming macro FE interpolation $I_h^{\text{div}}$, i.e., without solving any local problems. It is evident from the numerical results that this simplified version of the method fails to deliver the optimal convergence, specially for the velocity approximation. Therefore, for the higher-order methods ($k \ge 3$), the solution of local problem is necessary in the reconstruction in order to obtain the pressure-robust scheme with the optimal convergence.

<figure><img src="image_7.png" /><figcaption>Figure 5: The velocity error $\|u - u_h\|_{1,h}$ for $k = 1$ on mesh $\mathcal{T}_h^1$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: The velocity error $||u - u_h||_{1,h}$ for $k = 2$ on mesh $T_h^1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: The velocity error $||u - u_h||_{1,h}$ for $k = 3$ on mesh $T_h^1$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: The pressure error $||p - p_h||_0$ for $k = 1$ on mesh $T_h^1$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: The pressure error $||p - p_h||_0$ for $k = 2$ on mesh $\mathcal{T}_h^1$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: The pressure error $||p - p_h||_0$ for $k = 3$ on mesh $\mathcal{T}_h^1$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: The velocity error $||u - u_h||_{1,h}$ for $k = 1$ on mesh $\mathcal{T}_h^2$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12: The velocity error $||u - u_h||_{1,h}$ for $k = 2$ on mesh $\mathcal{T}_h^2$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: The velocity error $||u - u_h||_{1,h}$ for $k = 3$ on mesh $\mathcal{T}_h^2$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: The pressure error $||p - p_h||_0$ for $k = 1$ on mesh $\mathcal{T}_h^2$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15: The pressure error $||p - p_h||_0$ for $k = 2$ on mesh $\mathcal{T}_h^2$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16: The pressure error $||p - p_h||_0$ for $k = 3$ on mesh $\mathcal{T}_h^2$.</figcaption></figure>

<table><caption>Table 1: Numerical results for $k = 3$ on mesh $\mathcal{T}_h^1$ without solving local problems.</caption><caption>the velocity error ||u - u_h||<sub>1,h</sub></caption><thead><tr><th>h</th><th>0.355305</th><th>0.197791</th><th>0.102512</th><th>0.049521</th><th>0.023756</th><th>Rate</th></tr></thead><tbody><tr><th>1</th><td>5.9224e-01</td><td>8.3492e-02</td><td>1.3707e-02</td><td>2.7713e-03</td><td>6.0939e-04</td><td>2.06</td></tr><tr><th>10<sup>-3</sup></th><td>5.9224e-01</td><td>8.3492e-02</td><td>1.3707e-02</td><td>2.7713e-03</td><td>6.0939e-04</td><td>2.06</td></tr><tr><th>10<sup>-6</sup></th><td>5.9224e-01</td><td>8.3492e-02</td><td>1.3707e-02</td><td>2.7713e-03</td><td>6.0939e-04</td><td>2.06</td></tr><tr><th>10<sup>-9</sup></th><td>5.9224e-01</td><td>8.3492e-02</td><td>1.3707e-02</td><td>2.7713e-03</td><td>6.0939e-04</td><td>2.06</td></tr></tbody></table>

<table><caption>the pressure error ||p - p_h||<sub>0</sub></caption><thead><tr><th>h</th><th>0.355305</th><th>0.197791</th><th>0.102512</th><th>0.049521</th><th>0.023756</th><th>Rate</th></tr></thead><tbody><tr><th>1</th><td>3.7672e-01</td><td>5.7717e-02</td><td>9.6027e-03</td><td>1.8684e-03</td><td>3.9469e-04</td><td>2.12</td></tr><tr><th>10<sup>-3</sup></th><td>4.9629e-04</td><td>6.8882e-05</td><td>1.0637e-05</td><td>1.9469e-06</td><td>4.0043e-07</td><td>2.15</td></tr><tr><th>10<sup>-6</sup></th><td>3.2308e-04</td><td>3.7597e-05</td><td>4.5744e-06</td><td>5.4738e-07</td><td>6.7542e-08</td><td>2.85</td></tr><tr><th>10<sup>-9</sup></th><td>3.2308e-04</td><td>3.7597e-05</td><td>4.5744e-06</td><td>5.4737e-07</td><td>6.7541e-08</td><td>2.85</td></tr></tbody></table>
<table><caption>Table 2: Numerical results for $k = 3$ on mesh $\mathcal{T}_h^2$ without solving local problems.</caption><thead><tr><th colspan="7">the velocity error ||u - u<sub>h</sub>||<sub>1,h</sub></th></tr><tr><th>h</th><th>0.364434</th><th>0.182217</th><th>0.091109</th><th>0.045554</th><th>0.022777</th><th>Rate</th></tr></thead><tbody><tr><th>1</th><td>6.3263e-01</td><td>1.1898e-01</td><td>2.2227e-02</td><td>4.1897e-03</td><td>8.5817e-04</td><td>2.23</td></tr><tr><th>10<sup>-3</sup></th><td>6.3263e-01</td><td>1.1898e-01</td><td>2.2227e-02</td><td>4.1897e-03</td><td>8.5817e-04</td><td>2.23</td></tr><tr><th>10<sup>-6</sup></th><td>6.3263e-01</td><td>1.1898e-01</td><td>2.2227e-02</td><td>4.1897e-03</td><td>8.5817e-04</td><td>2.23</td></tr><tr><th>10<sup>-9</sup></th><td>6.3263e-01</td><td>1.1898e-01</td><td>2.2227e-02</td><td>4.1897e-03</td><td>8.5817e-04</td><td>2.23</td></tr></tbody></table><table><thead><tr><th colspan="7">the pressure error ||p - p<sub>h</sub>||<sub>0</sub></th></tr><tr><th>h</th><th>0.364434</th><th>0.182217</th><th>0.091109</th><th>0.045554</th><th>0.022777</th><th>Rate</th></tr></thead><tbody><tr><th>1</th><td>4.0843e-01</td><td>8.2562e-02</td><td>1.5460e-02</td><td>2.6288e-03</td><td>4.4783e-04</td><td>2.55</td></tr><tr><th>10<sup>-3</sup></th><td>5.3521e-04</td><td>9.2903e-05</td><td>1.6337e-05</td><td>2.7097e-06</td><td>4.5527e-07</td><td>2.57</td></tr><tr><th>10<sup>-6</sup></th><td>3.4588e-04</td><td>4.2597e-05</td><td>5.2808e-06</td><td>6.5723e-07</td><td>8.1971e-08</td><td>3.00</td></tr><tr><th>10<sup>-9</sup></th><td>3.4588e-04</td><td>4.2597e-05</td><td>5.2807e-06</td><td>6.5722e-07</td><td>8.1970e-08</td><td>3.00</td></tr></tbody></table>

## Acknowledgments

This work was supported by National Natural Science Foundation of China (No. 12371411). The author Jikun Zhao is funded by the Training Program for Young Backbone Teachers in Higher Education Institutions of Henan Province. The author Shipeng Mao was supported by National Key Research and Development Program of China (2024YFA1012502) and the National Natural Science Foundation of China (12271514).

## References

[1] B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo. Equivalent projectors for virtual element methods. *Comput. Math. Appl.*, 66(3):376-391, 2013.

[2] P. F. Antonietti, L. Beirão da Veiga, D. Mora, and M. Verani. A stream virtual element formulation of the Stokes problem on polygonal meshes. *SIAM J. Numer. Anal.*, 52(1):386-404, 2014.

[3] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus, homological techniques, and applications. *Acta Numer.*, 15:1-155, 2006.

[4] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. Basic principles of virtual element methods. *Math. Models Methods Appl. Sci.*, 23(1):199-214, 2013.

[5] L. Beirão da Veiga, F. Dassi, D. Di Pietro, and J. Droniou. Arbitrary-order pressure-robust DDR and VEM methods for the Stokes problem on polyhedral meshes. *Comput. Methods Appl. Mech. Engrg.*, 397:Paper No. 115061, 31, 2022.

[6] L. Beirão da Veiga, F. Dassi, and G. Vacca. The Stokes complex for virtual elements in three dimensions. *Math. Models Methods Appl. Sci.*, 30(3):477-512, 2020.

[7] L. Beirão da Veiga, C. Lovadina, and G. Vacca. Divergence free virtual elements for the Stokes problem on polygonal meshes. *ESAIM Math. Model. Numer. Anal.*, 51(2):509-535, 2017.

[8] S. C. Brenner, Q. Guan, and L.-Y. Sung. Some estimates for virtual element methods. *Comput. Methods Appl. Math.*, 17(4):553-574, 2017.

[9] A. Cangiani, Z. Dong, and E. H. Georgoulis. hp-version discontinuous Galerkin methods on essentially arbitrarily-shaped elements. *Math. Comp.*, 91(333):1-35, 2021.
[10] A. Cangiani, V. Gyrya, and G. Manzini. The nonconforming virtual element method for the Stokes equations. *SIAM J. Numer. Anal.*, 54(6):3411–3435, 2016.

[11] C. Chen, X. Huang, and H. Wei. Virtual element methods without extrinsic stabilization. *SIAM J. Numer. Anal.*, 62(1):567–591, 2024.

[12] M. Crouzeix and P.-A. Raviart. Conforming and nonconforming finite element methods for solving the stationary Stokes equations. I. *Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge*, 7:33–75, 1973.

[13] D. A. Di Pietro, J. Droniou, and J. J. Qian. A pressure-robust discrete de Rham scheme for the Navier-Stokes equations. *Comput. Methods Appl. Mech. Engrg.*, 421:Paper No. 116765, 21, 2024.

[14] D. Frerichs and C. Merdon. Divergence-preserving reconstructions on polygons and a really pressure-robust virtual element method for the Stokes problem. *IMA J. Numer. Anal.*, 42(1):597–619, 2022.

[15] V. Girault and P.-A. Raviart. *Finite element methods for Navier-Stokes equations: Theory and algorithms, volume 5* of *Springer Series in Computational Mathematics*. Springer-Verlag, Berlin, 1986.

[16] V. John, A. Linke, C. Merdon, M. Neilan, and L. G. Rebholz. On the divergence constraint in mixed finite element methods for incompressible flows. *SIAM Rev.*, 59(3):492–544, 2017.

[17] P. L. Lederer, A. Linke, C. Merdon, and J. Schöberl. Divergence-free reconstruction operators for pressure-robust Stokes discretizations with continuous pressure finite elements. *SIAM J. Numer. Anal.*, 55(3):1291–1314, 2017.

[18] H. Liang and H. Rui. A parameter robust reconstruction nonconforming virtual element method for the incompressible poroelasticity model. *Appl. Numer. Math.*, 202:127–142, 2024.

[19] 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.

[20] A. Linke and C. Merdon. Guaranteed energy error estimators for a modified robust Crouzeix-Raviart Stokes element. *J. Sci. Comput.*, 64(2):541–558, 2015.

[21] D. Mora, C. Reales, and A. Silgado. A $C^1$-virtual element method of high order for the Brinkman equations in stream function formulation with pressure recovery. *IMA J. Numer. Anal.*, 42(4):3632–3674, 2022.

[22] G. Vacca. An $H^1$-conforming virtual element for Darcy and Brinkman equations. *Math. Models Methods Appl. Sci.*, 28(1):159–194, 2018.

[23] G. Wang, L. Mu, Y. Wang, and Y. He. A pressure-robust virtual element method for the Stokes problem. *Comput. Methods Appl. Mech. Engrg.*, 382:Paper No. 113879, 19, 2021.

[24] G. Wang, Y. Wang, and Y. He. Two robust virtual element methods for the Brinkman equations. *Calcolo*, 58(4):Paper No. 49, 27, 2021.

[25] Y. Wang, G. Wang, and Y. Shen. A pressure-robust virtual element method for the Navier-Stokes problem on polygonal mesh. *Comput. Math. Appl.*, 131:124–137, 2023.

[26] H. Wei, X. Huang, and A. Li. Piecewise divergence-free nonconforming virtual elements for Stokes problem in any dimensions. *SIAM J. Numer. Anal.*, 59(3):1835–1856, 2021.

[27] J. Zhao, B. Zhang, S. Mao, and S. Chen. The divergence-free nonconforming virtual element for the Stokes problem. *SIAM J. Numer. Anal.*, 57(6):2730–2759, 2019.

[28] J. Zhao, B. Zhang, S. Mao, and S. Chen. The nonconforming virtual element method for the Darcy-Stokes problem. *Comput. Methods Appl. Mech. Engrg.*, 370:113251, 23, 2020.

[29] J. Zhao, W. Zhu, B. Zhang, and Y. Yang. The stabilized nonconforming virtual element method for the Darcy-Stokes problem. *Commun. Nonlinear Sci. Numer. Simul.*, 138:Paper No. 108252, 14, 2024.