The non-conforming Virtual Element Method forthe Stokes equations

The non-conforming Virtual Element Method for
the Stokes equations

Andrea Cangiani1    Vitaliy Gyrya2       Gianmarco Manzini2 3 {gyrya,gmanzini}
11 Department of Mathematics, University of Leicester, University Road - Leicester LE1 7RH, United Kingdom
22 Los Alamos National Laboratory, Theoretical Division, Group T-5, MS B284, Los Alamos, NM-87545, USA
33 Istituto di Matematica Applicata e Tecnologie Informatiche (IMATI) – CNR, via Ferrata 1, I – 27100 Pavia, Italy

We present the non-conforming 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 non-polynomial 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 non-polynomial functions is not required. This approach makes it possible to construct nonconforming (virtual) spaces for any polynomial degree regardless of the parity, for two-and three-dimensional problems, and for meshes with very general polygonal and polyhedral elements. We show that the non-conforming VEM is inf-sup 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 high-order accurate approximations.


irtual element method, finite element method, polygonal and polyehdral mesh, high-order discretization, Stokes equations


65N30, 65N12, 65G99, 76R99

1 Introduction

We are concerned with the development of the non-conforming virtual element method (VEM) for the Stokes problem in the unknown fields and satisfying


where is a polygonal or polyhedral domain in , with boundary . We will refer to and as velocity and pressure, respectively.

Historically, the first non-conforming finite element space dates back to the work of Crouzeix and Raviart in [27]. Their method provides a low-order accurate approximation of the velocity field of the Stokes equations on triangular meshes based on linear polynomials. Later on, higher-order 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 one-dimensional -order Legendre polynomials defined over the edges and can be used as the nodes of the Gauss-Legendre quadrature rule. This minimal continuity requirement ensures the optimal convergence rate; see, for instance, [27]. The construction of the non-conforming elements for the Stokes problem has been recently generalized on triangles in [37, 6] to consider polynomials of any degree , thus resulting in the so-called family of Gauss-Legendre non-conforming methods. Robust a posteriori estimates for such schemes can be found in [3]. A major drawback of the non-conforming Gauss-Legendre elements is that the space construction for even differs from that of odd . This feature also affects the classical low-order cases for , i.e., the formulation of the non-conforming spaces for and in [27, 26] is not the same as for in [29].

The generalization of the non-conforming 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 non-conforming 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 non-conforming 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 sub-space 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 non-conforming VEM in [5] is found in [32].

The objective of this paper is to develop the non-conforming 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 inf-sup condition, see [16]. The major features of this VEM are: each component of the velocity is locally approximated by the non-conforming virtual element space of order that contains the subspace of polynomials of degree at most and is globally non-conforming 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 divergence-free nature of the Stokes velocity is reproduced in the virtual framework by enforcing the divergence-free 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 well-posedness 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 pseudo-stress 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 non-conforming virtual element framework used to approximate the velocity field, while implementation details can be found in Section 5. In Section 4 we prove the well-posedness 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

The primary velocity-pressure formulation of the Stokes problem (1)-(3) takes the variational form:

Find with on and such that for it holds:


where is the boundary of and the bilinear forms and are defined by:


The well-posedness of (4)-(5) follows from the coercivity of the form on the kernel of the form and the inf-sup 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 non-negative 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 non-conforming 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


with and , respectively. The vector field in the right-hand side integral of (7) is a suitable approximation of the vector field . The well-posedness of problem (7)-(8) will follow from a discrete inf-sup 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 non-overlapping 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 non-intersecting. 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 co-planar edges/faces, such as those typical of locally refined meshes with hanging nodes and non-convex 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 star-shaped with respect to a ball of radius ;

  • for , every face of the mesh is star-shaped 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


and for the tensor quantity we may consider the jump operator that is such that


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:


For vector fields, i.e., , definition (11) is applied component-wise, 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.

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

  2. 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.

Proof.  This theorem can be proven using the theory in [19] for star-shaped domains and its extension to more general shaped elements presented in, e.g., [28].     

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 sub-panels 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 non-conforming virtual element space

The scalar non-conforming virtual element space of order on the element is defined for as [5, 23]


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 multi-index 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 non-conforming space are [5, 23]:

  • for , the moments of degree on each edge/face :

  • for , the moments of degree inside the element :


The degrees of freedom for a pentagonal cell and the polynomial orders are illustrated by the left sub-panels 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


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 zero-th order moment, which is zero by hypothesis.

The definition of the Virtual Element Method relies on the availability of elemental projection operators. The non-conforming VEM for elliptic problems introduced in [5] is based on the Ritz-Galerkin 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


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 non-conforming space by [23]. The availability of the -projection becomes essential when low-order terms are present. Since in the present work we do not need this projection, we will not consider this issue anymore.     

The global scalar non-conforming virtual element space is defined as a finite dimensional subspace of the non-conforming 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 non-conforming 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 non-conforming 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 sub-panels in Fig. 1 for a pentagonal cell and the polynomial orders .

Figure 1: Illustration of the degrees of freedom for the velocity and pressure solving the two-dimensional Stokes problem.

This choice of degrees of freedom ensures that the discrete differential operators


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 right-hand side of the last equation above is computable using only the internal and the edge/face degrees of freedom of the (scalar) components .

Finally, in the virtual element formulation (7)-(8), we use the virtual element space , whose definition requires the boundary function . This function is such that on every edge/face is the -orthogonal projection of on the polynomial space .

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:


where represents the Ritz-Galerkin projection operator introduced in Section 3.3 applied component-wisely. The term is the VEM stabilization term, cf. [7]. This can be any symmetric and coercive bilinear form satisfying


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


where is the vector-valued 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 right-hand side of (18) as , but the approach used in (18) is more suitable to generalisations to problems with non-constant 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:

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)

  1. Polynomial consistency: If or , or both, belong to , the bilinear form satisfies

  2. Stability: There exist two positive constants and independent of and the mesh element such that, for all , the bilinear form satisfies


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 right-hand 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 Mesh-dependent energy norms

Hereafter, we shall use the energy semi-norm 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 semi-norm is a norm in , and for this reason throughout the paper we prefer to use the notation and instead of and . To prove the inf-sup stability, we use the mesh-dependent energy seminorm on given by




Also, we will consider the affine subspace of vector-valued functions where contains the scalar functions of whose trace on the boundary is equal to , i.e.,


and the linear subspace obtained for . On we have the following equivalence of norms.

Lemma 2

The seminorm defined by (24)-(25) is a norm on and


Proof.  As noted in the proof of the stability condition (23), implies that is a constant vector, and, from the definition of , it follows that constant. Finally, from it follows that . The equivalence between the two norms in (27) is a straightforward consequence of (23).     

3.7 Approximation of the right-hand side

We approximate the right-hand side term by the linear functional


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 first-order 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 right-hand 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 vector-valued forcing terms and for this reason is omitted.

Lemma 3 (Approximation of the right-hand 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 well-posedness of the discrete problem (7)-(8) is discussed in Section 4.1. The non-conformity 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

  • for , the degrees of freedom of associated with the mesh element are given by


The above relations must be interpreted component-wise. 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 non-conforming 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 inf-sup property that is proven in the following lemma by adapting a classical argument.

Lemma 4 (Inf-sup)

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