The nonconforming Virtual Element Method for
the Stokes
equations
Abstract
We present the nonconforming Virtual Element Method (VEM) for the numerical approximation of velocity and pressure in the steady Stokes problem. The pressure is approximated using discontinuous piecewise polynomials, while each component of the velocity is approximated using the nonconforming virtual element space. On each mesh element the local virtual space contains the space of polynomials of up to a given degree, plus suitable nonpolynomial functions. The virtual element functions are implicitly defined as the solution of local Poisson problems with polynomial Neumann boundary conditions. As typical in VEM approaches, the explicit evaluation of the nonpolynomial functions is not required. This approach makes it possible to construct nonconforming (virtual) spaces for any polynomial degree regardless of the parity, for twoand threedimensional problems, and for meshes with very general polygonal and polyhedral elements. We show that the nonconforming VEM is infsup stable and establish optimal a priori error estimates for the velocity and pressure approximations. Numerical examples confirm the convergence analysis and the effectiveness of the method in providing highorder accurate approximations.
irtual element method, finite element method, polygonal and polyehdral mesh, highorder discretization, Stokes equations
65N30, 65N12, 65G99, 76R99
1 Introduction
We are concerned with the development of the nonconforming virtual element method (VEM) for the Stokes problem in the unknown fields and satisfying
(1)  
(2)  
(3) 
where is a polygonal or polyhedral domain in , with boundary . We will refer to and as velocity and pressure, respectively.
Historically, the first nonconforming finite element space dates back to the work of Crouzeix and Raviart in [27]. Their method provides a loworder accurate approximation of the velocity field of the Stokes equations on triangular meshes based on linear polynomials. Later on, higherorder accurate methods were proposed by Fortin and Soulie [29], and Crouzeix and Falk [26], respectively, by using finite element spaces based on polynomials of degree and on triangles. The functions in these finite element spaces are continuous on a discrete set of points located at the internal mesh edges. These points are the roots of the onedimensional order Legendre polynomials defined over the edges and can be used as the nodes of the GaussLegendre quadrature rule. This minimal continuity requirement ensures the optimal convergence rate; see, for instance, [27]. The construction of the nonconforming elements for the Stokes problem has been recently generalized on triangles in [37, 6] to consider polynomials of any degree , thus resulting in the socalled family of GaussLegendre nonconforming methods. Robust a posteriori estimates for such schemes can be found in [3]. A major drawback of the nonconforming GaussLegendre elements is that the space construction for even differs from that of odd . This feature also affects the classical loworder cases for , i.e., the formulation of the nonconforming spaces for and in [27, 26] is not the same as for in [29].
The generalization of the nonconforming formulation to elements other than triangles in two and three dimensions is quite a hard task due to the difficulty of the construction of the shape functions for such elements. For example, successful attempts in this direction are found for quadrilaterals, tetrahedra and hexahedra in [36, 22, 35, 34]. Instead, the construction of the nonconforming virtual element space for the Stokes equations that we present in this paper is straightforward for any polynomial degree regardless of its parity, and is the same for elements with very general geometric shape in two and three space dimensions.
The VEM was first introduced as a conforming formulation for the Poisson equation with constant coefficients in [7]. The nonconforming formulation for the same problem was developed later in [5, 23]. In both formulations, the trial and test functions are defined implicitly on each mesh element as the solution of a boundary value problem and never explicitly constructed in practice, hence the name “virtual”. In view of obtaining a computable and accurate virtual formulation, two essential ingredients are sought for the virtual element space: it must contain a space of polynomials up to a given degree; orthogonal and projections of virtual functions onto the polynomial subspace must be computable just using the degrees of freedom. Properties  make possible to avoid the explicit construction of the shape functions and allows us to formulate and implement the method on very general polygonal and polyhedral meshes. These features are inherited from the Mimetic Finite Difference (MFD) method [33, 12], which can be seen as a precursor. A conforming low order MFD method for the steady Stokes problem on general polygonal and polyhedral meshes is found in [9, 10], and a higher order MFD method equivalent to the nonconforming VEM in [5] is found in [32].
The objective of this paper is to develop the nonconforming VEM for the weak form of (1)(3) (see (4)(5) in the next section) that is suitable for very general mesh partitioning of in polygons and polyhedra. As is standard in the finite element setting, the VEM approximation of (1)(3) proposed in this work is based on the construction of a pair of finite element spaces satisfying the infsup condition, see [16]. The major features of this VEM are: each component of the velocity is locally approximated by the nonconforming virtual element space of order that contains the subspace of polynomials of degree at most and is globally nonconforming in the sense specified in Section 3, see also [5, 23]. The pressure is locally approximated by polynomials of degree at most , and is globally discontinuous; gradient and divergence are approximated by their projection onto polynomials of degree . Both projections are computable exactly using only the degrees of freedom of the VEM. Therefore, the divergencefree nature of the Stokes velocity is reproduced in the virtual framework by enforcing the divergencefree condition on the velocity approximation in a weak sense on each element. Moreover, the degree of the polynomials determines the accuracy (convergence rate) of the VEM; the wellposedness of the VEM is ensured through an additional stabilization term in the discrete weak formulation, which is computable using only the degrees of freedom of the VEM and is zero when applied to polynomials; the VEM allows for the use of mesh partitionings of with polygonal elements in 2D or polyhedral elements in 3D of arbitrary shape provided that a few typical shape regularity conditions are satisfied. The formulation of the method is the same in 2D and 3D and for any cell shape. It is worth mentioning that this feature follows from the nonconforming nature of the formulation, since the virtual conforming formulation, which also holds for general meshes, has to be constructed hierarchically in the space dimensions.
A number of relevant numerical approaches for the Stokes problem have been proposed in recent years. The VEM framework has already been applied to the streamline formulation of the Stokes equation in [4], a pseudostress velocity formulation can be found in [21], and a divergence free virtual approach can be found in [13]. Among the recent developments in discontinuous Galerkin methods, it is worth mentioning [20], the Hybridized Discontinuous Galerkin [25, 24], and, concerning polygonal meshes, the Hybrid High Order method [1] and the Weak Galerkin method [38]. A comparison between these different approaches is surely worth of investigation and will be the subject of further investigations.
The paper is organized as follows. In Section 2 we state the steady Stokes problem in weak form and introduce the VEM as a Galerkin method. In Section 3 we review the nonconforming virtual element framework used to approximate the velocity field, while implementation details can be found in Section 5. In Section 4 we prove the wellposedness and convergence of the VEM and we derive the error estimates for the velocity and pressure approximation. In Section 6 we numerically assess the performance of the VEM by solving a set of representative problems. In Section 7 we offer our final remarks and conclusions.
2 Continuous Stokes problem and discrete formulation
Find with on and such that for it holds:
(4)  
(5) 
where is the boundary of and the bilinear forms and are defined by:
(6) 
The wellposedness of (4)(5) follows from the coercivity of the form on the kernel of the form and the infsup condition [16].
In (4)(5) and throughout the paper we use the standard definitions and notation of Sobolev spaces, inner products, seminorms and norms. In particular, if is an open bounded domain with Lipschitz boundary in for and a nonnegative integer, denotes the standard Sobolev space of order ; is the associated inner product; and are the induced norm and seminorm, respectively. When as in (4)(5) we drop the subscripted symbol.
Let be a fixed integer. A Virtual Element Method of order will be defined by two finite dimensional functional spaces and of discrete trial velocity and pressure fields and bilinear forms and discrete counterparts of and , respectively. Precise definition of the functional spaces and and the construction of the bilinear forms and will be the focus of most of the remainder of this paper. For the moment, we only anticipate that we shall not assume the inclusion as our main goal is the development of a nonconforming approximation. Moreover, let be a suitable piecewise polynomial approximation of on the mesh partitioning of . The precise definition is given at the end of section 3.4. The virtual element formulation for the approximate solution of (4)(5) reads as:
Find such that
(7)  
(8) 
with and , respectively. The vector field in the righthand side integral of (7) is a suitable approximation of the vector field . The wellposedness of problem (7)(8) will follow from a discrete infsup condition which shall be established under suitable coercivity and stability properties introduced in the following section.
3 Virtual element framework
3.1 Mesh regularity and polynomial approximation
To ease the exposition, we assume that is a polygonal domain for and a polyhedral domain for . For any fixed we have a finite decomposition (the mesh) of the domain into nonoverlapping simple polygonal/polyhedral elements with maximum size . The adjective “simple” refers to the fact that the boundary of each element in the decomposition must be nonintersecting. Moreover, the boundary of element is made of a uniformly bounded number of interfaces (edges/faces), which are either part of the boundary of , or shared with another element of the decomposition. The definition of simple polygons and simple polyhedra is general enough to include, for instance, elements with consecutive coplanar edges/faces, such as those typical of locally refined meshes with hanging nodes and nonconvex elements.
Below, we use to denote a dimensional mesh interface (either an edge when or a face when ), to denote its length, to denote its unit normal vector with orientation fixed once and for all, and to denote the set of all such mesh interfaces in . When referring to the boundary of a specific element (with edges/faces) we use the notation and will have the outward orientation.
Assumption 1 (Mesh regularity)
We assume that there exists a constant such that:

for every element of and every interface , it holds that ;

every element of is starshaped with respect to a ball of radius ;

for , every face of the mesh is starshaped with respect to a ball of radius .
If is an internal edge/face of , then, there exist two elements and such that . Consider a scalar function defined on . We denote by the trace of on from within and by the unit vector orthogonal to and pointing out of . Then, the jump of the scalar function across is defined as . If, on the other hand, is on the domain boundary , then , with representing the trace of from within the element having as an interface and is the unit vector orthogonal to and pointing out of . Similarly, the jump of the vector quantity at the internal interface is given by
(9) 
and for the tensor quantity we may consider the jump operator that is such that
(10) 
for every properly sized tensor quantity .
We denote by for the orthogonal projection onto the polynomial space , defined for any function as the unique solution of the problem:
(11) 
For vector fields, i.e., , definition (11) is applied componentwise, thus giving the vector of polynomials . The well known approximation property of orthogonal projection are summarised in the following theorem.
Theorem 1 (Approximation using polynomials)
Under Assumption 1, the two following propositions hold true.

Let and let , for , denote the orthogonal projection onto the polynomial space . Then, for any , with , it holds that

Let be an interface shared by and let , for , denote the orthogonal projector onto the polynomial space . Then, for every , with , it holds
In both instances and , the positive constant depends only on the polynomial degree and the mesh regularity.
3.2 Discrete pressure space
As discrete trial space for pressures we use the standard space of piecewise polynomials of degree up to with respect to the domain partition :
The local degrees of freedom of for a pentagonal cell and the polynomial orders are illustrated by the right subpanels in Fig. 1. Note that by definition all functions in are with global zero mean. Therefore, if is the dimension of , the total number of degrees of freedom that are required for the pressure approximation is equal to .
3.3 Scalar nonconforming virtual element space
The scalar nonconforming virtual element space of order on the element is defined for as [5, 23]
(12) 
with the usual convention that . The virtual element space contains the space of polynomials of degree up to on . The complement is made up of functions that are deemed expensive to evaluate, although they can be represented in a discrete form through their degrees of freedom. The choice of the degrees of freedom is crucial to ensure that it is possible to define the bilinear forms in (7)(8) that are computable just using the degrees of freedom and the polynomial component of space . In practice, the discrete representation is sufficient for the construction of the method. To characterize the degrees of freedom of the functions in , we first introduce an appropriately scaled basis for . Denote by , , the set of scaled monomials
where is a multiindex and the center of gravity of . Furthermore, we define , a basis of the polynomial space whose size is . Bases for polynomial spaces defined on an interface can be similarly constructed; the same notation will be used.
The degrees of freedom for the scalar nonconforming space are [5, 23]:

for , the moments of degree on each edge/face :
(13) 
for , the moments of degree inside the element :
(14)
The degrees of freedom for a pentagonal cell and the polynomial orders are illustrated by the left subpanels in Fig. 1.
A counting argument shows that the cardinality of the above sets of degrees of freedom is , where we recall that denotes the number of edges/faces of element . Moreover, they are unisolvent in [5]. Indeed, if all the degrees of freedom of are zero we find that
(15) 
To prove (15), note that for it holds that ; for we have that and the first term on the right of (15) is a linear combination of the internal degrees of freedom of , which are zero by hypothesis. The second term on the right is also zero because it is a linear combination of the edge/face degrees of freedom of , which are zero by hypothesis. From it follows that is constant on and it must be zero since its value is equal to its zeroth order moment, which is zero by hypothesis.
The definition of the Virtual Element Method relies on the availability of elemental projection operators. The nonconforming VEM for elliptic problems introduced in [5] is based on the RitzGalerkin projection operator that for gives as the solution of the problem
together with the condition
This is shown in [5] to be computable for any using only the degrees of freedom (13) and (14). Furthermore, we shall prove in the following section that the projector of Theorem 1 is also computable when applied to first order derivatives of virtual functions. Together, these projectors will permit us to define a virtual formulation for the Stokes problem.
Remark 1
Instead, we note that an projection onto is not available. In view of definition (11), to compute the projection we would need the solution of the finite dimensional variational problem: find such that
(16) 
which needs the internal moments of up to order . As the corresponding degrees of freedom are available only up to , we need to resort to the strategy originally devised in [2] for the conforming virtual element spaces and extended to the nonconforming space by [23]. The availability of the projection becomes essential when loworder terms are present. Since in the present work we do not need this projection, we will not consider this issue anymore.
The global scalar nonconforming virtual element space is defined as a finite dimensional subspace of the nonconforming Sobolev space . The latter is a subspace of the broken Sobolev space
and, for , is given by
where is the jump of across the mesh interface defined in subsection 3.1. The global scalar nonconforming virtual element space of order is now given by
The degrees of freedom of are the edge/face moments (13) and the internal moments (14), and the size of is clearly given by , where is the number of cells in and the number of edges/faces in . Note that the edge/face moments are the same for the two mesh cells sharing a given internal edge/face. Therefore, the weak continuity condition on the jumps in the definition of is automatically satisfied.
3.4 Discrete velocity space
The discrete velocity space is defined by using the scalar nonconforming virtual element space of the previous subsection for each component. Hence,
and, similarly, . The degrees of freedom of are those inherited from each component, and are illustrated by the left subpanels in Fig. 1 for a pentagonal cell and the polynomial orders .
This choice of degrees of freedom ensures that the discrete differential operators
(17) 
when applied to a vector field in are computable using only the degrees of freedom of on . The projection is defined as
Integration by parts yields:
Likewise, the projection is defined as
Integration by parts yields:
The components of and are polynomials of degree at most in and, when restricted to each , are polynomials of degree at most . Therefore, the righthand side of the last equation above is computable using only the internal and the edge/face degrees of freedom of the (scalar) components .
3.5 Approximation of the bilinear forms and
In view of the following analysis, it is useful to extend the definition of the continuous bilinear forms and , to the whole of as a sum of the elemental contributions and ,
where and are defined by restricting the integrals in (6) to .
We define the approximate bilinear forms and used in (7)(8) by splitting them into local contributions
for any and , where and are bilinear forms on and , respectively.
The former bilinear form is defined by:
(18) 
where represents the RitzGalerkin projection operator introduced in Section 3.3 applied componentwisely. The term is the VEM stabilization term, cf. [7]. This can be any symmetric and coercive bilinear form satisfying
(19) 
for two positive constants and independent of and the mesh element .
Following [7, 5], in all computations presented in Section 6 we used the choice
(20) 
where is the vectorvalued linear operator that associates any virtual function with the vector of its th local degrees of freedom .
Remark 2
Again following [7, 5], we may have defined the consistency term on the righthand side of (18) as , but the approach used in (18) is more suitable to generalisations to problems with nonconstant coefficients, see, e.g., [23]. Furthermore, the availability of the projection for all proven in the previous section is of its own interest.
The second bilinear form is defined by:
(21) 
Remark 3
We note that in for any .
Definition (18) with the VEM stabilisation term satisfying (19) guarantees that the following polynomial consistency and stability properties are satisfied by the bilinear form .
Lemma 1 (Consistency and Stability)

Polynomial consistency: If or , or both, belong to , the bilinear form satisfies
(22) 
Stability: There exist two positive constants and independent of and the mesh element such that, for all , the bilinear form satisfies
(23)
Proof. Property is a straightforward consequence of the fact that the stabilization term is zero on polynomial vectors. To prove Property we first show that implies that is a constant vector, i.e., and have the same kernel. Indeed, consider such that . We find that
The coercivity of and property imply that , i.e., , and, thus, . Property implies that . Therefore, it holds that ; hence, is a constant vector. Now, (23) follows from (19) as in the scalar case, see [23] for details.
Remark 4
As is usual in the virtual element methodology, in order to satisfy conditions (22) and (23), the bilinear form is built as the sum of a “consistency” and a “stabilising” term, corresponding to the first and the second term in the righthand side of (18), respectively. The consistency term is exactly computable using only the degrees of freedom of and and satisfies the consistency condition (22). However, the consistency term alone does not satisfy the stability condition (23) on the whole virtual element space due to a rank deficiency of the operator, or equivalently, to the existence of a spurious kernel. Hence, a stabilization term must be added to fix this issue and actually remove the spurious kernel. This latter term is designed to vanish on the polynomial subspace not to affect the consistency property. As such, choice of the stabilization term is not unique [23]. For example, for the closely related MFD method, a number of studies investigated the optimal choice of the stabilization term with respect to a given criterion (reduction of dispersion effects, existence of a discrete maximum/minimum principle), cf. [31, 30, 17]. By exploiting the strict relation existing between the MFD method and the VEM, these alternative constructions of the stabilization term could be optionally considered in the present context.
3.6 Meshdependent energy norms
Hereafter, we shall use the energy seminorm on the broken Sobolev space :
A standard application of the results in [18] shows that a Poincaré inequality holds for the functions in , . Therefore, the seminorm is a norm in , and for this reason throughout the paper we prefer to use the notation and instead of and . To prove the infsup stability, we use the meshdependent energy seminorm on given by
(24) 
with
(25) 
Also, we will consider the affine subspace of vectorvalued functions where contains the scalar functions of whose trace on the boundary is equal to , i.e.,
(26) 
and the linear subspace obtained for . On we have the following equivalence of norms.
3.7 Approximation of the righthand side
We approximate the righthand side term by the linear functional
(28) 
Since is a polynomial of degree at most , each local linear functional is bounded and, for , it is computable by using the internal degrees of freedom of . However, it is not computable for . Indeed, the computation of above requires the knowledge of the average value of the components of on each element and such information is only available for . Therefore, the case deserves a special treatment. We follow Reference [5] and approximate by the average of the th order moments of associated with the edge/face of cell . Namely we consider
and note that is a firstorder approximation to , i.e., we have that
Then, we use to compute through the approximation:
Therefore, for we take:
We collect the results for the approximation of the righthand side functional for and in the following lemma. The proof follows from a straightforward extension of the scalar case, which is found in [5], to sized vectorvalued forcing terms and for this reason is omitted.
Lemma 3 (Approximation of the righthand side )
Let be integer numbers and consider , , and defined as above. There exists a constant independent of such that
Proof. See [5].
4 Error analysis
The wellposedness of the discrete problem (7)(8) is discussed in Section 4.1. The nonconformity error is estimated in Section 4.2. The convergence analysis is carried out in Sections 4.3 and 4.4, where error estimates for the approximation of the velocity and pressure fields are respectively derived. For the convergence analysis, we assume that on in (3) through the whole section.
Define the virtual interpolant of the vector field as the unique vector field whose degrees of freedom are the internal and edge/face moments of . Formally,

for , the degrees of freedom of associated with the edge/face are given by
(29) 
for , the degrees of freedom of associated with the mesh element are given by
(30)
The above relations must be interpreted componentwise. The unisolvence of the degrees of freedom implies the uniqueness of . Moreover, if all its moments on each edge/face on are zero and, consequently, belongs to . We also have the following result regarding the approximation of sufficiently smooth functions by the virtual interpolant, which may be proven as in [5].
Theorem 2 (Approximation using virtual element functions)
Let the nonconforming virtual element space of Section 3.4 for any integer , a positive integer such that , and a closed subset of . Under Assumption 1 (mesh regularity), for any , there exists an element such that
where is a positive constant that depends only on the polynomial degree and the mesh regularity constant .
4.1 Existence and uniqueness of the virtual element solution
The main result of this section is the existence and uniqueness of the virtual element solution , which is stated in Theorem 3. The proof of this theorem is based on the infsup property that is proven in the following lemma by adapting a classical argument.
Lemma 4 (Infsup)
There exists a strictly positive constant independent of such that for every in there exists a vector in such that
Proof. From [16] we know that there exists a strictly positive constant independent of such that for every in there exists a vector in such that
We can restrict this inequality to and for any consider the corresponding vector