Stability Analysis for the Virtual Element Method

We analyse the Virtual Element Methods (VEM) on a simple elliptic model problem, allowing for more general meshes than the one typically considered in the VEM literature. For instance, meshes with arbitrarily small edges (with respect to the parent element diameter), can be dealt with. Our general approach applies to different choices of the stability form, including, for example, the"classical"one introduced in [L. Beirao da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (2013), no. 1, 199-214], and a recent one presented in [Wriggers, P., Rust, W.T., and Reddy, B.D., A virtual element method for contact, submitted for publication]. Finally, we show that the stabilization term can be simplified by dropping the contribution of the internal-to-the-element degrees of freedom. The resulting stabilization form, involving only the boundary degrees of freedom, can be used in the VEM scheme without affecting the stability and convergence properties. The numerical tests are in accordance with the theoretical predictions.


Introduction
The virtual element method (VEM) has been introduced recently in [3, 4, 11, 1] as a generalization of the finite element method that allows to make use of general polygonal/polyhedral meshes. The virtual element method, that enjoyed an increasing interest in the recent literature, has been developed in many aspects and applied to many different problems; we here cite only a few works [14,5,6,7,8,10,19,12,23,24,29,28] in addition to the ones above, without pretending to be exhaustive. We also note that VEM is not the only recent method that can make use of polytopal meshes: we refer, again as a minimal sample list of papers, to [13,15,16,20,25,26,27].
A VEM scheme may be seen as a Galerkin method built by means of two parts: 1. a first term strongly consistent on polynomials, which guarantees the accuracy; 2. a stabilization term, involving a suitably designed bilinear form, typically written as the sum of two contributions: s E (·, ·) = s • E (·, ·)+s ∂ E (·, ·). The form s • E (·, ·) uses the interior degrees of freedom, while the form s ∂ E (·, ·) uses the boundary degrees of freedom.
We remark that under the usual assumptions on the polygonal mesh (namely, shape regularity and the property that the length of each edge is uniformly comparable to the diameter of the parent element), devising and proving the stability features of the form s E (·, ·) is quite simple. This is the reason why, in the VEM literature, the focus is on describing explicit expression for s E (·, ·), while the proof of the corresponding stability result is often omitted. Instead, the stability analysis is more involved if one allows for more general mesh assumptions (for instance dropping the edge length condition mentioned above).
The present paper focuses on the stability properties of the bilinear form s E (·, ·). Although the approach we follow is quite general, we here consider the problem and notation of [3,5] in order to keep the presentation clearer. Our main results are the following.
• The development of a new strategy to prove the convergence of the VEM schemes, which requires weaker stability conditions on s E (·, ·) than the usual ones. Our approach is used to analyse the situations described below.
1. VEM schemes using a sequence of meshes with minor restrictions than the ones usually requested. In particular, our analysis covers some instance of shape regular meshes with edges arbitrarily short with respect to the diameters of the elements they belong to.
2. Different instances of stabilization forms s E (·, ·). Among them, we provide a detailed analysis of both the standard choice presented in [3,5], and a new one proposed in [29]. In addition, it is worth remarking that a stability analysis for this latter choice could be developed using the tools of [3]. However, the resulting error bound would be sub-optimal, in contrast with the numerical evidences. Our new approach, instead, leads to establish error bounds in perfect accordance with the numerical tests. We also show that the choice presented in [29], can have superior robustness properties in the presence of "small" edges.
• The development of a stability result concerning the choice of s E (·, ·) presented in [3,5] that is valid under more general mesh assumptions. Essentially, we prove that the stabilization term is equivalent to the H 1 seminorm, where one of the two equivalence constants logarithmically degenerates in presence of "small" edges.
• An interesting result regarding the structure of s E (·, ·). More precisely, we prove that the internal term s • E (·, ·) can be dropped without any detriment to the stability features of the underlying VEM scheme.
A brief outline of the paper is as follows. We present the continuous model problem and we review its virtual element discretization in Section 2. In Section 3 we develop a set of basic technical lemmas concerning virtual elements and polygonal elements. In Section 4 we present and develop our error analysis strategy. Afterwards, in Section 5 we apply such an approach in order to analyse some existing choices of the stability form, under more general mesh assumptions than the ones typically adopted in the VEM literature. Finally, we present some numerical tests in Section 6.

The continuous and discrete problems
In this section we briefly present the continuous problem and its discretization with the Virtual Element Method. More details can be found in [3,5].

The continuous problem
As a model elliptic problem we consider the diffusion problem in primal form. Defining (·, ·) as the scalar product in L 2 , and a(u, v) := (K∇u, ∇v), the variational formulation of the problem reads: Find where Ω ⊂ R 2 is a polygonal domain and the loading f ∈ L 2 (Ω). The diffusion symmetric tensor K = K(x, y) is assumed to satisfy: Above, | · | denotes the euclidean norm in R 2 . It is well known that problem (1) has a unique solution, because our assumptions on K and the Poincaré inequality yield: with 0 < α < M < 1. Note that the bilinear form a(·, ·) in (2) can obviously be split as for all v, w ∈ V .

The virtual element method
Let an integer k, equal or greater than 1, and let {Ω h } h denote a family of meshes, made of general simple polygons, on Ω. Given an element E ∈ Ω h of diameter h E and area |E|, its boundary ∂E is subdivided into N = N (E) straight segments, which are called edges, with a little abuse of terminology. Accordingly, the endpoints of the edges are called vertices of the element E. We remark that several consecutive edges of E may be collinear; as a consequence, the number of edges (and vertices) may be greater than the number of maximal straight segments of ∂E. Hence, a triangle may have ten edges, for instance. Furthermore, the length of an edge e ∈ ∂E is denoted by h e . Moreover, in the sequel we assume that the diffusion tensor K is piecewise constant with respect to the meshes {Ω h } h .
For each E ∈ Ω h we now introduce the local virtual space where P n , n ∈ N, denotes the polynomial space of degree n, n ∈ N, with the convention that P −1 = {0}. The associated set of local degrees of freedom Ξ (divided into boundary ones Ξ ∂ , and internal ones Ξ • ) are given by • point values at the vertexes of E; • for each edge, point values at (k −1) distinct points on the edge (this are typically taken as Gauss-Lobatto nodes, see [3,5]); • the internal moments against a scaled polynomial basis For future reference, we collect all the N k boundary degrees of freedom (the first two items above) and denote them with is obtained by gluing the above spaces, and the same holds for the global degrees of freedom. We refer to [3] for the explicit expression. On each element E we also define a projector Π ∇ E : V E → P k (E), orthogonal with respect to the bilinear form a E (·, ·). More explicitly, for all v ∈ V E : where ∂ s denotes the tangent derivative along the edges.
Another interesting point considered in this paper, is the possibility to completely neglect the internal part of the stability form. In other words, we will show that the choice does not spoil the stability feature of the numerical scheme.
Given any symmetric and coercive form s E (·, ·), one can define the local discrete bilinear forms on that are computable and approximate a E (·, ·). Given the global discrete form the discrete problem is: For a discussion about the approximated loading term < f h , v h >, we refer to [3,1]. In order to shorten the notation, and also to underline the generality of the proposed approach, in the following we will simply use Π E instead of Π ∇ E to denote the projector operator.
The following assumptions on the mesh will be considered in the present work. In all the presented results we will explicitly write which of those hypotheses are used, if any. We remark that assumptions A1 and A3 are those considered in [3]. Here, we want to consider also weaker assumptions in terms of the edge requirements, namely the combination of A1 and A2. It is easy to check that, provided A1 holds, assumption A3 implies A2. However, assumption A2 is much weaker than A3, as it allows for edges arbitrarily small with respect to the element diameter.
In the following the symbol will denote a bound up to a constant that is uniform for all E ∈ {Ω h } h (but may depend on the polynomial degree k). Moreover, in order to make the notation shorter, for any non negative real s we will denote by the standard H s Sobolev (semi)norm on the measurable open set ω.

Remark 2.
The stability form s E (·, ·) may also be scaled by a multiplicative factor τ E > 0, to take into account the magnitude of the material parameter K, for instance. With this respect, a possible choice could be to set τ E as the trace of K on each element. In this paper, we do not investigate on how to select τ E , but we address the reader to [4,19,14] for some study on such an issue.

Preliminary results
In this section we present some technical results that will be needed in the sequel of the paper. The H 1/2 boundary norm will have an important role in the following. We here use the (one-dimensional) double integral definition where, with a small abuse of notation, v stands for v| ∂E , and where s 1 , s 2 denote curvilinear abscissae along the boundary. We begin with the following Lemmas.
Moreover, for all E ∈ Ω h and all v ∈ H 1/2 (∂E), there exists an Proof. We only sketch the simple proof, based on a mapping argument. Up to a translation of the element E, we may assume that the ball center x E is the origin of the coordinate axes. Let then the function Ψ : [0, 2π) → [ρ E , h E ] describe the boundary of E, as follows. The boundary curve Γ = ∂E can be parametrized in a unique way as with θ representing the angle in radial coordinates. Note that property A1 implies Ψ ∈ W 1,∞ [0, 2π), uniformly with respect to E ∈ Ω h . We then introduce the radial mapping F : B E → E, associating a point expressed in polar coordinates x,ŷ = r cos (θ),r sin (θ) ,r ∈ [0, ρ E ],θ ∈ [0, 2π), with the point (x, y) = F (x,ŷ), whose coordinates are x, y = r cos (θ), r sin (θ) , r =r Ψ(θ) ρ E , θ =θ.
By recalling A1, it can be checked that F ∈ W 1,∞ (B E ), and the same holds for the inverse mapping, i.e. F −1 ∈ W 1,∞ (E). It is easy to see that |F | 1,∞,B E C and |F −1 | 1,∞,E C. As a consequence, bound (17) can be simply proved by a standard "pull-back and push-forward" argument: i) map v ∈ H 1 (E) from E into B E using F ; ii) notice that the trace bound analogous to (17) holds on the ball B E ; iii) map back to E using F −1 . Bound (18) is similarly proved: one only needs to map the boundary data into B E , to consider the harmonic extension inside B E , and finally to map back to E.

Corollary 3.2. Let assumption A1 hold. Then
with n E denoting the outward unit normal to the boundary of E.
Proof. By the definition of dual norm and using (18), we get An integration by parts, using div w = 0, and the Cauchy-Schwarz inequality lead to estimate (20):

Lemma 3.3.
Let assumption A1 hold true. Then we have: Proof. For E ∈ Ω h , let T E ⊂ E denote an equilateral triangle inscribed in the ball B E . We start observing that, due to assumption A1, for any polynomial p of given maximum degree it holds p 0,E p 0,T E . This follows from noting that the smallest ball containing E and the largest ball contained in T E have uniformly comparable radii. We now recall that ∆v h ∈ P k−2 . Let b ∈ P 3 (T E ) denote the standard cubic bubble in T E with unitary maximum value. Standard properties and inverse estimates of polynomial spaces on shape regular triangles yield

Remark 3. The same argument in the proof of Lemma 3.3 can be used to prove inverse estimates for polynomials of fixed maximum degree, on polygons satisfying assumption A1.
The next Lemma can be considered as a variant of Lemma 3.1 in [9], supposing that the number of edges is uniformly bounded.

Lemma 3.4. Let A1 and A2 hold. For all E ∈ Ω h and all
Proof. For w h ∈ V E , we recall that w h|∂E ∈ C 0 (∂E) and w h|∂E is a polynomial of degree at most k on each edge. In addition, we first suppose that ∂E w h = 0. Then, by definition (16) and by using a scaling argument on each edge of the mesh, we obtain where s denotes the curvilinear abscissae along the generic edge. By assumption A2 and recalling that w h is continuous on the boundary, from the above bound we have: where we also used that w h | ∂E has zero average and thus it vanishes at least at one point of ∂E. For a generic v h ∈ V E (not necessarily with vanishing mean value), the proof follows easily from (22) by adding and subtracting its average on the boundary v h and simple bounds: where |∂E| denotes the length of ∂E.
The following approximation result is an extension of the one in [23] to the case of higher order norms and more general mesh assumptions. Theorem 3.5. Let assumption A1 hold. Then there exists a real number σ > 3/2 such that for all u ∈ H s (Ω), 1 < s ≤ k + 1, and all 1 ≤ σ < min{σ, s}, it holds: Proof. For each element E, we can build a sub-triangulation by connecting all its vertexes with the center x E introduced in assumption A1. We denote by T h the global (conforming) triangular mesh obtained by applying such a procedure for all E ∈ Ω h . It is easy to check that, under assumption A1, the triangles in the sequence of meshes {T h } h have maximum angles that are uniformly bounded away from π (although shape regularity is not guaranteed). Let u r be the standard continuous and piecewise P k polynomial Lagrange interpolant of u over the triangulation T h . Then it holds: where we used the anisotropic approximation results in [2], also recalling the angle property above. In the following we denote by u π a piecewise discontinuous polynomial approximation of u over the mesh Ω h . For instance, one may think of the L 2 projection of u on P k (E) for each element E. We now introduce the function u I ∈ V h defined, on each element E, by Therefore, for all E ∈ Ω h , regularity results on Lipschitz domains (see [21]) guarantee that where σ E = 2 if E is convex, and σ E = 1 + π/ω E (with ω E the largest angle of E) otherwise. Let now σ = min E∈{Ω h } h σ E , where we stress that the minimum is taken among all elements of the whole mesh sequence. Due to assumption A1, that yields a uniform bound on the maximum element angles, the number σ is strictly bigger than 3/2. First a triangle inequality and bound (26), then a trace inequality yield for all 1 ≤ σ ≤ σ and all E ∈ Ω h . The result follows combining the above bound with (24) and standard polynomial approximation estimates on shape regular polygons.
Regarding the operator R, we have the following Lemma.
1. For R defined by (6) or by (7): unless v is a polynomial, in which case it holds: Finally, if also assumption A2 holds, then Proof. Estimates (27) and (28) follow, recalling assumption A1, from standard approximation theory on shape regular polygons. Bound (29) follows immediately from (28) by using an inverse estimate for polynomials on polygons, see Remark 3.
To prove (30), take any v h ∈ V E and set v h : Furthermore, from (8) and recalling that Rv h is a constant, we get Since v h − v h has zero mean value on ∂E, using Lemma 3.4 we get Hence, from (33) and (17), we obtain Estimates (31), (32) and (34) give (30).

A general error analysis
In the present section we derive an error analysis which is more general than the standard one detailed in [3]. We remark that the present approach can be applied to any other linear symmetric elliptic problem. For the analysis, the following discrete semi-norm, induced by the stability term, will play an important role: Above, V E ⊆ V |E is a subspace of sufficiently regular functions in order for s E (·, ·) to make sense. We now introduce the following assumption, for all E ∈ T h .
Main assumption -We assume that it holds with C 1 (E), C 2 (E) positive constants which depend on the shape and possibly on the size of E.
and also the bound where Proof. We start by noting that, from definition (5) it is immediate to check that Using first (40), then noting that Π E (I − Π E ) = 0 and applying (37), we obtain (cf. Again using the first identity in (41), recalling definition (35), from the triangle inequality we get Since Π E is a projection with respect to a E and using (37) we obtain From (43) we immediately get and also, recalling (36), Combining the above bounds it follows As an immediate consequence of Lemma 4.1 and (36), the discrete bilinear form (14) associated to (13) satisfies where Therefore, due to (2), the discrete problem is positive definite and problem (15) has a unique solution.
We have moreover the following convergence result. For all sufficiently regular functions v, introduce the global semi-norms We notice that, by (36) and the Poincaré inequality, ||| · ||| is a norm on V h , not only a semi-norm. Furthermore, for any h, let F h denote the quantity We remark that, again by (36), it holds: Therefore, taking f h as in [3] and using the arguments in that paper, we infer: where C(f ) depends on suitable Sobolev norms of the source term f . We have the following result.

Theorem 4.2.
Let assumptions (36)-(37) hold and let the continuous solution of (1) satisfy u |E ∈ V E for all E ∈ T h . Then, for every u I ∈ V h and for every u π such that u π|E ∈ P k (E), the discrete solution u h of (15) with bilinear form (13) satisfies Proof. First using the coercivity property in Lemma 4.1, then with identical calculations as in Theorem 3.1, equation (3.11), of [3], we get where , C 2 (E)}, and the terms T i are given by For term T 1 , definition (47) and assumption (36) yield where Term T 2 is treated using both the bounds (38) and (39), that easily lead to the estimate where C (h) = max From (52), using the bounds (53), (54) and (55), then dividing by |||u h − u I |||, we get The triangle inequality and (36) give Combining (56) and (57), we get (50) with Remark 4. By using (56) and the triangle inequality it is immediate to check that, as a corollary of the above result, it also holds

Reduction to the boundary
In the present section we derive a result that allows to focus the analysis of assumptions (36) and (37) only on the boundary of the element. We here consider two cases for the internal stabilization form s • E (·, ·): 1. as in the standard VEM (e.g. [3]), we put (see (10)) 2. we completely neglect the internal contribution (see (12)), i.e.
Both cases will lead to the same results. Instead, the boundary bilinear form s ∂ E (·, ·) is left completely general for the moment, and different choices will be made and analysed in Section 5.
We start by showing the following Lemma.

Lemma 4.3. For all
Proof. We only sketch the very simple proof. Since for all v h ∈ V E it holds ∆v h ∈ P k−2 (E), there are (infinitely many) polynomials of degree k that satisfy ∆ p = ∆v h (cf. [?], for instance). In order to derive the bound (58), we first note that, thanks to assumption A1 and since ∆ p = ∆v h , inequality (58) is equivalent to The above bound, that is now restricted on balls, can be easily deduced by choosing p in the subspace and by a scaling argument.
Concerning assumption (36), we have the following result.
(60) Again an integration by parts and (60) yield with n E denoting the outward unit normal to the boundary of E. The first term above is bounded by the Cauchy-Schwarz inequality, Lemma 4.3 and Lemma 3.3. We obtain For the second term, we first note that div(∇(v h − p)) = ∆(v h − p) = 0. Therefore, after applying a (scaled) duality bound on the boundary of E, we can use Corollary 3.2 with w = ∇(v h − p), and obtain Moreover, by standard approximation estimates in one dimension, it holds h Therefore, using (59), the triangle inequality and again Lemmas 4.3 and 3.3, bound (63) yields The result follows by combining equations (61), (62), (64) and recalling that |v h | 2 Furthermore, concerning assumption (37), we have the following result.
Proof. We first note that the second term in (35) is immediately bounded: Therefore, by the definition of s E (·, ·) and using (65), it is sufficient to show that Clearly, the above bound is trivial for the choice (12). Hence, we can focus on the choice (10). By definition of s • E (·, ·) and recalling that m i L ∞ (E) 1, i = 1, 2, ..., n k−2 , we have Using property (27) for the operator R (or using (29) if the choice (8) is being used), we now have

Analysis of some choices for the boundary stabilization
In the present section we apply Propositions 4.4 and 4.5 for a couple of standard choices of the boundary stability term s ∂ E (·, ·). This allows to relax the mesh assumptions (with respect to the theory presented in [3]) in establishing stability and convergence properties of the proposed methods.

Identity matrix choice
This is the more standard, and simpler to code, choice for virtual elements. We recall it here again, for convenience: We may call it the identity matrix choice since in the implementation procedure of the method, the bilinear form (70) is clearly associated with an identity matrix of dimension N k.
Proof. We first recall that ∂E is meshed by means of its edges, so that ∂E = ∪ N j=1 e j . We also define h j := |e j |. Moreover, in the proof we will make use of the space H 1/2 00 (Γ), where Γ is a connected part of ∂E with |Γ| > 0. This space is defined by, see [22]: where Ext(v) denotes the extension by zero of v to the whole ∂E. Its norm where ρ(x) denotes the distance of x from ∂Γ, is equivalent to |Ext(v)| 1/2,∂E .
Given v h ∈ V E , we set v L ∈ V E as the usual piecewise linear Lagrange interpolant of v h , relative to the edge mesh. We have We now define w j = χ e j (v h − v L ) and we notice that, since v h − v L vanishes at all the nodes, we have Exploiting that w j is a polynomial of degrees ≤ k on e j , a scaling argument shows that Therefore, recalling assumption A2 and using that ||v by which It remains to estimate |v L | 2 1/2,∂E . We denote by ϕ i the usual hat function with support σ i := e i−1 ∪ e i (here i − 1 and i are intended modulo N). We write where v i ∈ R is the value of v L at the i-th node. We have, using assumption A2: Recalling (73), direct computations show that by which we obtain Therefore, using again assumption A2 and noting that from (78) and (79) we get Combining (74), (77) and (81), we get (71).
We now have the following stability result.
Proof. Standard results for polynomials in one dimension immediately give A combination of (83) and Lemma 5.1 yields: i.e. condition (59) holds with C 1 (E) (log (1 + h E /h m(E) )). We now prove that estimate (65) holds. Recalling assumption A2, it is immediate to check that Take any p ∈ P k (E). We get, using bound (85) , an inverse estimate for polynomials (cf. Remark 3), and recalling either (27) or (29) (depending on the choice of the operator R): i.e. condition (65) holds with C 2 (E) 1.
The following corollary shows that, even in the presence of arbitrarily small edges (provided the number of edges are uniformly bounded), the convergence rate of the Virtual Element Method is quasi-optimal, in the sense that only a logarithmic factor is lost. (1), assumed to be in H s (Ω), s > 1. Let u h be the solution of the discrete problem (15).

Corollary 5.3. Let assumptions A1 and A2 hold. Let u be the solution of problem
If the stronger assumption A3 holds, then clearly c(h) 1.

A "classical" stability bound
We close this part on the identity matrix choice by showing that the classical stability result of [3], see equation (3.7) of [3], can also be proved under the more general mesh assumptions considered in this paper. This result could be used to prove the same error estimate as in Corollary 5.3 without resorting to the approach described in this paper, but simply applying the standard theory of [3]. We need an additional preliminary lemma.
Lemma 5.4. Let assumption A1 hold. We have and for all v in H 1 (E).
Proof. The simple proof is based on an anisotropic scaling argument. Take an edge e ∈ ∂E, and let T ∈ T h be the associated triangle (see the proof of Theorem 3.5). By a rotation and translation of the cartesian (x, y)-coordinates, it is not restrictive to assume that e = {0} × [−h e /2, h e /2], and that the center of the ball, see assumption A1, x E = (x E , y E ) satisfies x E ≥ 0. As a consequence of assumption A1, it is easy to check that x E h E , y E h E and that the ball B E is contained in the half plane {(x, y) ∈ R 2 : x ≥ 0}. Therefore, we also have x E h E . Let nowT be the triangle of vertexes (0, h e /2), (0, −h e /2), (h e /2, 0). We now consider the unique affine mapping F : T →T that leaves the edge e (and its orientation) unchanged:ê := F (e) = e. By an explicit computation of F and its inverse F −1 , we get the Jacobian matrices The proof now follows by a scaling argument. Indeed, denotingv = v •F −1 , well known (scaled) trace estimates onT and a simple change of variables give By recalling the upper and lower bounds on (x E , y E ), the above estimate yields that immediately implies (102) by summing over all e ∈ ∂E.
Proposition 5.5. Let assumptions A1 and A2 hold. Then, for any of the choices (6), (7), (8) of the operator R, it holds where Note that if the stronger assumption A3 holds, then clearly c(h) 1.
We now show the other bound. We first consider s • E (v h , v h ) and notice that there is nothing to estimate if the choice (12) has been done. Therefore, we focus on the choice (10). By definition, the Holder inequality and recalling m i L ∞ 1, we get the estimate: that, using Lemmas 3.4, 3.1 and 5.4 yields Since Πv h = 0 implies Rv h = 0, either bound (27) (for the choices (6) and (7)), or bound (30) (for the choice (8)), yields The first bound in (103) now follows from combining (104), (105) and (106) and noting that |v h | 2

A stabilization based on boundary derivatives
We now analyze a different choice for the boundary part of the stabilization term, namely the one given by (cf. [29]): We highlight that, contrary to the identity matrix stabilization presented in section 5.1, the standard approach of [3] applied to (107), would lead to a strongly suboptimal result in the presence of small edges. Indeed, the term (107) Proof. We first prove that condition (65) is fulfilled. Take any p ∈ P k (E). Using assumption A1 and an inverse inequality for polynomials (cf. Remark 3), we get i.e. condition (65) holds with C 2 (E) 1.
To prove that condition (59) is fulfilled, we simply notice that (110) and we obtain that condition (59) holds with C 1 (E) 1.
Corollary 5.7. Let assumption A1 hold. Let u be the solution of problem (1), assumed to be in H s (Ω), s > 3/2. Let u h be the solution of the discrete problem (15), with the choice (107). Then it holds Proof. Theorem 5.6 allows to apply Propositions 4.4 and 4.5. Therefore, assumptions (36) and (37) hold with C 1 (E) 1 and C 2 (E) 1, respectively. Then, Theorem 4.2 can be invoked with C err (h) satisfying C err (h) 1. We now estimate the terms in the right-hand side of (50). Using exactly the same arguments of Corollary 5.3, and focusing again only on the non-trivial choice (10), we get: and (114) Therefore, we only need to estimate the boundary part To this end, take s > 3/2 and σ such that 3/2 < σ < s. We have, using a scaled trace inequality (use a similar argument to that in Lemma 3.1 for the function ∇(u − u I )), Hence, Theorem 3.5 gives: Combining (114) and (116), we get Similarly, using also standard approximation results on polygons (see [18]), we get We conclude by collecting estimates (111), (112), (113), (117) and (118).

Remark 5.
The same analysis can be employed to prove error estimates for many other choices of the stabilization. We here spend some word on the following variants of choice (107).
The first variant is an "L 2 -version" of (107): Since it is easy to check that, under assumptions A1 and A2, it holds

Numerical tests
For all numerical tests we will consider Laplace equation on the unit square Ω :=]0, 1[ 2 : with right hand side f and Dirichlet boundary condition g defined in such a way that the exact solution is u ex (x, y) := x 3 −xy 2 +x 2 y +x 2 −xy −x+y −1+sin(5x) sin(7y)+log(1+x 2 +y 4 ). (121) In the following brief experiments we will address some of the issues considered in the paper.

Small edges
In the first numerical experiment we consider the issue of the presence of very small edges. On the one hand, we show that the classical VEM stabilization (70) can generate small oscillations, that are of the order of the approximation error. In the case k = 1 these oscillations are visible and, depending on the application, may be preferable to avoid. However, already for k = 2 the oscillations become so small to be practically negligible. On the other hand, we show that the stabilization (107) eliminates this oscillations already for the k = 1 case. We consider a mesh obtained by gluing together two distinct meshes along x = 0.5; this case can happen for instance in contact problems, see [29]. The mesh is shown in Fig. 2, while in Fig. 3 we show the section of the mesh at x = 0.5. Note that around y = 0.4 there is a very small edge of length 3.21 × 10 −4 .
In Figs. 4 and 5 we plot the section at x = 0.5 of the VEM solution for the classical stabilization (70) (thick line) together with the exact solution (thinner line) for k = 1 and k = 2 respectively. A careful inspection shows that that in Fig. 4 there are some small oscillations in correspondence of the small edges. In Fig. 5 the oscillations are no more visible but are still present. We reproduce the same experiments in Figs. 6 and 7 with the boundary stabilization (107). Now the oscillation have disappeared also for k = 1 but in this case the VEM solution seems to be less accurate. The motivation is that the boundary stabilization (107) is too strong. The situation can be improved by taking a smaller stabilization parameter (see Remark 2); in Figs. 8 and 9 we show the same experiments with τ E = τ = 0.1. We have developed several further experiments (here not shown) using different meshes and loading, and choosing τ = 0.1 for s E (·, ·) as in (107): the obtained results were always accurate. This is a general property of VEM: the sensitivity of the method with respect to the stabilization parameter is very mild when considering different meshes and loading/boundary data. Nevertheless, a detailed study on such an issue is beyond the scope of the present paper.

Convergence in H 1
We will show, in a loglog scale, the convergence curves of the error in the H 1 seminorm between the exact solution u ex and the solution u h given by the Virtual Element Method. As the VEM solution u h is not explicitly known inside the elements, we compare ∇u ex with the elementwise L 2 −projection of ∇u h onto P k−1 , that is, with Π 0 k−1 ∇u h . It is easy to see that this latter quantity can indeed be computed starting from the degrees of freedom of u h . For the convergence test we consider four sequences of meshes.
The first sequence of meshes (labelled square) is simply a decomposition of the  Fig. 10 and in Fig. 11 respectively. The third sequence of meshes (labelled Lloyd-0) is a random Voronoi polygonal tessellation of the unit square in 25, 100, 400 and 1600 polygons. The fourth sequence (labelled Lloyd-100) is obtained starting from the previous one and performing 100 Lloyd iterations leading to a Centroidal Voronoi Tessellation (CVT) (see e.g. [17]). The 100-polygon mesh of each family is shown in Fig. 12 (Lloyd-0) and in Fig. 13 (Lloyd-100) respectively.
In the figures from 14 to 17 we plot for k = 1 (low order case) the H 1 error on each mesh family as a function of the mean diameter h of the polygons. We consider the classical stabilization (70) (solid line), the boundary stabilization (107) with τ = 1 (dotted line), and the boundary stabilization (107) with τ = 0.1 (dashed line). In the figures from 18 to 21 we finally plot the same values for k = 5 (as a sample high order case). We observe that as h goes to zero all stabilizations behave very similarly, namely as O(h k ), as predicted by the theory.