# Stability and interpolation properties for Stokes-like virtual element
spaces

J. Meng , L. Beirão da Veiga , L. Mascotto22footnotemark: 2 33footnotemark: 3
School of Mathematics and Statistics, Xi’an Jiaotong University, 710049,
Shaanxi, P.R. China<EMAIL_ADDRESS>di Matematica e
Applicazioni, Università di Milano Bicocca, 20125 Milan, Italy,
<EMAIL_ADDRESS><EMAIL_ADDRESS>27100, Pavia,
ItalyFakultät für Mathematik, Universität Wien, 1090 Vienna, Austria,
<EMAIL_ADDRESS>
###### Abstract

We prove stability bounds for Stokes-like virtual element spaces in two and
three dimensions. Such bounds are also instrumental in deriving optimal
interpolation estimates. Furthermore, we develop some numerical tests in order
to investigate the behaviour of the stability constants also from the
practical side.

AMS subject classification: 65N12; 65N30; 65N50.

Keywords: virtual element method; optimal convergence; stability; Stokes
problem.

##  1 Introduction

In recent years, due to their flexibility in handling complex data features
and adaptive mesh refinements, Galerkin methods based on polytopal meshes
received an increasing attention. The virtual element method (VEM) [5] is one
amid the most successful of such polytopal methods.

Amongst the various problems that have been tackled with the VEM, fluid static
and dynamic problems have a prominent role. The first paper coping with a
lowest order VEM for the Stokes problem is [4]. Later, its general order
conforming [11] and nonconforming versions [22, 35] have been discussed. Based
on that, conforming [12] and nonconforming VEMs for the Navier-Stokes problem
[34] were also introduced.

In addition to such works that represent the backbone of the VEM for fluid-
type problems, other similar topics have been studied as well, an incomplete
and short list being: mixed VEMs for the pseudo-stress-velocity formulation of
the Stokes problem [19]; mixed VEMs for quasi-Newtonian flows [21]; mixed VEMs
for the Navier-Stokes problem [32]; other variants of the VEM for the Darcy
problem [42, 40, 20, 37]; the analysis of the Stokes complex in the VEM
framework [13, 9]; a stabilized VEM for the unsteady incompressible Navier-
Stokes problem [33]; implementation details [29]; a pressure robust variant of
the VEM for the Stokes problem [31]; the magneto-hydrodynamic problem [8]; the
$hp$-version of the standard VEM for the Stokes problem [25]; stationary
quasi-geostrophic equations of the ocean [38]; the unsteady Navier-Stokes
problem [1].

Needless to write, other polytopal methods have been used to approximate the
above problems. For instance, we recall the local discontinuous Galerkin
method [26]; hybrid discontinuous Galerkin schemes for the Stokes flow [28,
27]; hybrid discontinuous Galerkin schemes for the Navier-Stokes problem [39];
hybrid high-order methods [2, 30, 16, 18].

The analysis of the VEM is based on showing optimal a priori error estimates,
which are proved by means of certain _stability_ and _polynomial consistency_
properties. Optimal convergence is then derived based on using best polynomial
and interpolation estimates.

All in all, the main difference with respect to the finite element setting
resides in further employing

  * •

interpolation estimates in virtual element spaces;

  * •

stability properties on a discrete bilinear form.

Several interpolation estimates for Stokes-type virtual element spaces are
available in the literature; see, e.g., [11, Proposition $4.2$] and [12,
Theorem $4.1$]. Such interpolation estimates are rather technical to prove,
are strongly hinge upon the definition of the local virtual element spaces,
and have been proved in two dimensions only.

On the other hand, to the best of our knowledge, stability properties for
Stokes-like virtual element spaces have never been explicitly proved. We point
out that several works coping with explicit stability estimates for standard
Poisson-like virtual element spaces are available; see, e.g., [17, 10, 24, 23,
7].

This paper aims at closing this theoretical gap. Notably, we contribute to the
current state of the art along the three following avenues:

  1. 1.

We prove explicit stability properties for Stokes-like virtual element spaces
in two and three dimensions. To this aim, we employ two different
stabilizations, namely one given by the inner product of the degrees of
freedom and one in integral form, which is independent of the chosen degrees
of freedom;

  2. 2.

Based on such stability estimates, we prove novel interpolation estimates,
which deliver the same convergence as those already available in the
literature but are much easier to prove.

  3. 3.

The above achievements are proven for regular polytopal meshes. Thus, we also
exhibit numerical results investigating the stability constants for sequences
of meshes with degenerating geometry as well as with respect to the degree of
accuracy of the method.

In the remainder of the introduction, we pinpoint some basic notation of the
paper, discuss the model problem we aim to approximate, introduce sequences of
regular polytopal meshes, and detail the structure of the paper.

#### Notation.

Throughout, we employ standard notation for Lebesgue and Sobolev spaces.
Notably, given a domain $D\subset\mathbb{R}^{d}$, $d=1,2,3$, $L^{2}(D)$
denotes the space of measurable and integrable squared functions and
$L^{2}_{0}(D)$ its subspace consisting of functions with zero average over
$D$. Given $r\in\mathbb{N}$, $H^{r}(D)$ denotes the Sobolev space of order
$r$, i.e., the subspace of $L^{2}(D)$ consisting of functions with integrable
squared weak derivatives $D^{r}\cdot$ up to order $r$ (we conventionally set
$H^{0}(D)=L^{2}(D)$). We endow the above spaces with the standard bilinear
forms $(\cdot,\cdot)_{r,D}:=(D^{r}\cdot,D^{r}\cdot)_{0,D}$ and (semi)norms
$\|\cdot\|_{r,D}^{2}:=\sum_{\ell=0}^{r}(\cdot,\cdot)_{\ell,D}$ and
$|\cdot|_{r,D}^{2}:=(\cdot,\cdot)_{r,D}$. We also introduce $H^{1}_{0}(D)$ as
the space of $H^{1}$ functions with zero trace over the boundary $\partial D$
of $D$. Noninteger order Sobolev spaces can be constructed by interpolation.

Given $\ell\in\mathbb{N}$, $\mathbb{P}_{\ell}(D)$ denotes the space of
polynomials of degree at most $\ell$ over $D$ and
$\mathbb{P}_{\ell}(D)\setminus\mathbb{R}:=\mathbb{P}_{\ell}(D)\cap
L^{2}_{0}(D)$.

Finally, given two positive quantities $a$ and $b$, we use the short-hand-
notation $a\lesssim b$ instead of _there exists a positive constant $c$
independent of the mesh such that_ $a\leq c\ b$. We further write $a\approx b$
if we have $a\lesssim b$ and $b\lesssim a$.

#### The model problem.

Let $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, be an open domain and
$\mathbf{f}\in[L^{2}(\Omega)]^{d}$.

As a model problem, we consider the Stokes problem

$\begin{cases}\text{find }(\mathbf{u},s)\text{ such that}\\\
-\boldsymbol{\Delta}\mathbf{u}-\nabla s=\mathbf{f}&\text{in }\Omega\\\
-\operatorname{div}\mathbf{u}=0&\text{in }\Omega\\\
\mathbf{u}=\mathbf{0}&\text{on }\partial\Omega,\end{cases}$

which in weak formulation reads as follows:

$\begin{cases}\text{find }(\mathbf{u},s)\in[H^{1}_{0}(\Omega)]^{d}\times
L^{2}_{0}(\Omega):=\mathbf{V}\times Q\text{ such that}\\\
(\boldsymbol{\nabla}\mathbf{u},\boldsymbol{\nabla}\mathbf{v})_{0,\Omega}+(\operatorname{div}\mathbf{v},s)_{0,\Omega}=(\mathbf{f},\mathbf{v})_{0,\Omega}&\forall\mathbf{v}\in\mathbf{V}\\\
(\operatorname{div}\mathbf{u},t)_{0,\Omega}=0&\forall t\in Q.\end{cases}$

The well posedness of this problem is standard [15].

#### Regular polytopal meshes

Throughout, we are given sequences $\\{\mathcal{T}_{n}\\}$ of polytopal meshes
over the domain $\Omega$. In $d=2$, $\mathcal{T}_{n}$ consists of conforming
polygons; in $d=3$, $\mathcal{T}_{n}$ consists of conforming polyhedra. We
denote a generic element of $\mathcal{T}_{n}$ by $K$; $\partial K$ denotes the
boundary of $K$ with outward unit vector $\mathbf{n}_{K}$. For any geometric
object $D\subset\mathbb{R}^{d}$, $d=1,2,3$, we denote its barycenter, measure,
and diameter by $\mathbf{x}_{D}$, $|D|$, and $h_{D}$, respectively.

Given an element $K$ in three dimensions, $\partial K$ is the union of its
faces $F$. Given an element in two dimensions or a face $F$ in three
dimensions, $\partial F$ is the union of its edges $e$.

We demand standard regularity assumptions on $\mathcal{T}_{n}$: there exists
$\rho>0$ such that

  * •

for $d=2$,

    * –

every polygon $F$ is star-shaped with respect to a disk of diameter greater
than or equal to $\rho h_{F}$;

    * –

every edge $e$ satisfies $h_{e}\geq\rho h_{F}$;

  * •

for $d=3$,

    * –

every polyhedron $K$ is star-shaped with respect to a disk of diameter greater
than or equal to $\rho h_{F}$;

    * –

every face $F$ of $K$ is star-shaped with respect to a disk of diameter
greater than or equal to $\rho h_{F}$;

    * –

for every face $F$ of $K$ and edge $e$ of $F$, we have $h_{e}\geq\rho
h_{F}\geq\rho^{2}h_{e}$.

Given $\mathbf{x}=(x_{1},y_{1})$, we define
$\mathbf{x}^{\\!\perp}:=(x_{2},-x_{1})$. We denote the set of edges of a
polytope $K$ by $\mathcal{E}^{K}$ and the set of faces of a polyhedron $K$ by
$\mathcal{E}^{F}$.

#### Outline of the paper.

In Sections 2 and 3, after recalling the definition of Stokes-type virtual
element spaces [11, 9], we prove stability and interpolation properties in two
and three dimensions, respectively. We perform the analysis for two explicit
stabilizations. In Section 4, we numerically check the stability properties
(in 2D only) on sequences of meshes with degenerating geometry and degree of
accuracy. We draw some conclusions in Section 5.

##  2 The two dimensional case

This section is devoted at proving stability and interpolation properties in
two dimensions. In Section 2.1, we recall the definition of the Stokes-like
virtual element space [11]. Stability properties are derived in Sections 2.2
and 2.3 for a projection based and a degrees of freedom based stabilizations,
respectively. In Section 2.4, we provide a novel and shorter proof of
interpolation estimates based on the previously proven stability properties.

In what follows, we fix $p\in\mathbb{N}$, which will denote the degree of
accuracy of the space.

###  2.1 Virtual element spaces in two dimensions

Given a polygon $K\in\mathcal{T}_{n}$, introduce the space

$\mathbf{V}_{n}(K):=\\{\mathbf{v}_{n}\in[H^{1}(K)]^{2}\mid\mathbf{v}_{n}\text{
satisfies }\eqref{local-auxiliary-problem}\\},$

where, for some $s\in L^{2}_{0}(K)$,

$\begin{cases}-\boldsymbol{\Delta}\mathbf{v}_{n}-\nabla
s=\mathbf{x}^{\\!\perp}q_{p-3}&q_{p-3}\in\mathbb{P}_{p-3}(K)\\\\[4.0pt]
\operatorname{div}\mathbf{v}_{n}=q_{p-1}&q_{p-1}\in\mathbb{P}_{p-1}(K)\\\\[4.0pt]
\mathbf{v}_{n}{}_{|\partial K}\in[\mathcal{C}^{0}(\partial
K)]^{2},\quad\mathbf{v}_{n}{}_{|e}\in[\mathbb{P}_{p}(e)]^{2}&\forall
e\in\mathcal{E}^{K},\end{cases}$ (1)

with all equations to be intended in a weak sense.

We endow the space $\mathbf{V}_{n}(K)$ with the following set of unisolvent
degrees of freedom (DoFs) [11]: given $\mathbf{v}_{n}\in\mathbf{V}_{n}(K)$

  * •

the vector values Dv1($\mathbf{v}_{n}$) at the vertices of $K$;

  * •

the vector values Dv2($\mathbf{v}_{n}$) at the $p-1$ internal Gauß-Lobatto
nodes on each edge $e$ of $\mathcal{E}^{K}$;

  * •

for $p\geq 3$, given $\\{\mathbf{q}_{\boldsymbol{\alpha}}\\}$ a basis of
$\mathbf{x}^{\\!\perp}\mathbb{P}_{p-3}(K)$, the “orthogonal” moments

$\textbf{Dv}_{3}(\mathbf{v}_{n}):=\frac{1}{|K|}\int_{K}\mathbf{v}_{n}\cdot\mathbf{q}_{\boldsymbol{\alpha}};$
(2)

  * •

given $\\{m_{\boldsymbol{\alpha}}\\}$ a basis of
$\mathbb{P}_{p-1}(K)\setminus\mathbb{R}$, the “divergence” moments

$\textbf{Dv}_{4}(\mathbf{v}_{n}):=\frac{h_{K}}{|K|}\int_{K}\operatorname{div}\mathbf{v}_{n}\
m_{\boldsymbol{\alpha}}.$ (3)

As usual, we require that the bases $\\{\mathbf{q}_{\boldsymbol{\alpha}}\\}$
and $\\{m_{\boldsymbol{\alpha}}\\}$ are invariant with respect to translations
and dilations; see [3, 29].

It is known [11] that, given $\mathbf{v}_{n}\in\mathbf{V}_{n}(K)$ with known
DoFs, $\operatorname{div}\ \mathbf{v}_{n}$ is explicitly computable. Further,
we can compute the two orthogonal projectors
$\boldsymbol{\Pi}^{0}_{p-2}:[L^{2}(K)]^{2}\to[\mathbb{P}_{p-2}(K)]^{2}$ and
$\boldsymbol{\Pi}^{0}_{\perp,p-3}:[L^{2}(K)]^{2}\to\mathbf{x}^{\\!\perp}\mathbb{P}_{p-3}(K)$
defined as follows: for all $\mathbf{v}_{n}$ in $\mathbf{V}_{n}(K)$,

$\begin{split}&(\mathbf{q}_{p-2},\mathbf{v}_{n}-\boldsymbol{\Pi}^{0}_{p-2}\mathbf{v}_{n})_{0,K}=0\qquad\quad\;\;\;\forall\mathbf{q}_{p-2}\in[\mathbb{P}_{p-2}(K)]^{2},\\\\[4.0pt]
&(\mathbf{x}^{\\!\perp}q_{p-3},\mathbf{v}_{n}-\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{v}_{n})_{0,K}=0\qquad\forall
q_{p-3}\in\mathbb{P}_{p-3}(K).\end{split}$ (4)

We can also compute the $H^{1}$ projector
$\boldsymbol{\Pi}^{\nabla}_{p}:[L^{2}(K)]^{2}\to[\mathbb{P}_{p}(K)]^{2}$
defined as

$(\boldsymbol{\nabla}\mathbf{q}_{p},\boldsymbol{\nabla}(\mathbf{v}_{n}-\boldsymbol{\Pi}^{\nabla}_{p}\mathbf{v}_{n}))_{0,K}=0,\qquad\int_{\partial
K}\mathbf{v}_{n}=\mathbf{0}\qquad\forall\mathbf{v}_{n}\in\mathbf{V}_{n}(K),\;\forall\mathbf{q}_{p}\in[\mathbb{P}_{p}(K)]^{2}.$
(5)

Functions in the virtual element space $\mathbf{V}_{n}(K)$, as well as their
gradients, are not available in closed form.

For this reason, following the virtual element gospel [5], we discretize the
bilinear form $(\nabla\cdot,\nabla\cdot)_{0,\Omega}$ as follows: for all
$K\in\mathcal{T}_{n}$,

$\begin{split}a_{h}^{K}(\mathbf{u}_{n},\mathbf{v}_{n})&:=a^{K}(\boldsymbol{\Pi}^{\nabla}_{p}\mathbf{u}_{n},\boldsymbol{\Pi}^{\nabla}_{p}\mathbf{v}_{n})+S^{K}((\mathbf{I}-\boldsymbol{\Pi}^{\nabla}_{p})\mathbf{u}_{n},(\mathbf{I}-\boldsymbol{\Pi}^{\nabla}_{p})\mathbf{v}_{n})\\\
&:=(\boldsymbol{\nabla}\boldsymbol{\Pi}^{\nabla}_{p}\mathbf{u}_{n},\boldsymbol{\nabla}\boldsymbol{\Pi}^{\nabla}_{p}\mathbf{v}_{n})_{0,K}+S^{K}((\mathbf{I}-\boldsymbol{\Pi}^{\nabla}_{p})\mathbf{u}_{n},(\mathbf{I}-\boldsymbol{\Pi}^{\nabla}_{p})\mathbf{v}_{n}).\end{split}$
(6)

The bilinear form $S^{K}(\cdot,\cdot)$ is required to be coercive and
continuous on $\mathbf{V}_{n}(K)\cap\ker(\boldsymbol{\Pi}^{\nabla}_{p})$
uniformly in the mesh elements. More precisely, on this space, we require
$S^{K}(\cdot,\cdot)\approx|\cdot|^{2}_{1,K}$. Proving such an equivalence is
our goal in the forthcoming sections, for two explicit stabilizations.

The global counterpart of the space $\mathbf{V}_{n}(K)$ is constructed by a
standard $H^{1}$-conforming DoFs-coupling.

######  Remark 1.

Due to the mesh regularity assumptions, each $K\in\mathcal{T}_{n}$ can be
split into the union of shape regular simplices

$K=\cup_{T\in\widetilde{\mathcal{T}}_{n}^{K}}T.$

This applies both in the two and the three dimensional cases.

###  2.2 Stability estimates for a projection based stabilization

Given $K\in\mathcal{T}_{n}$, we consider the local stabilization

$S^{K}(\mathbf{u}_{n},\mathbf{v}_{n}):=h_{K}^{-2}(\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{u}_{n},\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{v}_{n})_{0,K}+(\operatorname{div}\mathbf{u}_{n},\operatorname{div}\mathbf{v}_{n})_{0,K}+h_{K}^{-1}(\mathbf{u}_{n},\mathbf{v}_{n})_{0,\partial
K},$ (7)

where we recall the polynomial bases in the definition of the DoFs (2) and (3)
are defined so as to be invariant with respect to translation and dilations.

######  Theorem 2.1.

The following stability bounds are valid: there exist
$0<\alpha_{*}<\alpha^{*}$ independent of $h_{K}$ such that

$\displaystyle\alpha_{*}|\mathbf{v}_{n}|_{1,K}^{2}\leq
S^{K}(\mathbf{v}_{n},\mathbf{v}_{n})\qquad\forall\mathbf{v}_{n}\in\mathbf{V}_{n}(K),$
(8) $\displaystyle
S^{K}(\mathbf{v}_{n},\mathbf{v}_{n})\leq\alpha^{*}|\mathbf{v}_{n}|_{1,K}^{2}\qquad\forall\mathbf{v}_{n}\in[H^{1}(K)]^{2}\text{
such that }\int_{\partial K}\mathbf{v}_{n}=\mathbf{0}.$ (9)

In particular, bounds (8) and (9) are valid for functions in
$\mathbf{V}_{n}(K)\cap\ker(\boldsymbol{\Pi}^{\nabla}_{p})$.

###### Proof.

We begin by proving (8) splitting its proof into four smaller steps. Let
$\mathbf{v}_{n}\in\mathbf{V}_{n}(K)$ solve (1) and $s\in L^{2}_{0}(K)$ be the
associated auxiliary pressure.

Preliminary fact $1$. We first observe that an integration by parts yields

$|\boldsymbol{\Delta}\mathbf{v}_{n}|_{-1,K}:=\sup_{\mathbf{0}\neq\boldsymbol{\Phi}\in
H^{1}_{0}(K)}\frac{(\boldsymbol{\Delta}\mathbf{v}_{n},\boldsymbol{\Phi})_{0,K}}{\|\boldsymbol{\Phi}\|_{1,K}}\lesssim|\mathbf{v}_{n}|_{1,K}.$
(10)

Preliminary fact $2$. If $\operatorname{rot}(\mathbf{x}^{\\!\perp}q_{p-3})=0$,
then $\mathbf{x}^{\\!\perp}q_{p-3}=\mathbf{0}$. Thence, a scaling argument and
equivalence of norms for polynomial spaces of finite maximum degree imply

$\|\mathbf{x}^{\\!\perp}q_{p-3}\|_{0,K}\lesssim
h_{K}\|\operatorname{rot}(\mathbf{x}^{\\!\perp}q_{p-3})\|_{0,K}.$ (11)

In light of this, we show an auxiliary bound on the $L^{2}$ norm of the right-
hand side in (1). Using (11) and the first line in (1), we write

$\|\mathbf{x}^{\\!\perp}q_{p-3}\|_{0,K}\lesssim
h_{K}\|\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n})\|_{0,K}.$

Next, let $b^{K}$ denote the piecewise cubic bubble function over the sub
tessellation $\widetilde{\mathcal{T}}_{n}^{K}$ of $K$ introduced in Remark 1
such that $\|b^{K}\|_{\infty,\widetilde{T}}=1$ for all $\widetilde{T}$ in
$\widetilde{\mathcal{T}}_{n}^{K}$. Since
$\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n})$ is a polynomial, the
following polynomial inverse estimate involving bubbles holds true:

$\|\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n})\|_{0,K}^{2}\lesssim\|b^{K}\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n})\|_{0,K}^{2}.$

Integrating by parts twice, observing that $(b^{K})^{2}$ and its normal
derivative annihilate on the boundary of each $\widetilde{T}$ in
$\widetilde{\mathcal{T}}_{n}^{K}$, and using the Cauchy-Schwarz inequality and
a polynomial inverse inequality twice, we arrive at

$\begin{split}\|\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n})\|_{0,K}^{2}&\lesssim\left(\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n}),(b^{K})^{2}\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n})\right)_{0,K}=\left(\boldsymbol{\nabla}\mathbf{v}_{n},\boldsymbol{\nabla}\operatorname{\textbf{curl}}((b^{K})^{2}\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n}))\right)_{0,K}\\\\[4.0pt]
&\lesssim|\mathbf{v}_{n}|_{1,K}\;h_{K}^{-2}\|(b^{K})^{2}\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n})\|_{0,K}\leq|\mathbf{v}_{n}|_{1,K}\;h_{K}^{-2}\|\operatorname{rot}(\boldsymbol{\Delta}\mathbf{v}_{n})\|_{0,K}.\end{split}$

Combining the two above bounds yields

$\|\mathbf{x}^{\\!\perp}q_{p-3}\|_{0,K}\lesssim
h_{K}^{-1}|\mathbf{v}_{n}|_{1,K}.$ (12)

Preliminary fact $3$. We show a technical bound on the auxiliary pressure $s$
in (1). To this aim, we first observe that a scaled Poincaré inequality
entails

$\|\mathbf{v}\|_{-1,K}\lesssim
h_{K}\|\mathbf{v}\|_{0,K}\qquad\forall\mathbf{v}\in[L^{2}(K)]^{2}.$ (13)

Next, resorting to the standard inf-sup condition [15] for the couple
$[H^{1}_{0}(K)]^{2}\times L^{2}_{0}(K)$, integrating by parts, and using the
triangle inequality, the first equation in (1), (13), (10), and (12), we
deduce

$\begin{split}\|s\|_{0,K}&\lesssim\sup_{\mathbf{0}\neq\boldsymbol{\Phi}\in[H^{1}_{0}(K)]^{2}}\frac{(s,\operatorname{div}\boldsymbol{\Phi})_{0,K}}{|\boldsymbol{\Phi}|_{1,K}}=\sup_{\mathbf{0}\neq\boldsymbol{\Phi}\in[H^{1}_{0}(K)]^{2}}\frac{(\nabla
s,\boldsymbol{\Phi})_{0,K}}{|\boldsymbol{\Phi}|_{1,K}}=\|\nabla
s\|_{-1,K}\\\\[4.0pt]
&\leq\|\mathbf{x}^{\\!\perp}q_{p-3}\|_{-1,K}+\|\boldsymbol{\Delta}\mathbf{v}_{n}\|_{-1,K}\lesssim
h_{K}\|\mathbf{x}^{\\!\perp}q_{p-3}\|_{0,K}+|\mathbf{v}_{n}|_{1,K}\lesssim|\mathbf{v}_{n}|_{1,K}.\end{split}$
(14)

Proving the lower bound (8). We integrate by parts, use the first equation in
(1), denote the $2\times 2$ identity matrix by $\mathbf{I}$, integrate by
parts again, and use the definition of $\boldsymbol{\Pi}^{0}_{\perp,p-3}$ in
(4), and deduce

$\begin{split}|\mathbf{v}_{n}|_{1,K}^{2}&=(\boldsymbol{\nabla}\mathbf{v}_{n},\boldsymbol{\nabla}\mathbf{v}_{n})_{0,K}=(\mathbf{v}_{n},(\boldsymbol{\nabla}\mathbf{v}_{n})\mathbf{n}_{K})_{0,\partial
K}+(\mathbf{v}_{n},\mathbf{x}^{\\!\perp}q_{p-3}+\nabla s)_{0,K}\\\\[4.0pt]
&=(\mathbf{v}_{n},(\boldsymbol{\nabla}\mathbf{v}_{n}+\mathbf{I}s)\mathbf{n}_{K})_{0,\partial
K}+(\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{v}_{n},\mathbf{x}^{\\!\perp}q_{p-3})_{0,K}-(\operatorname{div}\mathbf{v}_{n},s)_{0,K}.\end{split}$
(15)

Introduce
$\boldsymbol{\sigma}:=\boldsymbol{\nabla}\mathbf{v}_{n}+\mathbf{I}s$. Due to
the first equation in (1),
$\operatorname{div}\boldsymbol{\sigma}=-\mathbf{x}^{\\!\perp}q_{p-3}$.

Using (12) and (14) in (15) yields

$|\mathbf{v}_{n}|_{1,K}^{2}\lesssim\left(h_{K}^{-1}\|\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{v}_{n}\|_{0,K}+\|\operatorname{div}\mathbf{v}_{n}\|_{0,K}\right)|\mathbf{v}_{n}|_{1,K}+\|\mathbf{v}_{n}\|_{\frac{1}{2},\partial
K}\|(\boldsymbol{\sigma})\mathbf{n}_{K}\|_{-\frac{1}{2},\partial K}.$

Applying the divergence trace inequality, see [36, Section $3.5.2$], and a
polynomial inverse inequality on $\partial K$ (recall that $\mathbf{v}_{n}$ is
a piecewise polynomial over $\partial K$) gives

$\begin{split}|\mathbf{v}_{n}|_{1,K}^{2}&\lesssim\left(h_{K}^{-1}\|\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{v}_{n}\|_{0,K}+\|\operatorname{div}\mathbf{v}_{n}\|_{0,K}\right)|\mathbf{v}_{n}|_{1,K}\\!+\\!h_{K}^{-\frac{1}{2}}\|\mathbf{v}_{n}\|_{0,\partial
K}\left(\|\boldsymbol{\sigma}\|_{0,K}+h_{K}\|\operatorname{div}\boldsymbol{\sigma}\|_{0,K}\right)\\\\[4.0pt]
&\lesssim\left(h_{K}^{-1}\|\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{v}_{n}\|_{0,K}\\!+\\!\|\operatorname{div}\mathbf{v}_{n}\|_{0,K}\right)|\mathbf{v}_{n}|_{1,K}\\!+\\!h_{K}^{-\frac{1}{2}}\|\mathbf{v}_{n}\|_{0,\partial
K}\left(\|\boldsymbol{\sigma}\|_{0,K}\\!+\\!h_{K}\|\mathbf{x}^{\\!\perp}q_{p-3}\|_{0,K}\right)\\!.\end{split}$

On the other hand, the triangle inequality and (14) entail

$\|\boldsymbol{\sigma}\|_{0,K}\leq|\mathbf{v}_{n}|_{1,K}+\|s\|_{0,K}\lesssim|\mathbf{v}_{n}|_{1,K}.$

Combining the two above estimates and recalling (12) leads to (8).

Next, we prove the upper bound (9). Let $\mathbf{v}_{n}\in[H^{1}(K)]^{2}$ with
zero average over $\partial K$. We have to show an upper bound of the three
terms on the right-hand side of (7) in terms of $|\mathbf{v}_{n}|_{1,K}$.

As for
$h_{K}^{-2}\|\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{v}_{n}\|_{0,K}^{2}$, it
suffices to use the stability of the $L^{2}$ projector and the scaled Poincaré
inequality; we can estimate $\|\operatorname{div}\mathbf{v}_{n}\|_{0,K}^{2}$
by $|\mathbf{v}_{n}|_{1,K}$ directly; finally, we control
$h_{K}^{-1}\|\mathbf{v}_{n}\|_{0,\partial K}^{2}$ by using the trace and the
Poincaré inequality. ∎

######  Remark 2.

Let $\\{\text{dof}^{\text{B}}_{j}\\}$ denote the set of boundary DoFs in
$\mathbf{V}_{n}(K)$, i.e., the DoFs of type Dv1 and Dv2. Following, e.g., [14,
eq. $(2.14)$], the boundary contribution
$h_{K}^{-1}(\mathbf{u}_{n},\mathbf{v}_{n})_{0,\partial K}$ in (7) can be
replaced by the equivalent term

$\sum_{j}\text{dof}^{\text{B}}_{j}(\mathbf{u}_{n})\;\text{dof}^{\text{B}}_{j}(\mathbf{v}_{n}).$

Thus, the stabilization in (7) is spectrally equivalent to

$h_{K}^{-2}(\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{u}_{n},\boldsymbol{\Pi}^{0}_{\perp,p-3}\mathbf{v}_{n})_{0,K}+(\operatorname{div}\mathbf{u}_{n},\operatorname{div}\mathbf{v}_{n})_{0,K}+\sum_{j}\text{dof}^{\text{B}}_{j}(\mathbf{u}_{n})\text{dof}^{\text{B}}_{j}(\mathbf{v}_{n}).$

###  2.3 Stability estimates for the “dofi-dofi” stabilization

Here, we prove stability estimates for the classical “dofi-dofi” stabilization

$S^{K}_{D}(\mathbf{u}_{n},\mathbf{v}_{n}):=\sum_{j=1}^{\dim(\mathbf{V}_{n}(K))}\text{dof}_{j}(\mathbf{u}_{n})\text{dof}_{j}(\mathbf{v}_{n}).$
(16)

We recall the following technical result; see, e.g., [24, Lemma $4.1$].

######  Lemma 2.2.

Let $\\{\mathbf{m}_{\boldsymbol{\alpha}}\\}$ be a set of linearly independent
polynomials of maximum degree $p\in\mathbb{N}$ over a polygon $K$, that are
invariant with respect to translation and dilation. For every polynomial
$\mathbf{q}_{p}\in[\mathbb{P}_{p}(K)]^{d}$, $d=2,3$, consider the
decomposition
$\mathbf{q}_{p}=\sum_{\boldsymbol{\alpha}}(\overrightarrow{\mathbf{q}_{p}})_{\boldsymbol{\alpha}}\mathbf{m}_{\boldsymbol{\alpha}}$,
where $\overrightarrow{\mathbf{q}_{p}}$ denotes the coefficient vector of
$\mathbf{q}_{p}$ with respect to the basis
$\\{\mathbf{m}_{\boldsymbol{\alpha}}\\}$. Then, the following equivalence of
norms is valid:

$h_{K}\|\overrightarrow{\mathbf{q}_{p}}\|_{\ell^{2}}\lesssim\|\mathbf{q}_{p}\|_{0,K}\lesssim
h_{K}\|\overrightarrow{\mathbf{q}_{p}}\|_{\ell^{2}}.$ (17)

Under the mesh regularity assumption in Section 1, the hidden constant in (17)
are uniform with respect to the element $K$.

Next, we prove the following stability result, based on the techniques
developed in the proof of Theorem 2.1.

######  Theorem 2.3.

The following stability bounds are valid: there exist
$0<\alpha_{*}<\alpha^{*}$ independent of $h_{K}$ such that

$\displaystyle\alpha_{*}|\mathbf{v}_{n}|_{1,K}^{2}\leq
S^{K}_{D}(\mathbf{v}_{n},\mathbf{v}_{n})\qquad\forall\mathbf{v}_{n}\in\mathbf{V}_{n}(K),$
(18) $\displaystyle
S^{K}_{D}(\mathbf{v}_{n},\mathbf{v}_{n})\leq\alpha^{*}|\mathbf{v}_{n}|_{1,K}^{2}\qquad\forall\mathbf{v}_{n}\in\mathbf{V}_{n}(K)\text{
such that }\int_{\partial K}\mathbf{v}_{n}=\mathbf{0}.$ (19)

In particular, bounds (18) and (19) are valid for functions in
$\mathbf{V}_{n}(K)\cap\ker(\boldsymbol{\Pi}^{\nabla}_{p})$.

###### Proof.

We begin by proving the lower bound (18). Throughout, we use the same notation
as in the proof of Theorem 2.1. For the sake of presentation, we denote the
boundary ($\textbf{Dv}_{1}$ and $\textbf{Dv}_{2}$), “orthogonal”
($\textbf{Dv}_{3}$), and divergence ($\textbf{Dv}_{4}$) DoFs by
$\text{dof}^{\text{B}}(\cdot)$, $\text{dof}^{\perp}(\cdot)$, and
$\text{dof}^{\text{div}}(\cdot)$, respectively.

We have

$\begin{split}&|\mathbf{v}_{n}|_{1,K}^{2}\\\
&=-\int_{K}\boldsymbol{\Delta}\mathbf{v}_{n}\cdot\mathbf{v}_{n}+\int_{\partial
K}(\boldsymbol{\nabla}\mathbf{v}_{n})\mathbf{n}_{K}\cdot\mathbf{v}_{n}=\int_{K}(\nabla
s+\mathbf{x}^{\\!\perp}q_{p-3})\cdot\mathbf{v}_{n}+\int_{\partial
K}(\boldsymbol{\nabla}\mathbf{v}_{n})\mathbf{n}_{K}\cdot\mathbf{v}_{n}\\\
&=\int_{K}\mathbf{x}^{\\!\perp}q_{p-3}\cdot\mathbf{v}_{n}-\int_{K}s\
\operatorname{div}\mathbf{v}_{n}+\int_{\partial
K}(\mathbf{I}s+\boldsymbol{\nabla}\mathbf{v}_{n})\mathbf{n}_{K}\cdot\mathbf{v}_{n}=A+B+C.\end{split}$
(20)

We estimate the three terms on the right-hand side separately.

Denote $\mathbf{x}^{\\!\perp}q_{p-3}$ by $\mathbf{g}$ and consider the
expansion

$\mathbf{g}:=\sum_{\\!\boldsymbol{\alpha}}\overrightarrow{\mathbf{g}}_{\boldsymbol{\alpha}}\mathbf{q}_{\boldsymbol{\alpha}},$

where $\\{\mathbf{q}_{\boldsymbol{\alpha}}\\}$ is any basis of
$\mathbf{x}^{\\!\perp}\mathbb{P}_{p-3}(K)$ as in Lemma 2.2, and
$\overrightarrow{\mathbf{g}}$ is the vector of the coefficients of
$\mathbf{g}$ with respect to the basis
$\\{\mathbf{q}_{\boldsymbol{\alpha}}\\}$.

We obtain

$\begin{split}A&=\int_{K}\mathbf{x}^{\\!\perp}q_{p-3}\cdot\mathbf{v}_{n}=\sum_{\boldsymbol{\alpha}}\overrightarrow{\mathbf{g}}_{\\!\boldsymbol{\alpha}}(\mathbf{q}_{\boldsymbol{\alpha}},\mathbf{v}_{n})_{0,K}=\sum_{\boldsymbol{\alpha}}\overrightarrow{\mathbf{g}}_{\\!\boldsymbol{\alpha}}|K|\text{dof}^{\perp}_{\boldsymbol{\alpha}}(\mathbf{v}_{n})\\\
&\lesssim
h_{K}^{2}\|\overrightarrow{\mathbf{g}}\|_{\ell^{2}}\Big{(}\sum_{\boldsymbol{\alpha}}\text{dof}^{\perp}_{\boldsymbol{\alpha}}(\mathbf{v}_{n})^{2}\Big{)}^{\frac{1}{2}}\overset{\eqref{equivalence:polynomials:l2}}{\lesssim}h_{K}\|\mathbf{g}\|_{0,K}\Big{(}\sum_{\boldsymbol{\alpha}}\text{dof}^{\perp}_{\boldsymbol{\alpha}}(\mathbf{v}_{n})^{2}\Big{)}^{\frac{1}{2}}.\end{split}$

On the other hand, using the definition of $\mathbf{g}$, bound (12) gives

$h_{K}\|\mathbf{g}\|_{0,K}=h_{K}\|\mathbf{x}^{\\!\perp}q_{p-3}\|_{0,K}\lesssim|\mathbf{v}_{n}|_{1,K}.$

Combining the above bounds yields

$A\lesssim|\mathbf{v}_{n}|_{1,K}\Big{(}\sum_{\boldsymbol{\alpha}}\text{dof}^{\perp}_{\boldsymbol{\alpha}}(\mathbf{v}_{n})^{2}\Big{)}^{\frac{1}{2}}.$
(21)

Next, we focus on the term $B$. Recall that
$\operatorname{div}\mathbf{v}_{n}\in\mathbb{P}_{p-1}(K)\setminus\mathbb{R}$
and set

$\Pi^{0}_{p-1}s=:g=\sum_{\boldsymbol{\beta}}\overrightarrow{g}_{\\!\boldsymbol{\beta}}m_{\boldsymbol{\beta}},$

where $\\{m_{\boldsymbol{\beta}}\\}$ is any basis of
$\mathbb{P}_{p-1}(K)\setminus\mathbb{R}$ as in (the scalar version of) Lemma
2.2.

We deduce

$\begin{split}B&=\sum_{\boldsymbol{\beta}}\overrightarrow{g}_{\\!\boldsymbol{\beta}}\int_{K}m_{\boldsymbol{\beta}}\operatorname{div}\mathbf{v}_{n}\approx
h_{K}\sum_{\boldsymbol{\beta}}\overrightarrow{g}_{\\!\boldsymbol{\beta}}\text{dof}^{\text{div}}_{\boldsymbol{\beta}}(\mathbf{v}_{n})\leq
h_{K}\|\overrightarrow{g}\|_{\ell^{2}}\Big{(}\sum_{\boldsymbol{\beta}}\text{dof}^{\text{div}}_{\boldsymbol{\beta}}(\mathbf{v}_{n})^{2}\Big{)}^{\frac{1}{2}}\\\
&\overset{\eqref{equivalence:polynomials:l2}}{\lesssim}\|\Pi^{0}_{p-1}s\|_{0,K}\Big{(}\sum_{\boldsymbol{\beta}}\text{dof}^{\text{div}}_{\boldsymbol{\beta}}(\mathbf{v}_{n})^{2}\Big{)}^{\frac{1}{2}}\overset{\eqref{proof-
stab-
fact-3}}{\lesssim}|\mathbf{v}_{n}|_{1,K}\Big{(}\sum_{\boldsymbol{\beta}}\text{dof}^{\text{div}}_{\boldsymbol{\beta}}(\mathbf{v}_{n})^{2}\Big{)}^{\frac{1}{2}}.\end{split}$
(22)

The term $C$ can be estimated as in the proof of Theorem 2.1:

$C\lesssim h_{K}^{-\frac{1}{2}}\|\mathbf{v}_{n}\|_{0,\partial
K}|\mathbf{v}_{n}|_{1,K}.$ (23)

Inserting (21), (22), and (23) in (20), we deduce

$|\mathbf{v}_{n}|_{1,K}^{2}\lesssim\sum_{\boldsymbol{\alpha}}\text{dof}^{\perp}_{\boldsymbol{\alpha}}(\mathbf{v}_{n})^{2}+\sum_{\boldsymbol{\beta}}\text{dof}^{\text{div}}_{\boldsymbol{\beta}}(\mathbf{v}_{n})^{2}+h_{K}^{-1}\|\mathbf{v}_{n}\|_{0,\partial
K}^{2}.$

Finally, the boundary contribution is spectrally equivalent to the sum of the
boundary degrees of freedom squared; see Remark 2.

Next, we prove the upper bound (19). Notably, we need to estimate three types
of degrees of freedom. Note that Lemma 2.2 easily implies
$\|\mathbf{q}_{\boldsymbol{\alpha}}\|_{0,K}\lesssim h_{K}$ and
$\|m_{\boldsymbol{\alpha}}\|_{0,K}\lesssim h_{K}$. Since $\mathbf{v}_{n}$ has
zero average over $\partial K$ and we consider scaled polynomial functions in
the definition of the DoFs (2) and (3), we readily have a bound on the DoFs of
type $\textbf{Dv}_{3}$ and $\textbf{Dv}_{4}$:

$\frac{1}{|K|}\int_{K}\mathbf{v}_{n}\cdot\mathbf{q}_{\boldsymbol{\alpha}}\lesssim
h_{K}^{-2}\|\mathbf{v}_{n}\|_{0,K}\|\mathbf{q}_{\boldsymbol{\alpha}}\|_{0,K}\lesssim|\mathbf{v}_{n}|_{1,K}$

and

$\frac{h_{K}}{|K|}\int_{K}\operatorname{div}\mathbf{v}_{n}m_{\boldsymbol{\alpha}}\lesssim
h_{K}^{-1}\|\operatorname{div}\mathbf{v}_{n}\|_{0,K}\|m_{\boldsymbol{\alpha}}\|_{0,K}\lesssim|\mathbf{v}_{n}|_{1,K}.$

As for the boundary DoFs contribution $\textbf{Dv}_{1}$ and $\textbf{Dv}_{2}$,
we resort to Remark 2, write the sum of the boundary DoFs equivalently as
$\|\mathbf{v}_{n}\|_{0,\partial K}$, and use a trace inequality and a Poincaré
inequality. ∎

######  Remark 3.

Recalling the definition of the discrete bilinear form
$a_{h}^{K}(\cdot,\cdot)$ in (6), stability bounds involving
$a_{h}^{K}(\cdot,\cdot)$ instead of $S^{K}(\cdot,\cdot)$ are analogous to
those shown in Theorems 2.1 and 2.3 with $\alpha_{*}$ and $\alpha^{*}$
replaced by $\min(1,\alpha_{*})$ and $\max(1,\alpha^{*})$.

###  2.4 Interpolation estimates

Interpolation estimates for Stokes-type virtual element functions are well-
known; see [11, Proposition $4.2$] for the standard 2D case and [12, Theorem
$4.1$] for the enhanced 2D case. To the best of our knowledge, no explicit
interpolation estimates are available for 3D Stokes-type virtual element
functions.

Here, we prove interpolation estimates undertaking a different avenue, notably
using the stability estimates in Theorem 2.1. This novel approach is
interesting _per se_. In fact, it can be easily generalized to derive
interpolation properties for other virtual element spaces, once stability
estimates are available.

######  Theorem 2.4.

Given $\mathbf{u}\in[H^{s+1}(K)]^{2}$, $0<s\leq p$, and $\mathbf{u}_{I}$ its
DoFs interpolant, the following bound is valid:

$\|\mathbf{u}-\mathbf{u}_{I}\|_{0,K}+h_{K}|\mathbf{u}-\mathbf{u}_{I}|_{1,K}\lesssim
h_{K}^{s+1}|\mathbf{u}|_{s+1,K}.$

The hidden constant depends on the shape-regularity of the mesh and the degree
of accuracy $p$.

###### Proof.

Let $\mathbf{u}_{\pi}$ be $\boldsymbol{\Pi}^{\nabla}_{p}\mathbf{u}$, being
$\boldsymbol{\Pi}^{\nabla}_{p}$ defined in (5). The triangle inequality gives

$|\mathbf{u}-\mathbf{u}_{I}|_{1,K}\leq|\mathbf{u}_{\pi}-\mathbf{u}_{I}|_{1,K}+|\mathbf{u}-\mathbf{u}_{\pi}|_{1,K}.$

Thanks to standard polynomial approximation results, we only need to bound the
first term on the right-hand side. Observe that
$\mathbf{u}_{\pi}-\mathbf{u}_{I}$ belongs to $\mathbf{V}_{n}(K)$. Let
$S^{K}(\cdot,\cdot)$ be defined in (7). Then, we use bound (8) and write

$\begin{split}|\mathbf{u}_{\pi}-\mathbf{u}_{I}|_{1,K}^{2}&\lesssim
S^{K}(\mathbf{u}_{\pi}-\mathbf{u}_{I},\mathbf{u}_{\pi}-\mathbf{u}_{I})\\\\[4.0pt]
&=h_{K}^{-2}\|\boldsymbol{\Pi}^{0}_{\perp,p-3}(\mathbf{u}_{\pi}-\mathbf{u}_{I})\|_{0,K}^{2}+\|\operatorname{div}(\mathbf{u}_{\pi}-\mathbf{u}_{I})\|_{0,K}^{2}+\sum_{j}\text{dof}^{\text{B}}_{j}(\mathbf{u}_{\pi}-\mathbf{u}_{I})^{2}.\end{split}$

We prove that each of the three terms above can be bounded by an error term
involving the difference $\mathbf{u}\\!-\\!\mathbf{u}_{\pi}$. To this aim, we
preliminary observe that

$\boldsymbol{\Pi}^{0}_{\perp,p-3}(\mathbf{u}-\mathbf{u}_{I})=0,\qquad\operatorname{div}\mathbf{u}_{I}=\Pi^{0}_{p-1}(\operatorname{div}\mathbf{u}),\qquad\text{dof}^{\text{B}}_{j}(\mathbf{u}-\mathbf{u}_{I})=0,$
(24)

where $\Pi^{0}_{p-1}:L^{2}(K)\to\mathbb{P}_{p-1}(K)$ is the $L^{2}$ scalar
orthogonal projection onto $\mathbb{P}_{p-1}(K)$.

Thus, as for the bulk $L^{2}$ term, we use (24), the stability of orthogonal
projections, the fact that $\mathbf{u}-\mathbf{u}_{\pi}$ has zero average over
$K$, and the Poincaré inequality. Instead, we estimate the divergence
contribution by using (24) and bounding the $L^{2}$ norm of the divergence by
the $H^{1}$ seminorm. Eventually, the boundary contribution is bounded using
(24) again and the Sobolev embedding
$[H^{1+\varepsilon}(K)]^{2}\hookrightarrow[L^{\infty}(K)]^{2}$,
$\varepsilon>0$, as follows: for each boundary degree of freedom
$\text{dof}^{\text{B}}_{j}$,

$\text{dof}^{\text{B}}_{j}(\mathbf{u}_{\pi}-\mathbf{u}_{I})^{2}=\text{dof}^{\text{B}}_{j}(\mathbf{u}-\mathbf{u}_{\pi})^{2}\leq\|\mathbf{u}-\mathbf{u}_{\pi}\|_{\infty,K}^{2}\lesssim
h_{K}^{-2}\|\mathbf{u}-\mathbf{u}_{\pi}\|_{0,E}^{2}+h_{K}^{2\varepsilon}|\mathbf{u}-\mathbf{u}_{\pi}|_{1+\varepsilon,E}^{2}.$

Estimates in the $H^{1}$ norm follow using standard polynomial approximation
properties as in [41]. Instead, estimates in the $L^{2}$ norm are a simple
consequence of Poincaré-type arguments and the estimates in the $H^{1}$
seminorm. ∎

##  3 The three dimensional case

This section is devoted at proving stability and interpolation properties in
three dimensions. In Section 3.1, we recall the definition of the Stokes-like
virtual element space [9]. Stability properties are derived in Section 3.2 for
a projection based stabilization (comments on the degrees of freedom based
stabilization are discussed in Remark 4). In Section 3.3, we provide
interpolation estimates based on the previously proven stability properties.

###  3.1 Virtual element spaces in three dimensions

#### Virtual element spaces on faces.

Given a polyhedron $K\in\mathcal{T}_{n}$, on each of its faces $F$, we define
the $H^{1}$ projector
$\boldsymbol{\Pi}^{\nabla,F}_{p}:[H^{1}(F)]^{3}\to[\mathbb{P}_{p}(F)]^{3}$ as

$(\nabla_{F}\mathbf{q}_{p}^{F},\nabla_{F}(\mathbf{v}-\boldsymbol{\Pi}^{\nabla,F}_{p}\mathbf{v}))_{0,F}=0,\qquad\int_{\partial
F}\mathbf{v}=\mathbf{0}\qquad\forall\mathbf{v}\in[H^{1}(F)]^{3},\;\forall\mathbf{q}_{p}^{F}\in[\mathbb{P}_{p}(F)]^{3}.$
(25)

Based on this, we define the nodal (enhanced) virtual element space

$\mathbf{V}_{n}(F):=\\{\mathbf{v}_{n}\in[\mathcal{C}^{0}(F)]^{3}\mid\Delta_{F}\mathbf{v}_{n}\in[\mathbb{P}_{p}(F)]^{3},\;\mathbf{v}_{n}{}_{|e}\in[\mathbb{P}_{p}(e)]^{3}\;\forall
e\in\mathcal{E}^{F},\;\mathbf{v}_{n}\text{ satisfies }\eqref{enhancing-
constraint-face}\\},$ (26)

where, given $[\mathbb{P}_{p,p-2}(K)]^{3}$ the space of homogeneous vector
polynomials of degree larger than $p-2$ and smaller than or equal to $p$,

$\int_{F}(\mathbf{v}_{n}-\boldsymbol{\Pi}^{\nabla,F}_{p}\mathbf{v}_{n})\mathbf{q}_{p,p-2}^{F}=0\quad\forall\mathbf{q}_{p,p-2}^{F}\in[\mathbb{P}_{p,p-2}(K)]^{3}.$
(27)

We endow the space $\mathbf{V}_{n}(F)$ with the following set of unisolvent
DoFs [3]:

  * •

the vector values Dv${}_{1}^{F}$($\mathbf{v}_{n}$) at the vertices of $F$;

  * •

the vector values Dv${}_{2}^{F}$($\mathbf{v}_{n}$) at the $p-1$ internal Gauß-
Lobatto on each edge $e$ of $\mathcal{E}^{F}$;

  * •

for $p\geq 2$, given $\\{\mathbf{q}_{\boldsymbol{\alpha}}^{F}\\}$ a basis of
$[\mathbb{P}_{p-2}(F)]^{3}$, the moments

$\textbf{Dv}^{F}_{3}(\mathbf{v}_{n}):=\frac{1}{|F|}\int_{F}\mathbf{v}_{n}\cdot\mathbf{q}_{\boldsymbol{\alpha}}^{F}.$

Such DoFs allow for the computation of the projector
$\boldsymbol{\Pi}^{\nabla,F}_{p}$ in (25). Thus, see [5], the enhancing
constraint (27) allows for the computation of the orthogonal projector
$\boldsymbol{\Pi}^{0,F}_{p}:[L^{2}(F)]^{3}\to[\mathbb{P}_{p}(F)]^{3}$ defined
as

$\begin{split}(\mathbf{q}_{p}^{F},\mathbf{v}_{n}-\boldsymbol{\Pi}^{0,F}_{p}\mathbf{v}_{n})_{0,F}=0\qquad\forall\mathbf{v}_{n}\in\mathbf{V}_{n}(F),\quad\forall\mathbf{q}_{p}^{F}\in[\mathbb{P}_{p}(F)]^{3}.\end{split}$

#### Virtual element spaces on polyhedra.

Given a polyhedron $K$, we define the space

$\mathbf{V}_{n}(K):=\\{\mathbf{v}_{n}\in[H^{1}(K)]^{3}\mid\mathbf{v}_{n}\text{
satisfies }\eqref{local-auxiliary-problem-3D}\\},$

where, for some $s\in L^{2}_{0}(K)$,

$\begin{cases}-\boldsymbol{\Delta}\mathbf{v}_{n}-\nabla
s=\mathbf{x}\wedge\mathbf{q}_{p-3}&\mathbf{q}_{p-3}\in[\mathbb{P}_{p-3}(K)]^{3}\\\\[4.0pt]
\operatorname{div}\mathbf{v}_{n}=q_{p-1}&q_{p-1}\in\mathbb{P}_{p-1}(K)\\\\[4.0pt]
\mathbf{v}_{n}{}_{|\partial K}\in[\mathcal{C}^{0}(\partial
K)]^{3},\mathbf{v}_{n}{}_{|F}\in\mathbf{V}_{n}(F)&\forall
F\in\mathcal{E}^{K},\end{cases}$ (28)

with all the equations to be understood in a weak sense.

We endow the space $\mathbf{V}_{n}(K)$ with the following set of unisolvent
DoFs [9]:

  * •

the vector values Dv1($\mathbf{v}_{n}$) at the vertices of $K$;

  * •

the vector values Dv2($\mathbf{v}_{n}$) at the $p-1$ internal Gauß-Lobatto on
each edge $e$ of $K$;

  * •

for all faces $F\in\mathcal{E}^{F}$, given
$\\{\mathbf{q}_{\boldsymbol{\alpha}}^{F}\\}$ a basis of
$[\mathbb{P}_{p}(F)]^{3}$, the moments

$\textbf{Dv}^{F}_{3}(\mathbf{v}_{n}):=\frac{1}{|F|}\int_{F}\mathbf{v}_{n}\cdot\mathbf{q}_{\boldsymbol{\alpha}}^{F}.$
(29)

  * •

for $p\geq 3$, given $\\{\mathbf{q}_{\boldsymbol{\alpha}}\\}$ a basis of
$\mathbf{x}\wedge[\mathbb{P}_{p-3}(K)]^{3}$, the bulk “orthogonal” moments

$\textbf{Dv}_{4}(\mathbf{v}_{n}):=\frac{1}{|K|}\int_{K}\mathbf{v}_{n}\cdot\mathbf{q}_{\boldsymbol{\alpha}};$
(30)

  * •

given $\\{m_{\boldsymbol{\alpha}}\\}$ a basis of
$\mathbb{P}_{p-1}(K)\setminus\mathbb{R}$, the bulk “divergence” moments

$\textbf{Dv}_{5}(\mathbf{v}_{n}):=\frac{h_{K}}{|K|}\int_{K}\operatorname{div}\mathbf{v}_{n}\
m_{\boldsymbol{\alpha}}.$ (31)

As for the two dimensional case, we require that the bases
$\\{\mathbf{q}_{\boldsymbol{\alpha}}^{F}\\}$,
$\\{\mathbf{q}_{\boldsymbol{\alpha}}\\}$, and $\\{m_{\boldsymbol{\alpha}}\\}$
are invariant with respect to translations and dilations; see [3, 29].

The unisolvence of the above DoFs is proved, e.g., as in [9]. Thanks to the
enhancement in the definition of the virtual element spaces on faces, we can
compute the three dimensional version of the projectors in (4) and (5).

As in the two dimensional case, the global counterpart of the space
$\mathbf{V}_{n}(K)$ is constructed by a standard $H^{1}$-conforming DoFs-
coupling.

###  3.2 Stability estimates

The definition of the DoFs in Section 3.1 allows us to compute the orthogonal
projector
$\boldsymbol{\Pi}^{0}_{\wedge,p-3}:[L^{2}(K)]^{3}\to\mathbf{x}\wedge[\mathbb{P}_{p-3}(K)]^{3}$
defined as

$\begin{split}(\mathbf{q}_{\boldsymbol{\alpha}},\mathbf{v}_{n}-\boldsymbol{\Pi}^{0}_{\wedge,p-3}\mathbf{v}_{n})_{0,K}=0\qquad\forall\mathbf{v}_{n}\in\mathbf{V}_{n}(K),\quad\forall\mathbf{q}_{\boldsymbol{\alpha}}\in\mathbf{x}\wedge[\mathbb{P}_{p-3}(K)]^{3}.\end{split}$

We consider the local stabilization

$\begin{split}S^{K}(\mathbf{u}_{n},\mathbf{v}_{n})&:=h_{K}^{-2}(\boldsymbol{\Pi}^{0}_{\wedge,p-3}\mathbf{u}_{n},\boldsymbol{\Pi}^{0}_{\wedge,p-3}\mathbf{v}_{n})_{0,K}+(\operatorname{div}\mathbf{u}_{n},\operatorname{div}\mathbf{v}_{n})_{0,K}\\\\[4.0pt]
&\quad+\sum_{F\in\mathcal{E}^{K}}\left[h_{F}^{-1}(\boldsymbol{\Pi}^{0,F}_{p}\mathbf{u}_{n},\boldsymbol{\Pi}^{0,F}_{p}\mathbf{v}_{n})_{0,F}+(\mathbf{u}_{n},\mathbf{v}_{n})_{0,\partial
F}\right],\end{split}$ (32)

where we recall the polynomial bases in the definition of the DoFs (29), (30),
and (31), are defined so as to be invariant with respect to translation and
dilations. We prove the following stability estimates.

######  Theorem 3.1.

The following stability bounds are valid: there exists
$0<\alpha_{*}<\alpha^{*}$ independent of $h_{K}$ such that, for all
$\mathbf{v}_{n}$ in $\mathbf{V}_{n}(K)$ such that $\int_{\partial
K}\mathbf{v}_{n}=\mathbf{0}$,

$\displaystyle\alpha_{*}|\mathbf{v}_{n}|_{1,K}^{2}\leq
S^{K}(\mathbf{v}_{n},\mathbf{v}_{n}),$ (33) $\displaystyle
S^{K}(\mathbf{v}_{n},\mathbf{v}_{n})\leq\alpha^{*}|\mathbf{v}_{n}|_{1,K}^{2}.$
(34)

###### Proof.

We begin by proving the lower bound (33). As in the proof of Theorem 2.1, we
can prove the following bounds: for all $\mathbf{v}_{n}\in\mathbf{V}_{n}(K)$,
with “auxiliary pressure” $s$ and right-hand side
$\mathbf{x}\wedge\mathbf{q}_{p-3}$ in (28),

$|\boldsymbol{\Delta}\mathbf{v}_{n}|_{-1,K}\lesssim|\mathbf{v}_{n}|_{1,K},\qquad\|\mathbf{x}\wedge\mathbf{q}_{p-3}\|_{0,K}\lesssim
h_{K}^{-1}|\mathbf{v}_{n}|_{1,K},\qquad\|s\|_{0,K}\lesssim|\mathbf{v}_{n}|_{1,K}.$

With this at hand, as in the 2D case, we deduce

$|\mathbf{v}_{n}|_{1,K}\lesssim
h_{K}^{-1}\|\boldsymbol{\Pi}^{0}_{\wedge,p-3}\mathbf{v}_{n}\|_{0,K}+\|\operatorname{div}\mathbf{v}_{n}\|_{0,K}+\|\mathbf{v}_{n}\|_{\frac{1}{2},\partial
K}.$

Thus, we only have to estimate the last term on the right-hand side. After
recalling that $\mathbf{v}_{n}$ has zero vector average over $\partial K$, we
observe the obvious bound

$\|\mathbf{v}_{n}\|_{\frac{1}{2},\partial K}^{2}\lesssim
h_{K}\sum_{F\in\mathcal{E}^{K}}\|\nabla_{F}\mathbf{v}_{n}\|_{0,F}^{2}.$

We can estimate each face contribution by means of standard nodal virtual
element inverse estimates; see, e.g., [7, Theorem $2$] and [24, Section $3$]:

$\|\nabla_{F}\mathbf{v}_{n}\|_{0,F}^{2}\lesssim
h_{F}^{-2}\|\boldsymbol{\Pi}^{0,F}_{p}\mathbf{v}_{n}\|_{0,F}^{2}+h_{F}^{-1}\|\mathbf{v}_{n}\|_{0,\partial
F}^{2}.$

The assertion follows summing over all the faces and collecting the above
estimates.

Next, we focus on the upper bound (34). Notably, we estimate from above the
four terms on the right-hand side of (32). Since $\mathbf{v}_{n}$ has zero
average over $\partial K$, we readily have

$h_{K}^{-2}\|\boldsymbol{\Pi}^{0}_{\wedge,p-3}\mathbf{v}_{n}\|_{0,K}^{2}+\|\operatorname{div}\mathbf{v}_{n}\|_{0,K}^{2}\lesssim|\mathbf{v}_{n}|_{1,K}^{2}.$

The trace and the Poincaré inequality also yield

$\sum_{F\in\mathcal{E}^{K}}h_{F}^{-1}\|\boldsymbol{\Pi}^{0,F}_{p}\mathbf{v}_{n}\|_{0,F}^{2}\lesssim|\mathbf{v}_{n}|_{1,K}^{2}.$

We are left with estimating the fourth term on the right-hand side of (32).
Using a trace inequality on each face $F$ yields

$\|\mathbf{v}_{n}\|_{0,\partial F}\lesssim
h_{F}^{-\frac{1}{2}}\|\mathbf{v}_{n}\|_{0,F}+h_{F}^{\frac{1}{2}}|\mathbf{v}_{n}|_{1,F}.$

Inverse estimates for nodal virtual element functions, see, e.g., [24, Theorem
$3.6$], entail

$\|\mathbf{v}_{n}\|_{0,\partial F}\lesssim
h_{F}^{-\frac{1}{2}}\|\mathbf{v}_{n}\|_{0,F}.$

Taking the square on both sides and summing over the faces, then using another
trace inequality, and eventually a Poincaré inequality on $K$ give the
assertion. ∎

######  Remark 4.

Also the 3D “dofi-dofi” stabilization

$S^{K}_{D}(\mathbf{u}_{n},\mathbf{v}_{n}):=\sum_{j=1}^{\dim(\mathbf{V}_{n}(K))}\text{dof}_{j}(\mathbf{u}_{n})\text{dof}_{j}(\mathbf{v}_{n})$

satisfies the bounds (33) and (34). To see this it suffices to take the steps
from Theorem 3.1 and argue similarly as in the proof of Theorem 2.3. Moreover,
considerations analogous to those in Remark 3 are valid for the 3D case.

###  3.3 Interpolation estimates

Here, we prove interpolation estimates for 3D Stokes-like virtual element
spaces.

######  Theorem 3.2.

Given $\mathbf{u}\in[H^{s+1}(K)]^{3}$, $0<s\leq p$, and $\mathbf{u}_{I}$ its
DoFs interpolant, the following bound is valid:

$\|\mathbf{u}-\mathbf{u}_{I}\|_{0,K}+h_{K}|\mathbf{u}-\mathbf{u}_{I}|_{1,K}\lesssim
h_{K}^{s+1}|\mathbf{u}|_{s+1,K}.$

The hidden constant depends on the shape-regularity of the mesh and the degree
of accuracy $p$.

###### Proof.

Let $\mathbf{u}_{\pi}$ the best vector polynomial approximation of
$\mathbf{u}$ in $[H^{1}(K)]^{3}$; see, e.g., [41]. As in the proof of Theorem
2.4, we only need to estimate the energy of $\mathbf{u}_{\pi}-\mathbf{u}_{I}$.

In addition to the three dimensional counterpart of (24), we also have

$\boldsymbol{\Pi}^{0,F}_{p}(\mathbf{u}-\mathbf{u}_{I})=0\qquad\forall
F\in\mathcal{E}^{K}.$ (35)

Let $S^{K}(\cdot,\cdot)$ be defined in (32). Using (33), we write

$\begin{split}|\mathbf{u}_{\pi}-\mathbf{u}_{I}|_{1,K}^{2}&\lesssim
S^{K}(\mathbf{u}_{\pi}-\mathbf{u}_{I},\mathbf{u}_{\pi}-\mathbf{u}_{I})\\\
&=h_{K}^{-2}\|\boldsymbol{\Pi}^{0}_{\wedge,p-3}(\mathbf{u}_{\pi}-\mathbf{u}_{I})\|_{0,K}^{2}+\|\operatorname{div}(\mathbf{u}_{\pi}-\mathbf{u}_{I})\|_{0,K}^{2}\\\
&\quad+\sum_{F\in\mathcal{E}^{K}}\left[h_{F}^{-1}\|\boldsymbol{\Pi}^{0,F}_{p}(\mathbf{u}_{\pi}-\mathbf{u}_{I})\|_{0,F}^{2}\right]+\sum_{j}\text{dof}^{\text{B}}_{j}(\mathbf{u}_{\pi}-\mathbf{u}_{I})^{2}.\end{split}$

We have to bound the four terms on the right-hand side by
$|\mathbf{u}_{\pi}-\mathbf{u}_{I}|_{1,K}$. The first, second, and fourth are
dealt with by using (24), as in the proof of Theorem 2.4. As for the “new”
third term, it suffices to resort to (35), the continuity of
$\boldsymbol{\Pi}^{0,F}_{p}$, the trace inequality, and standard polynomial
approximation properties.

Estimates in the $L^{2}$ norm are a consequence of Poincaré type estimates on
$\mathbf{u}-\mathbf{u}_{I}$ and the energy estimates. ∎

######  Remark 5.

So far, we derived stability and interpolation properties for standard Stokes-
like spaces. Following, e.g., [12], we may also consider the enhanced version
of such spaces, in the same spirit as we defined the nodal virtual elements
spaces in (26). It is apparent that no essential modifications take place in
the proof of the stability estimates. The only difference between our setting
and the enhanced one is that in the latter we should employ polynomial inverse
estimates for slightly larger polynomial degrees.

##  4 Numerical validation of the stability estimates in 2D

In the foregoing sections, we proved stability bounds for fixed degree of
accuracy and regular polytopal meshes. Here, we numerically investigate the
behaviour of the stability constants from the practical side in two different
scenarios:

  * •

while keeping a mesh fixed, increase the degree of accuracy $p$;

  * •

while keeping the degree of accuracy $p$ fixed, consider sequences of elements
with degenerating geometry.

We focus on the two dimensional case, despite the arguments we discuss below
can be generalized to three dimensions.

Since we are interested in approximating the stability constants discussed in
Remark 3, we only need to investigate the behaviour of the minimum and maximum
generalized eigenvalues of problem

$\mathbf{A}\mathbf{v}=\lambda\mathbf{B}\mathbf{v},$ (36)

where the square symmetric matrices $\mathbf{A}$ and $\mathbf{B}$ are defined
as follows: for a given element $K\in\mathcal{T}_{n}$ and the canonical basis
$\\{\boldsymbol{\varphi}_{j}\\}$ of $\mathbf{V}_{n}(K)$, $j=1,\dots,N_{dof}$,
$N_{dof}$ being the dimension of $\mathbf{V}_{n}(K)$, associated with the
degrees of freedom Dv1($\cdot$), Dv2($\cdot$), Dv3($\cdot$), and Dv4($\cdot$),

$\mathbf{A}_{i,j}=a_{h}^{K}(\boldsymbol{\varphi}_{j},\boldsymbol{\varphi}_{i}),\qquad\mathbf{B}_{i,j}=(\boldsymbol{\nabla}\boldsymbol{\varphi}_{j},\boldsymbol{\nabla}\boldsymbol{\varphi}_{i})_{0,K}\qquad\forall
i,j=1,\dots,N_{dof}.$

The bilinear form $a_{h}^{K}(\cdot,\cdot)$ is computable via the degrees of
freedom following definition (6) and standard VEM arguments in order to
compute $\boldsymbol{\Pi}^{\nabla}_{p}$. Therefore, we only have to compute
the entries of $\mathbf{A}$. This is not immediate: the canonical basis
functions are not available in closed form since are solutions to local Stokes
problems with an unknown datum for the first equation. Thus, we need to detail
how to approximate them.

We split the canonical basis $\\{\boldsymbol{\varphi}\\}$ into three sets

$\\{\boldsymbol{\varphi}^{B}\\},\qquad\qquad\\{\boldsymbol{\varphi}^{\perp}\\},\qquad\qquad\\{\boldsymbol{\varphi}^{\text{div}}\\},$
(37)

which denote the basis functions associated with boundary degrees of freedom,
the “orthogonal” moments (2), and the divergence moments (3), respectively.

The elements of this basis are defined implicitly via the degrees of freedom,
so we cannot directly approximate them by means of any Galerkin methods
(unless resorting to some mixed formulation approach) For this reason, we
introduce a different basis $\\{\boldsymbol{\psi}\\}$ of the space
$\mathbf{V}_{n}(K)$ and split it into sets as those in (37):

$\\{\boldsymbol{\psi}^{B}\\},\qquad\qquad\\{\boldsymbol{\psi}^{\perp}\\},\qquad\qquad\\{\boldsymbol{\psi}^{\text{div}}\\}.$

The elements of this preliminary basis are constructed so that they solve
Stokes problems with given polynomial data, and can therefore be approximated
at any precision by means, e.g., of a finite element method on a
subtriangulation of the element.

We define the three type of basis functions as follows: given
$\\{m_{\boldsymbol{\alpha}}\\}$ a basis of $\mathbb{P}_{p-3}(K)$ and
$\\{m_{\boldsymbol{\gamma}}\\}$ a basis of
$\mathbb{P}_{p-1}(K)\setminus\mathbb{R}$ such that each
$m_{\boldsymbol{\gamma}}$ has zero average over $K$,

$\displaystyle\begin{cases}-\boldsymbol{\Delta}\boldsymbol{\psi}^{B}_{i}-\nabla
s=\mathbf{0}&\text{in }K\\\
\operatorname{div}\boldsymbol{\psi}^{B}_{i}=|K|^{-1}\int_{\partial
K}\boldsymbol{\varphi}^{B}_{i}\cdot\mathbf{n}_{K}&\text{in }K\\\
\boldsymbol{\psi}^{B}_{i}=\boldsymbol{\varphi}^{B}_{i}&\text{on }\partial
K\end{cases}$ $\displaystyle\qquad\forall i=1,\dots,2p\cdot(\\#\text{edges
of~{}$K$}),$
$\displaystyle\begin{cases}-\boldsymbol{\Delta}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\alpha}}-\nabla
s=\mathbf{x}^{\\!\perp}m_{\boldsymbol{\alpha}}&\text{in }K\\\
\operatorname{div}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\alpha}}=0&\text{in
}K\\\ \boldsymbol{\psi}^{\perp}_{\boldsymbol{\alpha}}=\mathbf{0}&\text{on
}\partial K\end{cases}$
$\displaystyle\qquad\forall|\boldsymbol{\alpha}|=0,\dots,p-3,$
$\displaystyle\begin{cases}-\boldsymbol{\Delta}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\gamma}}-\nabla
s=\mathbf{0}&\text{in }K\\\
\operatorname{div}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\gamma}}=m_{\boldsymbol{\gamma}}&\text{in
}K\\\ \boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\gamma}}=\mathbf{0}&\text{on
}\partial K\end{cases}$
$\displaystyle\qquad\forall|\boldsymbol{\gamma}|=1,\dots,p-1.$

It is possible to check that the span of the three above set of functions is
the space $\mathbf{V}_{n}(K)$. Therefore, we can write each basis function
$\boldsymbol{\varphi}$ as a linear combination of the $\boldsymbol{\psi}$
functions. Furthermore, the $\boldsymbol{\psi}$ functions can be approximated
at any precision by employing finite elements on sufficiently fine
triangulations of $K$. In other words, if we have at hand any finite element
approximation of the $\boldsymbol{\psi}$ basis functions, then we only have to
write the $\boldsymbol{\varphi}$ functions in terms of the $\boldsymbol{\psi}$
functions, and then compute the resulting matrix $\mathbf{B}$. This is what we
detail in Sections 4.1, 4.2, and 4.3. The algorithm reads as follows:

  1. 1.

approximate the $\boldsymbol{\psi}$ basis functions using a FEM triangulation
on $K$;

  2. 2.

find the $\boldsymbol{\varphi}$ basis functions as a linear combination of the
$\boldsymbol{\psi}$ basis functions, described in Sections 4.1, 4.2, and 4.3
below;

  3. 3.

compute the matrices $\mathbf{A}$ and $\mathbf{B}$;

  4. 4.

solve the generalized eigenvalue problem (36).

In Section 4.4, we check the behaviour of the minimum and maximum generalized
eigenvalues of (36) on a fixed element and increasing the degree of accuracy
$p$ of the scheme; instead, in Section 4.5, we keep fixed the degree of
accuracy, focus on two different elements, deform them in different ways, and
check the behaviour of the corresponding stability constants.

###  4.1 Expanding the boundary-type functions

For all $i=1,\dots,2p\cdot(\\#\text{edges of~{}$K$})$, we write

$\boldsymbol{\varphi}^{B}_{i}=\sum_{j}A^{(i)}_{j}\boldsymbol{\psi}^{B}_{j}+\sum_{\boldsymbol{\beta}}B^{(i)}_{\boldsymbol{\beta}}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}+\sum_{\boldsymbol{\delta}}C^{(i)}_{\boldsymbol{\delta}}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}.$

We have to determine the A, B, and C-type coefficients by imposing the DoFs
definition.

First, we use the boundary DoFs. For any vertex or Gauß-Lobatto node $N_{K}$
on any edge $e$ of $K$, we have

$\mathcal{\delta}_{i,k}=\boldsymbol{\varphi}^{B}_{i}(N_{k})=\sum_{j}A^{(i)}_{j}\boldsymbol{\psi}^{B}_{j}(N_{k})=A^{(i)}_{k}\qquad\forall
k=1,\dots,2p\cdot(\\#\text{edges of~{}$K$}).$

Therefore, we get the simplified expression

$\boldsymbol{\varphi}^{B}_{i}=\boldsymbol{\psi}^{B}_{i}+\sum_{\boldsymbol{\beta}}B^{(i)}_{\boldsymbol{\beta}}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}+\sum_{\boldsymbol{\delta}}C^{(i)}_{\boldsymbol{\delta}}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}.$

To find the B and C-type coefficients, we impose the divergence and
“orthogonal” DoFs definition. First, we impose the “orthogonal” moments
definition and write

$0=\int_{K}\boldsymbol{\varphi}^{B}_{i}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})=\int_{K}\boldsymbol{\psi}^{B}_{i}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})+\sum_{\boldsymbol{\beta}}B^{(i)}_{\boldsymbol{\beta}}\int_{K}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})+\sum_{\boldsymbol{\delta}}C^{(i)}_{\boldsymbol{\delta}}\int_{K}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}}),$

whence we deduce the conditions, for all
$|\widetilde{\boldsymbol{\beta}}|=0,\dots,p-3$,

$\sum_{\boldsymbol{\beta}}\left(\int_{K}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})\right)B^{(i)}_{\boldsymbol{\beta}}+\sum_{\boldsymbol{\delta}}\left(\int_{K}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})\right)C^{(i)}_{\boldsymbol{\delta}}=-\int_{K}\boldsymbol{\psi}^{B}_{i}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}}).$

Next, we impose the divergence moments definition, recall that the test
polynomial $m_{\widetilde{\boldsymbol{\delta}}}$ has zero average over $K$,
and get

$\begin{split}0&=\int_{K}\operatorname{div}\boldsymbol{\varphi}^{B}_{i}m_{\widetilde{\boldsymbol{\delta}}}=\int_{K}\operatorname{div}\boldsymbol{\psi}^{B}_{i}m_{\widetilde{\boldsymbol{\delta}}}+\sum_{\boldsymbol{\delta}}C^{(i)}_{\boldsymbol{\delta}}\int_{K}\operatorname{div}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}m_{\widetilde{\boldsymbol{\delta}}}\\\
&=|K|^{-1}\int_{\partial
K}\boldsymbol{\psi}^{B}_{i}\cdot\mathbf{n}_{K}\int_{K}m_{\widetilde{\boldsymbol{\delta}}}+\sum_{\boldsymbol{\delta}}C^{(i)}_{\boldsymbol{\delta}}\int_{K}m_{\boldsymbol{\delta}}\
m_{\widetilde{\boldsymbol{\delta}}}=\sum_{\boldsymbol{\delta}}C^{(i)}_{\boldsymbol{\delta}}\int_{K}m_{\boldsymbol{\delta}}\
m_{\widetilde{\boldsymbol{\delta}}},\end{split}$

whence we deduce the conditions, for all
$|\widetilde{\boldsymbol{\delta}}|=1,\dots,p-1$,

$\sum_{\boldsymbol{\delta}}\left(\int_{K}m_{\boldsymbol{\delta}}\
m_{\widetilde{\boldsymbol{\delta}}}\right)C^{(i)}_{\boldsymbol{\delta}}=0.$

Using the coercivity of any polynomial mass matrix, we deduce
$C^{(i)}_{\boldsymbol{\delta}}=0$ for all $i$ and $\boldsymbol{\delta}$.

Thus, the B-type coefficients are obtained by solving the linear system

$\sum_{\boldsymbol{\beta}}\left(\int_{K}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})\right)B^{(i)}_{\boldsymbol{\beta}}=-\int_{K}\boldsymbol{\psi}^{B}_{i}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}}).$

###  4.2 Expanding the orthogonal-type functions

For all $|\boldsymbol{\alpha}|=0,\dots,p-3$, we write

$\boldsymbol{\varphi}^{\perp}_{\boldsymbol{\alpha}}=\sum_{j}A^{(\boldsymbol{\alpha})}_{j}\boldsymbol{\psi}^{B}_{j}+\sum_{\boldsymbol{\beta}}B^{(\boldsymbol{\alpha})}_{\boldsymbol{\beta}}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}+\sum_{\boldsymbol{\delta}}C^{(\boldsymbol{\alpha})}_{\boldsymbol{\delta}}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}.$

We have to determine the A, B, and C-type coefficients by imposing the DoFs
definition. Imposing the boundary DoFs definition, we readily obtain that
$A^{(\boldsymbol{\alpha})}_{j}=0$ for all $j=1,\dots,2p\cdot(\\#\text{edges
of~{}$K$})$.

Thus, we focus on the other coefficients. First, we impose the “orthogonal”
moments definition and write

$\mathcal{\delta}_{\boldsymbol{\alpha},\widetilde{\boldsymbol{\beta}}}|K|=\int_{K}\boldsymbol{\varphi}^{\perp}_{\boldsymbol{\alpha}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})=\sum_{\boldsymbol{\beta}}B^{(\boldsymbol{\alpha})}_{\boldsymbol{\beta}}\int_{K}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})+\sum_{\boldsymbol{\delta}}C^{(\boldsymbol{\alpha})}_{\boldsymbol{\delta}}\int_{K}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}}),$

whence we deduce the conditions, for all
$|\widetilde{\boldsymbol{\beta}}|=0,\dots,p-3$,

$\sum_{\boldsymbol{\beta}}\left(\int_{K}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})\right)B^{(\boldsymbol{\alpha})}_{\boldsymbol{\beta}}+\sum_{\boldsymbol{\delta}}\left(\int_{K}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})\right)C^{(\boldsymbol{\alpha})}_{\boldsymbol{\delta}}=\mathcal{\delta}_{\boldsymbol{\alpha},\widetilde{\boldsymbol{\beta}}}|K|.$

Next, we impose the divergence moments definition and get

$0=\int_{K}\operatorname{div}\boldsymbol{\varphi}^{\perp}_{\boldsymbol{\alpha}}m_{\widetilde{\boldsymbol{\delta}}}=\sum_{\boldsymbol{\delta}}C^{(\boldsymbol{\alpha})}_{\boldsymbol{\delta}}\int_{K}\operatorname{div}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}m_{\widetilde{\boldsymbol{\delta}}}=\sum_{\boldsymbol{\delta}}C^{(\boldsymbol{\alpha})}_{\boldsymbol{\delta}}\int_{K}m_{\boldsymbol{\delta}}\
m_{\widetilde{\boldsymbol{\delta}}},$

whence we deduce the conditions, for all
$|\widetilde{\boldsymbol{\delta}}|=1,\dots,p-1$,

$\sum_{\boldsymbol{\delta}}\left(\int_{K}m_{\boldsymbol{\delta}}\
m_{\widetilde{\boldsymbol{\delta}}}\right)C^{(\boldsymbol{\alpha})}_{\boldsymbol{\delta}}=0.$

Using the coercivity of any polynomial mass matrix, we deduce
$C^{(\boldsymbol{\alpha})}_{\boldsymbol{\delta}}$ for all
$\boldsymbol{\alpha}$ and $\boldsymbol{\delta}$.

Thus, the B-type coefficients are obtained by solving the linear system

$\sum_{\boldsymbol{\beta}}\left(\int_{K}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})\right)B^{(\boldsymbol{\alpha})}_{\boldsymbol{\beta}}=\mathcal{\delta}_{\boldsymbol{\alpha},\widetilde{\boldsymbol{\beta}}}|K|.$

###  4.3 Expanding the divergence-type functions

For all $|\boldsymbol{\gamma}|=1,\dots,p-1$, we write

$\boldsymbol{\varphi}^{\text{div}}_{\boldsymbol{\gamma}}=\sum_{j}A^{(\boldsymbol{\gamma})}_{j}\boldsymbol{\psi}^{B}_{j}+\sum_{\boldsymbol{\beta}}B^{(\boldsymbol{\gamma})}_{\boldsymbol{\beta}}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}+\sum_{\boldsymbol{\delta}}C^{(\boldsymbol{\gamma})}_{\boldsymbol{\delta}}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}.$

We have to determine the A, B, and C-type coefficients by imposing the DoFs
definition. Imposing the boundary DoFs definition, we readily obtain that
$A^{(\boldsymbol{\gamma})}_{j}=0$ for all $j=1,\dots,2p\cdot(\\#\text{edges
of~{}$K$})$.

Thus, we focus on the other coefficients. First, we impose the “orthogonal”
moments definition and write

$0=\int_{K}\boldsymbol{\varphi}^{\text{div}}_{\boldsymbol{\gamma}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})=\sum_{\boldsymbol{\beta}}B^{(\boldsymbol{\gamma})}_{\boldsymbol{\beta}}\int_{K}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})+\sum_{\boldsymbol{\delta}}C^{(\boldsymbol{\gamma})}_{\boldsymbol{\delta}}\int_{K}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}}),$

whence we deduce the conditions, for all
$|\widetilde{\boldsymbol{\beta}}|=0,\dots,p-3$,

$\sum_{\boldsymbol{\beta}}\left(\int_{K}\boldsymbol{\psi}^{\perp}_{\boldsymbol{\beta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})\right)B^{(\boldsymbol{\gamma})}_{\boldsymbol{\gamma}}+\sum_{\boldsymbol{\delta}}\left(\int_{K}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}\cdot(\mathbf{x}^{\\!\perp}m_{\widetilde{\boldsymbol{\beta}}})\right)C^{(\boldsymbol{\gamma})}_{\boldsymbol{\delta}}=0.$
(38)

Next, we impose the divergence moments definition and get

$\frac{|K|}{h_{K}}\mathcal{\delta}_{\boldsymbol{\gamma},\widetilde{\boldsymbol{\delta}}}=\int_{K}\operatorname{div}\boldsymbol{\varphi}^{\text{div}}_{\boldsymbol{\alpha}}m_{\widetilde{\boldsymbol{\delta}}}=\sum_{\boldsymbol{\delta}}C^{(\boldsymbol{\gamma})}_{\boldsymbol{\delta}}\int_{K}\operatorname{div}\boldsymbol{\psi}^{\text{div}}_{\boldsymbol{\delta}}m_{\widetilde{\boldsymbol{\delta}}}=\sum_{\boldsymbol{\delta}}C^{(\boldsymbol{\gamma})}_{\boldsymbol{\delta}}\int_{K}m_{\boldsymbol{\delta}}\
m_{\widetilde{\boldsymbol{\delta}}},$

whence we deduce the conditions, for all
$|\widetilde{\boldsymbol{\delta}}|=1,\dots,p-1$,

$\sum_{\boldsymbol{\delta}}\left(\int_{K}m_{\boldsymbol{\delta}}\
m_{\widetilde{\boldsymbol{\delta}}}\right)C^{(\boldsymbol{\gamma})}_{\boldsymbol{\delta}}=\frac{|K|}{h_{K}}\mathcal{\delta}_{\boldsymbol{\gamma},\widetilde{\boldsymbol{\delta}}}.$
(39)

The B and C-type coefficients are obtained by solving the linear system
resulting from (38) and (39).

###  4.4 Stability constants when increasing the degree of accuracy

In light of the above approximation of the virtual element basis functions, we
provide here the minimum (nonzero) and maximum eigenvalues of the generalized
eigenvalue problem (36) on a given pentagon. Such constants clearly correspond
to $\lambda_{min}=c_{*}$ and $\lambda_{max}=c^{*}$ in the stability bound

$c_{*}a^{K}(\mathbf{v}_{n},\mathbf{v}_{n})\leq
a_{h}^{K}(\mathbf{v}_{n},\mathbf{v}_{n})\leq
c^{*}a^{K}(\mathbf{v}_{n},\mathbf{v}_{n})\qquad\forall\mathbf{v}_{n}\in\mathbf{V}_{n}(K)\setminus\mathbb{R}.$

Notably, we analyze the eigenvalues while increasing the degree of accuracy of
the scheme.

This is relevant to check as all the stability estimates are proved via
inverse estimates, which typically depend on the polynomial degree and the
degree of accuracy of the scheme. This is also interesting to check the value
of the “hidden” constants in the theoretical bounds. We employ the theoretical
(7) and the dofi-dofi (16) stabilizations and show the results in Table 1.

Table 1: Minimum and maximum eigenvalues of the generalized eigenvalue problem (36) employing the theoretical (7) and the dofi-dofi (16) stabilizations $S^{K}(\cdot,\cdot)$ and $S^{K}_{D}(\cdot,\cdot)$. We fix a pentagonal element and increase the degree of accuracy $p$. | $S^{K}(\cdot,\cdot)$ | $S_{D}^{K}(\cdot,\cdot)$
---|---|---
$p$ | $\lambda_{\min}$ | $\lambda_{\max}$ | $\lambda_{\min}$ | $\lambda_{\max}$
2 | 1.7265e-01 | 1.0037e+00 | 1.7098e-01 | 1.0037e+00
3 | 1.5570e-01 | 2.4381e+01 | 1.3578e-01 | 2.4269e+01
4 | 1.4606e-01 | 2.5338e+01 | 1.2649e-01 | 2.5206e+01
5 | 1.5364e-01 | 8.0281e+01 | 1.3667e-01 | 8.0154e+01
6 | 1.4909e-01 | 9.8691e+01 | 1.2966e-01 | 9.8519e+01
7 | 1.3999e-01 | 2.2056e+02 | 1.1612e-01 | 2.1990e+02
8 | 1.2007e-01 | 2.9917e+02 | 1.0388e-01 | 2.9877e+02
9 | 1.1133e-01 | 4.8713e+02 | 8.7618e-02 | 4.8637e+02

From Table 1, we observe that the stability constants depend on the degree of
accuracy $p$ only in a moderate way. This is not surprising as a similar
behaviour was observed for the Poisson-type virtual element method in [6,
Table 1] and [7, Table 1].

We are further interested in checking the condition number of the matrices
$\mathbf{A}$ and $\mathbf{B}$ appearing in (36), i.e., the condition number of
the approximate exact local stiffness matrix and its virtual element
counterpart; the details are detailed in Table 2 below.

Table 2: Condition numbers of the virtual element matrices $\mathbf{A}$ and $\mathbf{A}_{D}$ (computed with respect to the theoretical (7) and dofi-dofi (16) stabilizations) and the approximate exact stiffness matrix $\mathbf{B}$ for increasing degree of accuracy $p$. $p$ | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---
$\mathbf{A}$ | 8.76e+01 | 4.92e+03 | 4.84e+05 | 3.72e+07 | 2.35e+09 | 1.41e+11 | 8.32e+12 | 4.93e+14
$\mathbf{A}_{D}$ | 8.78e+01 | 4.91e+03 | 4.80e+05 | 3.66e+07 | 2.30e+09 | 1.37e+11 | 8.05e+12 | 4.75e+14
$\mathbf{B}$ | 9.43e+01 | 5.40e+03 | 5.75e+05 | 4.42e+07 | 2.79e+09 | 1.67e+11 | 9.80e+12 | 5.68e+14

From Table 2, we see that the condition number of the approximate exact
stiffness matrix is always larger than that of the corresponding virtual
element matrices. In a sense, this suggests that approximating the exact
stiffness matrix by a consistency and a stabilization term as in (6) is not
leading to larger condition numbers but rather has a slight beneficial effect.

###  4.5 Stability constants on sequences of badly shaped elements

Here, we provide the minimum (nonzero) and maximum eigenvalues of the
generalized eigenvalue problem (36) for a fixed degree of accuracy. Notably,
we analyze the behaviour of the eigenvalues on sequences of meshes with
degenerating geometry.

The first sequence we consider is constructed as in Figure 1: the first
element is a square with a hanging node; the other elements are obtained by
moving the hanging node towards the opposite edge. It is easy to check that
the star-shapedness constant of the sequence goes to zero.

$K_{1}$

$K_{2}$

$K_{3}$

$K_{4}$

$K_{5}$

Figure 1: First sequence of badly shaped elements. The first element is a
square with a hanging node; the other elements are obtained by moving the
hanging node towards the opposite edge.

Next, we consider a sequence of elements obtained starting with a regular
pentagon and halving the height of the previous element in the sequence; see
Figure 2.

$K_{1}$

$K_{2}$

$K_{3}$

$K_{4}$

$K_{5}$

Figure 2: First sequence of anisotropic elements. The first element is a
regular pentagon; the other elements are obtained by halving the height of the
previous polygon in the sequence.

Note that both sequences do not satisfy the regularity assumptions in Section
1. For this reason, we cannot guarantee theoretically that the stability
constants are robust with respect to the deformation of the elements.

In Table 3, we show the minimum (nonzero) and maximum generalized eigenvalues
of (36) for $p=3$, the theoretical (7) and the dofi-dofi (16) stabilizations,
and the sequences of meshes in Figures 1 and 2.

Table 3: Minimum and maximum eigenvalues of the generalized eigenvalue problem (36) employing the theoretical (7) and the dofi-dofi (16) stabilizations $S^{K}(\cdot,\cdot)$ and $S^{K}_{D}(\cdot,\cdot)$. We fix $p=3$ and consider the sequences of elements in Figures 1 and 2. | $S^{K}(\cdot,\cdot)$ | $S_{D}^{K}(\cdot,\cdot)$
---|---|---
$p=3$ | $\lambda_{\min}$ | $\lambda_{\max}$ | $\lambda_{\min}$ | $\lambda_{\max}$
Figure 1 | 1.7245e-01 | 2.4405e+01 | 1.5507e-01 | 2.4077e+01
| 2.4645e-02 | 2.6380e+01 | 2.2646e-02 | 2.5184e+01
| 2.0023e-02 | 5.2481e+01 | 1.9591e-02 | 5.0439e+01
| 1.1318e-02 | 8.3766e+01 | 1.1064e-02 | 8.0796e+01
| 6.2083e-03 | 1.1166e+02 | 6.0253e-03 | 1.0795e+02
Figure 2 | 1.1024e-01 | 2.9102e+01 | 1.0747e-01 | 2.9077e+01
| 3.5699e-02 | 5.7804e+01 | 3.5492e-02 | 5.7774e+01
| 8.6714e-03 | 1.7000e+02 | 8.6613e-03 | 1.6996e+02
| 1.9184e-03 | 6.1036e+02 | 1.9181e-03 | 6.1033e+02
| 4.9894e-04 | 2.3554e+03 | 4.9898e-04 | 2.3554e+03

From Table 3, we deduce that the stability constants are indeed deteriorating
together with the shape-regularity of the elements.

As in Section 4.4, we are interested in checking the condition number of the
matrices $\mathbf{A}$ and $\mathbf{B}$ appearing in (36), i.e., the condition
number of the approximate exact local stiffness matrix and its virtual element
counterpart; the results are given in Table 4 below.

Table 4: Condition numbers of the virtual element matrices $\mathbf{A}$ and $\mathbf{A}_{D}$ (computed with respect to the theoretical (7) and dofi-dofi (16) stabilizations) and the approximate exact stiffness matrix $\mathbf{B}$ for degree of accuracy $p=3$ on the two sequences of elements in Figures 1 and 2. $\mathbf{A}$ Fig. 1 | 2.9439e+03 | 1.3337e+04 | 1.8254e+04 | 2.1036e+04 | 2.2459e+04
---|---|---|---|---|---
$\mathbf{A}$ Fig. 2 | 2.2872e+04 | 3.1029e+05 | 4.6531e+06 | 7.2463e+07 | 1.1456e+09
$\mathbf{A}_{D}$ Fig. 1 | 2.9272e+03 | 1.3352e+04 | 1.8291e+04 | 2.1082e+04 | 2.2501e+04
$\mathbf{A}_{D}$ Fig. 2 | 2.2866e+04 | 3.1020e+05 | 4.6527e+06 | 7.2463e+07 | 1.1455e+09
$\mathbf{B}$ Fig. 1 | 3.6700e+03 | 3.2448e+04 | 1.9277e+05 | 3.2232e+05 | 4.1319e+05
$\mathbf{B}$ Fig. 2 | 2.4815e+04 | 3.2299e+05 | 9.1850e+06 | 4.0054e+08 | 1.4112e+10

From Table 4, the condition number of the approximate exact stiffness matrix
is always larger than that of the corresponding virtual element matrices.
Hence, similar comments as those for Table 2 apply.

##  5 Conclusions

We investigated some open issues in the analysis of Stokes-like virtual
element spaces. Notably, we proved stability properties in two and three
dimensions, and furthermore derived interpolation estimates by simplifying the
current state-of-the-art proofs. Numerical experiments seem to indicate that
the stability constants only depend moderately on $p$ but can degenerate more
rapidly for nonregular element geometries. On the other hand, such a
discrepancy of the discrete form with respect to the “exact” one may be
beneficial and explain why in this degenerate mesh conditions the VEM often
performs remarkably well.

#### Acknowledgments.

J. Meng has been supported by the China Scholarship Council (No.
$202106280167$) and the Fundamental Research Funds for the Central
Universities (No. xzy $022019040$). L. Beirão da Veiga was partially supported
by the Italian MIUR through the PRIN grants n. $905$ $201744$KLJL. L. Mascotto
acknowledges support from the Austrian Science Fund (FWF) Project P33477.

## References

  * [1] D. Adak, D. Mora, S. Natarajan, and A. Silgado.  A virtual element discretization for the time dependent Navier–Stokes equations in stream-function formulation.  ESAIM Math. Model. Numer. Anal., 55(5):2535–2566, 2021.
  * [2] J. Aghili, S. Boyaval, and D. A. Di Pietro.  Hybridization of mixed high-order methods on general meshes and application to the Stokes equations.  Comput. Methods Appl. Math., 15(2):111–134, 2015.
  * [3] 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.
  * [4] 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.
  * [5] 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(01):199–214, 2013.
  * [6] L. Beirão da Veiga, A. Chernov, L. Mascotto, and A. Russo.  Basic principles of $hp$ virtual elements on quasiuniform meshes.  Math. Models Methods Appl. Sci., 26(8):1567–1598, 2016.
  * [7] L. Beirão da Veiga, A. Chernov, L. Mascotto, and A. Russo.  Exponential convergence of the $hp$ virtual element method with corner singularity.  Numer. Math., 138(3):581–613, 2018.
  * [8] L. Beirão da Veiga, F. Dassi, G. Manzini, and L. Mascotto.  The virtual element method for the 3D resistive magnetohydrodynamic model.  https://arxiv.org/abs/2201.04417, 2022.
  * [9] L. Beirão da Veiga, F. Dassi, and G. Vacca.  The Stokes complex for virtual elements in three dimensions.  Math. Models Meth. Appl. Sci., 30(03):477–512, 2020.
  * [10] L. Beirão da Veiga, C. Lovadina, and A. Russo.  Stability analysis for the virtual element method.  Math. Models Methods Appl. Sci., 27(13):2557–2594, 2017.
  * [11] 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.
  * [12] L. Beirão da Veiga, C. Lovadina, and G. Vacca.  Virtual elements for the Navier–Stokes problem on polygonal meshes.  SIAM J. Numer. Anal., 56(3):1210–1242, 2018.
  * [13] L. Beirão da Veiga, D. Mora, and G. Vacca.  The Stokes complex for virtual elements with application to Navier–Stokes flows.  J. Sci. Comput., 81(2):990–1018, 2019.
  * [14] C. Bernardi and Y. Maday.  Polynomial interpolation results in Sobolev spaces.  J. Comput. Appl. Math., 43(1):53–80, 1992.
  * [15] D. Boffi, F. Brezzi, and M. Fortin.  Mixed Finite Element Methods and Applications, volume 44.  Springer Series in Computational Mathematics, 2013.
  * [16] L. Botti, D. A. Di Pietro, and J. Droniou.  A Hybrid High-Order method for the incompressible Navier–Stokes equations based on Temam’s device.  J. Comput. Phys., 376:786–816, 2019.
  * [17] S. C. Brenner and L.-Y. Sung.  Virtual element methods on meshes with small edges or faces.  Math. Models Methods Appl. Sci., 268(07):1291–1336, 2018.
  * [18] E. Burman, G. Delay, and A. Ern.  An unfitted hybrid high-order method for the Stokes interface problem.  IMA J. Numer. Anal., 41(4):2362–2387, 2021.
  * [19] E. Cáceres and G. N. Gatica.  A mixed virtual element method for the pseudostress-velocity formulation of the Stokes problem.  IMA J. Numer. Anal., 37(1):296–331, 2017.
  * [20] E. Cáceres, G. N. Gatica, and F. A. Sequeira.  A mixed virtual element method for the Brinkman problem.  Math. Models Meth. Appl. Sci., 27(04):707–743, 2017.
  * [21] E. Cáceres, G. N. Gatica, and F. A. Sequeira.  A mixed virtual element method for quasi-Newtonian Stokes flows.  SIAM J. Numer. Anal., 56(1):317–343, 2018.
  * [22] A. Cangiani, V. Gyrya, and G. Manzini.  The non-conforming virtual element method for the Stokes equations.  SIAM J. Numer. Anal., 54(6):3411–3435, 2016.
  * [23] S. Cao and L. Chen.  Anisotropic error estimates of the linear virtual element method on polygonal meshes.  SIAM J. Numer. Anal., 56(5):2913–2939, 2018.
  * [24] L. Chen and J. Huang.  Some error analysis on virtual element methods.  Calcolo, 55(1):1–23, 2018.
  * [25] A. Chernov, C. Marcati, and L. Mascotto.  $p$\- and $hp$-virtual elements for the Stokes problem.  Adv. Comp. Math., 47(24), 2021.
  * [26] B. Cockburn, G. Kanschat, D. Schötzau, and Ch. Schwab.  Local discontinuous Galerkin methods for the Stokes system.  SIAM J. Numer. Anal., 40(1):319–343, 2002.
  * [27] B. Cockburn and F.-J. Sayas.  Divergence-conforming HDG methods for Stokes flows.  Math. Comp., 83(288):1571–1598, 2014.
  * [28] B. Cockburn and K. Shi.  Devising HDG methods for Stokes flow: an overview.  Comput. & Fluids, 98:221–229, 2014.
  * [29] F. Dassi and G. Vacca.  Bricks for the mixed high-order virtual element method: Projectors and differential operators.  Appl. Numer. Math., 155:140–159, 2020.
  * [30] D. A. Di Pietro and S. Krell.  A hybrid high-order method for the steady incompressible Navier–Stokes problem.  J. Sci. Comput., 74(3):1677–1705, 2018.
  * [31] 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.
  * [32] G. N. Gatica, M. Munar, and F. A. Sequeira.  A mixed virtual element method for the Navier-Stokes equations.  Math. Models Methods Appl. Sci, 28(14):2719–2762, 2018.
  * [33] D. Irisarri and G. Hauke.  Stabilized virtual element methods for the unsteady incompressible Navier–Stokes equations.  Calcolo, 56(4):38, 2019.
  * [34] X. Liu and Z. Chen.  The nonconforming virtual element method for the Navier-Stokes equations.  Adv. Comput. Math., 45(1):51–74, 2019.
  * [35] X. Liu, J. Li, and Z. Chen.  A nonconforming virtual element method for the Stokes problem on general meshes.  Comput. Methods Appl. Mech. Engrg., 320:694–711, 2017.
  * [36] P. Monk.  Finite Element Methods for Maxwell’s Equations.  Oxford University Press, 2003.
  * [37] D. Mora, C. Reales, and A. Silgado.  A $\mathcal{C}^{1}$ virtual element method of high order for the Brinkman equations in stream function formulation with pressure recovery.  IMA J. Numer. Anal., 2021.  https://doi.org/10.1093/imanum/drab078.
  * [38] D. Mora and A. Silgado.  A $\mathcal{C}^{1}$ virtual element method for the stationary quasi-geostrophic equations of the ocean.  Comput. Math. Appl., 116:212–228, 2022.
  * [39] W. Qiu and K. Shi.  A superconvergent HDG method for the incompressible Navier–Stokes equations on general polyhedral meshes.  IMA J. Numer. Anal., 36(4):1943–1967, 2016.
  * [40] G. Vacca.  An ${H}^{1}$-conforming virtual element for Darcy and Brinkman equations.  Math. Models Methods Appl. Sci., 28(01):159–194, 2018.
  * [41] R. Verfürth.  A note on polynomial approximation in Sobolev spaces.  Math. Model. Numer. Anal., 33(4):715–719, 1999.
  * [42] G. Wang, F. Wang, L. Chen, and Y. He.  A divergence free weak virtual element method for the Stokes–Darcy problem on general meshes.  Comput. Methods Appl. Mech. Engrg., 344:998–1020, 2019.