{'fraction_non_alphanumeric': 0.04872766648619383, 'fraction_numerical': 0.01493412741382422, 'mean_word_length': 3.9358646030508853, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The non-conforming virtual element method for the Stokes equations

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.

1. Introduction. We are concerned with the development of the non-conforming virtual element method (VEM) for the Stokes problem in the unknown fields u and p satisfying −∆u + ∇p = f in Ω, (1.1) div u = 0 in Ω, (1.2) u = g on Γ, (1.3) where Ω is a polygonal or polyhedral domain in R d , d = 2, 3 with boundary Γ. We will refer to u and p 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, 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 k = 2 and 3 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 k th -order Legendre polynomials defined over the edges and can be used as the nodes of the k th 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 [6,37] to consider polynomials of any degree k > 3, 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 k differs from that of odd k. This feature also affects the classical low-order cases for k ≤ 3, i.e., the formulation of the non-conforming spaces for k = 1 and k = 3 in [26,27] is not the same as for k = 2 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 [22,34,35,36]. 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 C 0 -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: (i) it must contain a space of polynomials up to a given degree; (ii) orthogonal L 2 and H 1 projections of virtual functions onto the polynomial sub-space must be computable just using the degrees of freedom. Properties (i)-(ii) 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 [12,33], 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.1)-(1.3) (see (2.1)- (2.2) 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.1)- (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: (i) each component of the velocity is locally approximated by the non-conforming virtual element space of order k that contains the subspace of polynomials of degree at most k 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 k − 1, and is globally discontinuous; (ii) gradient and divergence are approximated by their projection onto polynomials of degree k − 1. 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; (iii) 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; (iv) 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 [24,25], 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.1)-(1.3) takes the variational form: where Γ is the boundary of Ω and the bilinear forms a and b are defined by: The well-posedness of (2.1)-(2.2) follows from the coercivity of the form a on the kernel of the form b and the inf-sup condition [16]. In (2.1)-(2.2) and throughout the paper we use the standard definitions and notation of Sobolev spaces, inner products, seminorms and norms. In particular, if D is an open bounded domain with Lipschitz boundary in R d for d = 2, 3 and m a non-negative integer, H m (D) denotes the standard Sobolev space of order m; (·, ·) m,D is the associated inner product; · m,D and | · | m,D are the induced norm and seminorm, respectively. When D = Ω as in (2.1)-(2.2) we drop the subscripted symbol.
Let k ≥ 1 be a fixed integer. A Virtual Element Method of order k will be defined by two finite dimensional functional spaces V k h and Φ k−1 h of discrete trial velocity and pressure fields and bilinear forms a h : → R discrete counterparts of a and b, respectively. Precise definition of the functional spaces V k h and Φ k h and the construction of the bilinear forms a h and b h 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 V k h ⊂ H 1 (Ω) d as our main goal is the development of a non-conforming approximation. Moreover, let g h be a suitable piecewise polynomial approximation of g 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 (2.1)-(2.2) reads as: The vector field f h in the right-hand side integral of (2.4) is a suitable approximation of the vector field f . The well-posedness of problem (2.4)-(2.5) will follow from a discrete inf-sup condition which shall be established under suitable coercivity and stability properties introduced in the following section.

Virtual element framework.
3.1. Mesh regularity and polynomial approximation. To ease the exposition, we assume that Ω is a polygonal domain for d = 2 and a polyhedral domain for d = 3. For any fixed h > 0 we have a finite decomposition (the mesh) T h of the domain Ω into non-overlapping simple polygonal/polyhedral elements with maximum size h. The adjective "simple" refers to the fact that the boundary of each element in the decomposition must be non-intersecting. Moreover, the boundary ∂E of element E 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 s to denote a d − 1 dimensional mesh interface (either an edge when d = 2 or a face when d = 3), |s| to denote its length, n s to denote its unit normal vector with orientation fixed once and for all, and E h to denote the set of all such mesh interfaces in T h . When referring to the boundary of a specific element E (with ν E edges/faces) we use the notation s ∈ ∂E and n s will have the outward orientation.
Assumption 1 (Mesh regularity). We assume that there exists a constant ρ > 0 such that: • for every element E of T h and every interface s ∈ E, it holds that h s ≥ ρh E ; • every element E of T h is star-shaped with respect to a ball of radius ρh E ; • for d = 3, every face s of the mesh is star-shaped with respect to a ball of radius ρh s .
If s is an internal edge/face of T h , then, there exist two elements E + and E − such that s ⊂ ∂E + ∩ ∂E − . Consider a scalar function v defined on Ω. We denote by v ± the trace of v | E ± on s from within E ± and by n ± s the unit vector orthogonal to s and pointing out of E ± . Then, the jump of the scalar function v across s is defined as [[v]] := v + n + s + v − n − s . If, on the other hand, s is on the domain boundary Γ, then [[v]] := vn s , with v representing the trace of v from within the element E having s as an interface and n s is the unit vector orthogonal to s and pointing out of Ω. Similarly, the jump of the vector quantity v h at the internal interface s is given by and for the tensor quantity nv h we may consider the jump operator [[ · ]] that is such that for every properly sized tensor quantity τ .
We denote by Π E l : L 2 (E) → P l (E) for l ≥ 0 the L 2 (E)-orthogonal projection onto the polynomial space P l (E), defined for any function v ∈ L 2 (E) as the unique solution of the problem: Theorem 1 (Approximation using polynomials). Under Assumption 1, the two following propositions hold true.
(i) Let E ∈ T h and let Π E l : L 2 (E) → P l (E), for l ≥ 0, denote the L 2 (E)-orthogonal projection onto the polynomial space P l (E). Then, for any w ∈ H m (E), with 1 ≤ m ≤ l + 1, it holds that (ii) Let s be an interface shared by E + , E − ∈ T h and let Π E l : L 2 (s) → P l (s), for l ≥ 0, denote the L 2 (s)-orthogonal projector onto the polynomial space P l (s). Then, for every In both instances (i) and (ii), the positive constant C depends only on the polynomial degree l 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 k − 1 with respect to the domain partition T h : The local degrees of freedom of Φ k−1 h for a pentagonal cell and the polynomial orders k = 1, . . . , 4 are illustrated by the right sub-panels in Fig. 3.1. Note that by definition all functions in Φ k−1 h are with global zero mean. Therefore, if N E d,k−1 is the dimension of P k−1 (E), the total number of degrees of freedom that are required for the pressure approximation is equal to E∈T h N E d,k−1 − 1. 3.3. Scalar non-conforming virtual element space. The scalar non-conforming virtual element space of order k ≥ 1 on the element E is defined for d = 2, 3 as [5,23] with the usual convention that 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 (2.4)-(2.5) that are computable just using the degrees of freedom and the polynomial component of space V k h (E). In practice, the discrete representation is sufficient for the construction of the method. To characterize the degrees of freedom of the functions in V k h (E), we first introduce an appropriately scaled basis for P k (E). Denote by M l (E), l ∈ N, the set of scaled monomials where α is a multi-index and x E the center of gravity of E. Furthermore, we define M k (E) := α=1 , a basis of the polynomial space P k (E) whose size is N d,k . Bases for polynomial spaces defined on an interface s can be similarly constructed; the same notation will be used.
The degrees of freedom for the scalar non-conforming space V k h (E) are [5,23]: • for k ≥ 1, the moments of degree (k − 1) on each edge/face s ∈ ∂E: • for k ≥ 2, the moments of degree (k − 2) inside the element E: The degrees of freedom for a pentagonal cell and the polynomial orders k = 1, . . . , 4 are illustrated by the left sub-panels in Fig. 3.1.
A counting argument shows that the cardinality of the above sets of degrees of freedom is where we recall that ν E denotes the number of edges/faces of element E. Moreover, they are unisolvent in V k h (E) [5]. Indeed, if all the degrees of freedom of v h are zero we find that To prove (3.7), note that for k = 1 it holds that ∆v h = 0; for k > 1 we have that ∆v h ∈ P k−2 (E) and the first term on the right of (3.7) is a linear combination of the internal degrees of freedom of v h , 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 v h , which are zero by hypothesis. From ∇v h 0,E = 0 it follows that v h is constant on E 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 Π ∇ k : together with the condition This is shown in [5] to be computable for any v h ∈ V k E using only the degrees of freedom (3.5) and (3.6). Furthermore, we shall prove in the following section that the L 2 -projector Π E k−1 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 which needs the internal moments of v h up to order k. As the corresponding degrees of freedom are available only up to k−2, 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 L 2 -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 H 1,nc k (T h ). The latter is a subspace of the broken Sobolev space and, for k ≥ 1, is given by ] is the jump of v across the mesh interface s ∈ E h defined in subsection 3.1. The global scalar non-conforming virtual element space of order k ≥ 1 is now given by The degrees of freedom of V k h are the edge/face moments (3.5) and the internal moments (3.6), and the size of V k h is clearly given by where N T h is the number of cells in T h and N E h the number of edges/faces in E h . 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 H 1,nc k (T h ) 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, The degrees of freedom of V k h are those inherited from each component, and are illustrated by the left sub-panels in Fig. 3.1 for a pentagonal cell and the polynomial orders k = 1, . . . , 4. This choice of degrees of freedom ensures that the discrete differential operators when applied to a vector field v h in V k h (E) are computable using only the degrees of freedom of v h on E. The projection Π E k−1 ∇v h is defined as Integration by parts yields: Integration by parts yields: in E and, when restricted to each s ∈ ∂E, are polynomials of degree at most k − 1. 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 v h .
Finally, in the virtual element formulation (2.4)-(2.5), we use the virtual element space V k h,g , whose definition requires the boundary function g h . This function is such that g h|s on every edge/face s ∈ Γ is the L 2 -orthogonal projection of g on the polynomial space P k−1 (s).
3.5. Approximation of the bilinear forms a(·, ·) and b(·, ·). In view of the following analysis, it is useful to extend the definition of the continuous bilinear forms a and b, to the whole of H 1 (T h ) as a sum of the elemental contributions a E and b E , where a E and b E are defined by restricting the integrals in (2.3) to E. We define the approximate bilinear forms a h and b h used in (2.4)-(2.5) by splitting them into local contributions The former bilinear form is defined by: where Π ∇ k represents the Ritz-Galerkin projection operator introduced in Section 3.3 applied component-wisely. The term S E h is the VEM stabilization term, cf. [7]. This can be any symmetric and coercive bilinear form satisfying for two positive constants c * and c * independent of h and the mesh element E. Following [5,7], in all computations presented in Section 6 we used the choice where χ i is the vector-valued linear operator that associates any virtual function φ with the vector of its i-th local degrees of freedom χ i (φ) ∈ R d . Remark 2. Again following [5,7], we may have defined the consistency term on the righthand side of (3.10) as E ∇Π ∇ k u h : ∇Π ∇ k v h dx, but the approach used in (3.10) is more suitable to generalisations to problems with non-constant coefficients, see, e.g., [23]. Furthermore, the availability of the projection Π E k−1 (∇v h ) for all v h ∈ V k h proven in the previous section is of its own interest.
The second bilinear form is defined by: Definition (3.10) with the VEM stabilisation term satisfying (3.11) guarantees that the following polynomial consistency and stability properties are satisfied by the bilinear form a E h . Lemma 1 (Consistency and Stability). (3.14) (ii) Stability: There exist two positive constants α * and α * independent of h and the mesh element E such that, for all v h ∈ V k h (E), the bilinear form a E h satisfies Proof. Property (i) is a straightforward consequence of the fact that the stabilization term is zero on polynomial vectors. To prove Property (ii) we first show that a E h (v h , v h ) = 0 implies that v h is a constant vector, i.e., a E and a E h have the same kernel. Indeed, The coercivity of S E h and property (b) imply that v h|E = Π E k (v h ), i.e., v h|E ∈ P k (E) d , and, thus, Property (a) implies that Π E k−1 (∇v h ) = 0. Therefore, it holds that ∇v h|E = 0; hence, v h|E is a constant vector. Now, (3.15) follows from (3.11) as in the scalar case, see [23] for details.
Remark 4. As is usual in the virtual element methodology, in order to satisfy conditions (3.14) and (3.15), the bilinear form a E h (u h , v h ) 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 (3.10), respectively.
The consistency term is exactly computable using only the degrees of freedom of u h and v h and satisfies the consistency condition (3.14). However, the consistency term alone does not satisfy the stability condition (3.15) 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. [17,30,31]. 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 H 1 (T h ) d : A standard application of the results in [18] shows that a Poincaré inequality holds for the functions in H 1,nc k (T h ), k ≥ 1. Therefore, the semi-norm |·| 1,h is a norm in H 1,nc k (T h ) d , and for this reason throughout the paper we prefer to use the notation v 1,h and v 1,E instead of |v| 1,h and |v| 1,E .
To prove the inf-sup stability, we use the mesh-dependent energy seminorm on V k h given by Also, we will consider the affine subspace of vector-valued functions contains the scalar functions of V k h whose trace on the boundary Γ is equal to g, i.e., 18) and the linear subspace V k h,0 obtained for g = 0. On V k h,0 we have the following equivalence of norms.
Proof. As noted in the proof of the stability condition (3.15), |||v h ||| E = 0 implies that v h|E is a constant vector, and, from the definition of H 1,nc The equivalence between the two norms in (3.19) is a straightforward consequence of (3.15).

3.7.
Approximation of the right-hand side (f , ·). We approximate the right-hand side term (f , · ) by the linear functional Since Π E max(k−2,0) (f ) is a polynomial of degree at most k −2, each local linear functional is bounded and, for k ≥ 2, it is computable by using the internal degrees of freedom of v h . However, it is not computable for k = 1. Indeed, the computation of (f h , v h ) E above requires the knowledge of the average value of the components of v h on each element E and such information is only available for k ≥ 2. Therefore, the case k = 1 deserves a special treatment. We follow Reference [5] and approximate Π E 0 (v h ) by the average of the 0-th order moments of v h associated with the edge/face of cell E. Namely we consider Then, we use v h to compute (f h , v h ) through the approximation: Therefore, for k = 1 we take: We collect the results for the approximation of the right-hand side functional (f , ·) for k = 1 and k > 1 in the following lemma. The proof follows from a straightforward extension of the scalar case, which is found in [5], to d-sized vector-valued forcing terms f and for this reason is omitted. Lemma 3 (Approximation of the right-hand side (f , ·)). Let m, k ≥ 1 be integer numbers and consider f ∈ H m−1 (Ω) d , r = min(k, m), and (f h , ·) ∈ (V k h ) d defined as above. There exists a constant C independent of h such that Proof. See [5]. • for k ≥ 2, the degrees of freedom of v I associated with the mesh element E are given by The above relations must be interpreted component-wise. The unisolvence of the degrees of freedom implies the uniqueness of v I . Moreover, if v ∈ H 1 0 (Ω) d all its moments on each edge/face s on Γ are zero and, consequently, v I belongs to V k h,0 . We also have the following result regarding the approximation of sufficiently smooth functions by the virtual interpolant, which may be proven as in [5].
where C is a positive constant that depends only on the polynomial degree k and the mesh regularity constant ρ.

Existence and uniqueness of the virtual element solution.
The main result of this section is the existence and uniqueness of the virtual element solution (u h , p h ) ∈ V k h × Φ k−1 h , 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 h such that for Proof. From [16] we know that there exists a strictly positive constant β independent of h such that for every q in L 2 (Ω)/R there exists a vectorv q in We can restrict this inequality to Φ k−1 h ⊂ L 2 (Ω)/R and for any q h ∈ Φ k−1 h consider the corresponding vectorv q . We will prove that wherev I q is the virtual interpolation ofv q defined in the previous section and γ = α * / √ α * . Such properties easily imply that from which the assertion of the lemma follows with v q =v I q and β = β/γ. To prove (i), consider the following development that starts from the definition of Property (i) readily follows from the elemental decomposition of b h and b.
To prove (ii), first note that the stability condition of a E h implies that Now, sincev I q ∈ V k h (E) it holds that ∆v I q is a vector of polynomials of degree k − 2 inside E and (n s · ∇)v I q is a vector of polynomials of degree k − 1 along each edge s ∈ ∂E and we have that: [use the left inequality of (4.3)]

the last inequality and (4.3) implies that |||v
and property (ii) follows from the continuity ofv q in H 1 0 (Ω) d by summing over all E ∈ T h and Theorem 3. The solution of the discrete problem (2.4)-(2.5) exists and is unique.
In view of the stability condition (3.15), the coercivity of a h (·, ·) on ker(b h ) follows from the coercivity of a(·, ·) on ker(b). Therefore, existence and uniqueness of the solution follow from the inf-sup property proved in Lemma 4, cf. [16].

4.2.
Estimate of the non-conformity error. The non-conformity error is controlled as in Lemma 6. The proof of this lemma requires a bound on the jumps of v h and nv h (which, we recall, are defined in (3.1) and (3.2)). This bound is provided by Lemma 5.
for the pressure and for the velocity. Using the Cauchy-Schwartz inequality and then applying component-wise the approximation estimates of Theorem 1 to bound each of the resulting terms, we obtain, cf. [5] or [27], where for each side s the symbols E + and E − denote the two elements sharing that side. As the number of edges/faces is assumed to be uniformly bounded, the required result follows with a positive constant C independent of h and u.
Lemma 6 (non-conformity error). Let k, m ≥ 1, r = min(k, m) be integer numbers. Let u ∈ H 1 0 (Ω) ∩ H m+1 (Ω) d and p ∈ L 2 (Ω)/R ∩ H m (Ω) be the velocity and pressure solution of problem (2.1)-(2.2), with source term f ∈ L 2 (Ω) d and homogeneous boundary condition g = 0 on Γ. Then, for every v h ∈ V k h,0 , it holds that and we can apply the Green's identity. By rearranging the summation on the internal sides and using the definition of the jump operators given in section 3.1, we obtain: The assertion of the lemma follows by applying the result of Lemma 5.

Error estimate for the velocity. Let
Using these definitions, problem (2.1)-(2.2) is equivalent to [27]: and problem (2.4)-(2.5) is equivalent to [27]: Theorem 4 (H 1 abstract a priori error bound for the velocity). Let u ∈ W be the solution of problem (4.6), and u h ∈ W k h the solution of problem (4.7) with k ≥ 1. Then, it holds that: The last term in the right-hand side of the above error estimate measures the non-conformity error, i.e. it is non-zero because W k h is a non-conforming space. Proof. Let v h be an arbitrary element of W k h and let w h = u h − v h . From the stability property of the virtual bilinear form (3.15) and equation (4.7) it follows that Then, we add and subtract (f , w h ) and a(u, w h ) and we obtain: The first term on the right-hand side characterizes the approximation of the source term f by f h ; the second term is the conformity error determined by choosing the test function w h in W k h instead of W; the third term simultaneously characterizes the approximation of u by v h and a(·, ·) by a h (·, ·) in W k h . To separate this linked dependence, we reformulate the last term as the summation of local contributions; then, we add and subtract a E (q, w h ) and a E h (q, w h ) to each summation argument, where q is a generic vector-valued function in Φ k h d , and we obtain: Since q |E is a polynomial vector of degree k, the intermediate term above is zero due to the polynomial consistency relation (3.14). Furthermore, we transform the last term above by using the continuity and stability of the virtual bilinear form: Hence, for all v h , w h ∈ W k h with w h = 0 and q ∈ Φ k h d , using the last inequality and the continuity of a E (·, ·), multiplying and dividing by w h 1,h , and using the triangle inequality The result now follows by applying the estimates above to the triangle inequality properly taking the supremum upper bound of the terms with w h and the infimum over the arbitrary vector functions q and v h . .
with f h|E defined as in section 3.7. Denote by u h ∈ W k h the virtual element solution to problem (4.7) under Assumption 1 (mesh regularity), where W k h is the divergence-free non-conforming virtual element space of vector-valued functions defined in (4.5). Then, there exists a constant C independent of h such that Proof. Since u is also solution of problem (4.6), we start the proof of the theorem by separately bounding the terms of the abstract bound of Theorem 4. The first term on the right-hand side of (4.8), i.e. inf v h ∈W k h u − v h 1,h , is easily bounded by introducing the virtual interpolant u I ∈ V k h of u as in Theorem 2 and noting that u I belongs to W k h . Indeed, for any element E, we have Π E k−1 (div u I ) ∈ P k−1 (E) and for any polynomial q ∈ P k−1 (E) we have that: Hence, Π E k−1 (div u I ) = 0. Likewise, the second term is bounded by using in each cell E the results of Theorem 1 since The bound on the third term is given by Lemma 3 since w h ∈ W k h ⊂ V k h . The bound on the fourth term is given by Lemma 6 and noting that b(w h , p) = 0 for every w h ∈ W k h . Finally, the assertion of the theorem follows by combining the bounds derived above and noting that the inequality constant C may depend only on the stability constants α * and α * , and the mesh regularity constant ρ.
Remark 5. Optimal order error estimates for the velocity approximation in the L 2 -norm can be derived by duality arguments, see always [27]. However, a more accurate approximation of the forcing terms than that provided by (3.20) would be needed for the case k = 1, 2. This can be obtained by following the approach in [8].

4.4.
Error estimate for the pressure. Let p I denote the piecewise polynomial function that is defined on each elements E of mesh T h by the orthogonal projection of p on the space of polynomials of degree k − 1; formally, (p I ) |E = Π E k−1 (p) ∈ P k−1 (E). The accuracy of this approximation is characterized by Theorem 1. To ease the notation we will also use the symbol p I E to denote the restriction (p I ) |E .
Theorem 6 (Abstract L 2 a priori error bound for the pressure). Let k ≥ 1 be an integer number. Let u ∈ H 1 0 (Ω) ∩ H 2 (Ω) d and p ∈ L 2 (Ω)/R ∩ H 1 (Ω) be the exact velocity and pressure solution to problem (2.1)-(2.2) with the homogeneous boundary condition g = 0 on Γ. Denote the projection of (4.10) and, in view of the norm equivalence (3.19), it holds that We use (2.4), and add and subtract (f , v h ) and a(u, v h ) to obtain: (4.12) We test (1.1) against v h ∈ V k h,0 , apply the Green's identity, rearrange the summation on edges/faces in E, and introduce the jump notation (cf. section 3.1) to obtain: Also, from Remark 3 we know that b h (v h , p I ) = b(v h , p I ), and using this relation and (4.13) in (4.12) we obtain: (4.14) The second term on the right-hand side can be further bounded by introducing a generic vector q ∈ [Φ k h ] d and reasoning as in the proof of Theorem 4, thus yielding The assertion of the theorem follows from (4.11) by taking the absolute value of (4.14) with (4.15) v h 1,h and the supremum on v h and q.
Theorem 7 (L 2 a priori error bound for the pressure). Under the assumptions and notations of Theorem 5, let u h ∈ V k h,0 and p h ∈ Φ k−1 h be the virtual element velocity and solution to problem (2.4)-(2.5). Then, there exists a constant C > 0 depending only on the stability constants α * and α * , the inf-sup constant β and the mesh regularity constant ρ such that Proof. The proof is just a matter of bounding the terms of the abstract bound of Theorem 6. Bounds for the first, second, and third term are already given in Theorem 5. The bound for the term containing b(v h , p − p I ) follows from the polynomial approximation results of Theorem 1.
The last two terms are bounded by using Lemma 5.
5. Implementation details. According to section 3.4, the key component of the VEM implementation is the construction of the projection operators Π E k−1 • ∇, Π E k−1 • div, and Π ∇ k . Once these projection operators are constructed, all terms in the local bilinear forms can be computed as integrals of polynomials just as in the standard FEM with the only exception of the stabilization term in a E h . This latter term does not require any integration but is directly defined through the action of the projector operator Π ∇ k on the degrees of freedom. The matrix representation of Π E k−1 • ∇ and Π E k−1 • div can be derived from the corresponding representation of the scalar projection operator Π E k−1 • ∂ ∂xi , i = 1, . . . , d, which has already been worked out in [23]. Similarly, we refer to [2] for the details on the computation of Π ∇ k . In the rest of this section we show how to compute the two terms of the bilinear form a E h , i.e., the consistency and stability term, cf. Remark 4, assuming a matrix representation of the projectors is given. The implementation formulas for the bilinear form b E h and the right-hand side linear functional (f h , ·) E can be derived similarly and are not shown here.
Consider the Lagrangian basis {φ i } for the scalar virtual element space V k h (E) associated with the degrees of freedom introduced in Section 3.3. We collect the coefficients of the expansion of each monomials m α on the basis φ i for α = 1, . . . , N d,k on the columns of matrix D, so that m α = N E j=1 φ j D j,α . Similarly, we collect the coefficients of the expansions of the polynomials Π ∇ k φ i and Π E k−1 ∂φ i /∂x i with respect to the monomial basis m α on the columns of matrices Π ∇ k and Π E,x l k−1 , respectively. Hence, these projections can be expressed by the formulas where ( ) i,j denotes the i, j-th element of a given matrix argument . Then, for the local discrete velocity space V k h (E), we consider the set of basis functions {φ l i } that are such that (φ l i ) l = φ i if l = l and (φ l i ) l = 0 otherwise. The generic entry of the consistency term of a E h in (3.10) is given by where H is the matrix with coefficients H α,β = E m α m β dx for α, β = 1, . . . , N d,k−1 . These coefficients may be computed exactly in special cases or by applying a sufficiently accurate integration rule in general. Similarly, the generic entry of the stabilisation term of a E h , which is defined in (3.12), is given by The numerical experiments that we present in this section are aimed at confirming the a priori analysis developed in the previous sections. In a preliminary stage, the consistency of non-conforming VEM, i.e. the exactness for polynomial solutions, has been tested numerically by solving the Stokes equation with boundary and source data determined by u(x, y) = (y m , x m ) and p = x m + y m on different set of polygonal meshes and for m = 1 to 4. In all the cases, we measure an error whose magnitude is of the order of the arithmetic precision, thus confirming this property.
To study the accuracy of the method we solve the problem with the following solution on the domain Ω =]0, 1[×]0, 1[: p(x, y) = sin(2πx) sin(2πy), (6.1) with f (x) = x 5 e −x . The forcing term and the Dirichlet boundary condition are set in accordance with (6.1). The performance of the VEM for k = 1, 2, 3, 4 is investigated by evaluating the rate of convergence on three different sequences of five meshes, labeled by M 1 , M 2 and M 3 , respectively. The top panels of Fig. 6.1 show the first mesh of each sequence and the bottom panels show the mesh of the first refinement. The meshes in M 1 are built by partitioning the domain Ω into square cells and relocating each interior node to a random position inside a square box centered at that node. The sides of this square box are aligned with the coordinate axis and their length is equal to 0.8 times the minimum distance between two adjacent nodes of the initial square mesh. The meshes in M 2 are built as follows. First, we determine a primal mesh by remapping the position ( x, y) of the nodes of a uniform square partition of Ω by the smooth coordinate transformation: x = x + (1/10) sin(2π x) sin(2π y), y = y + (1/10) sin(2π x) sin(2π y).

Randomised quadrilaterals
Remapped hexagons Non-convex octagons  Table 6.1 Mesh data for the meshes in M 1 , M 2 , and M 3 ; N E , Ne and Nv are the numbers of mesh elements, interfaces and vertices, respectively, and h is the mesh size parameter.
The corresponding mesh of M 2 is built from the primal mesh by splitting each quadrilateral cell into two triangles and connecting the barycenters of adjacent triangular cells by a straight segment. The mesh construction is completed at the boundary by connecting the barycenters of the triangular cells close to the boundary to the midpoints of the boundary edges and these latters to the boundary vertices of the primal mesh. The meshes in M 3 are obtained by filling the unit square with a suitably scaled non-convex octagonal reference cell.
All the meshes are parametrised by the number of partitions in each direction. The starting mesh of every sequence is built from a 5 × 5 regular grid, and the refined meshes are obtained by doubling this resolution. Mesh data for each refinement level, i.e., numbers of mesh elements, number of edges, number of vertices, are reported in Table 6.1. More details on these mesh constructions can be found in [11,14]. The mesh data structures are created and managed using the C++ mesh manager tool described in [15].
For the approximation of the velocity, we compare the polynomial quantities Π E k (u h ) and Π E k−1 (∇u h ) with the exact velocity u and the gradient ∇u. We recall that we can compute these  Relative error curves for the non-conforming VEM approximation of the velocity (left) and its gradient (right) on the three mesh families M 1 -M 3 with k from 1 to 4. The expected slopes are indicated by triangles and labels on the plots.
projections exactly using only the degrees of freedom of the vector and scalar field u h although we do not know this field. For the approximation of the pressure, we compare the piecewise polynomial fields p h and ∇p h with the exact pressure p and gradient ∇p.
The relative error curves for pressure and velocity versus the mesh size h are shown respectively in the log-log plots of Figures 6.2 and 6.3 for the three mesh sequences as indicated therein. The plots on the left show the relative errors for the approximation of the velocity or pressure field, while the plots on the right show the relative errors for the approximation of the field's gradient. The expected slopes are shown for each error curve directly on the plots and indicated by numerical labels.
These results are in very good agreement with the convergence rates that are predicted by the analysis of the previous sections.

Conclusions.
We presented the non-conforming formulation of the virtual element method for the steady Stokes problem. We have been able to construct approximations of any order in two and three space dimensions in a unified fashion, a feat which is still out of reach for standard non-conforming finite elements. Moreover, the method is naturally defined on general polygonal and polyhedral meshes. In particular non-convex polygons and polyhedra with parallel adjacent interfaces are allowed. The formulation of the method relies on the element-wise construction of virtual approximation spaces for velocity and pressure, which are characterized by a positive integer order k. The local approximation space for the velocity contains vectors of polynomials of order k plus other functions that are not computed explicitly and dealt only in terms of their degrees of freedom. The local approximation space for the pressure consists of all polynomials of order k −1. We proved that the velocity-pressure pair of global approximation spaces satisfies the inf-sup condition, from which the stability and well-posedness of the scheme follow. We also proved the optimal convergence of the numerical approximation to the velocity and pressure solution fields. The accuracy of the approximation is determined by the degree k of the polynomials and optimal a priori error estimates was derived for the velocity and pressure.