Interior estimates for the Virtual Element Method

We analyze the local accuracy of the virtual element method. More precisely, we prove an error bound similar to the one holding for the finite element method, namely, that the local $H^1$ error in a interior subdomain is bounded by a term behaving like the best approximation allowed by the local smoothness of the solution in a larger interior subdomain plus the global error measured in a negative norm.


Introduction
Besides its ability to handle complex geometries, one of the features that contributed to the success of the finite element method as a tool for solving second order elliptic equations is their local behavior.Considering, to fix the ideas, the Poisson equation in Ω u = 0, on ∂Ω, the standard, well known, error estimate provides, for the order k finite element method, an error bound of the form (u and u h denoting, respectively, the exact and approximate solutions, and h the mesh size).
When measured in a global norm, the error can be negatively affected by the presence of even a few isolated singularities.However, since the early days in the history of such a method it is well known that, if a solution with low overall regularity is locally smoother, say in a subdomain Ω 1 ⊂⊂ Ω, an asymptotically higher order of convergence can be expected in any domain Ω 0 ⊂⊂ Ω 1 .More precisely, there exists an h 0 (depending on Ω 0 and Ω 1 ) such that, for h < h 0 , the H 1 (Ω 0 ) norm of the error can be bounded ( [37], see also [39,40]) as which, combined with an Aubin-Nitsche argument to bound the negative norm on the right hand side, yields, if Ω is sufficiently smooth, an O(h k ) bound for the error in the H 1 (Ω 0 ) norm, provided u ∈ H k+1 (Ω 1 ), even when u ∈ H k+1 (Ω).This feature is particularly appealing, as it allows to take advantage of the local regularity of the solution, thus enabling the method to perform effectively.One can, for instance, avoid the need of refining the mesh, whenever possible singularities are localized far from a region of interest.
It is therefore clearly desirable that new methods, aimed at generalizing finite elements, retain this property.We focus here on the virtual element method, a discretization approach that generalizes finite elements to general polygonal space tessellations.Analogously to the finite element method, the virtual element discretization space is continuously assembled from local spaces, constructed element by element in such a way that polynomials up to order k are included in the local space.
Contrary to the finite element case, however, the functions in the space are not known in closed form but are themselves solution to a partial differential equation, which is, however, never solved in the implementation.In order to handle the discrete functions, these are instead split as the sum of an exactly computable polynomial part, and of a non polynomial part.Exact handling of the polynomial part alone turns out to be sufficient to guarantee good approximation properties: to this end, the bricks needed for the solution of the problem at hand by a Galerkin approach (e.g.local contribution to the bilinear form and right hand side) are computed as a function of a set of unisolvent degrees of freedom, in a way that is locally exact for polynomials.The non polynomial part is instead handled by means of a stabilization term, which only needs to be spectrally equivalent to the bilinear form considered, resulting in a non conforming approximation.Since its introduction in the early 2010s (see [5,7]), the virtual element method has gained the interest of the scientific community and has seen a rapid development, with numerous contributions aimed at the theoretical analysis of the method (see e.g.[14,5,24,23] ), its efficient implementation (see, e.g., [20,22,32,31,25]), its extensions in different directions (see, e.g., [21,9,11,10,12,8,27]), and applications in different fields, such as fluid dynamics [16,15,3,18], continuum mechanics [30,28,13,6,34,42,43], electromagnetism [29] and others ([4, 38, 2]).
In this paper, we aim at proving that the approximation by the virtual element method has good localization properties, similar to the ones displayed by the finite element method.More precisely, under suitable assumptions on the tessellation, we will prove that an estimate of the form (1.1) also holds for the virtual element solution (see Theorem 6.1).Under suitable assumptions on the domain Ω (the same needed for the analogous result in the finite element method), this will imply that, provided u ∈ H k+1 (Ω 1 ), we have that u − u h 1,Ω 0 = O(h k ), independently of the overall smoothness of the solution.
The paper is organized as follows.After presenting some notation and recalling how some inequalities do (or do not) depend on the shape and size of the elements (see Section 2), in Section 3 we present the virtual element formulation we will be focusing on, and we will study the equation satisfied by the error.In Section 4 we will study how different linear and bilinear operators commute with the multiplication by a smooth weighting function.In Sections 5.1, 5.2 and 5.3 we will provide bounds for the error in, global and local negative norms.In Section 6 we will prove the main result, namely Theorem 6.1, and leverage it to obtain local error bounds (see Corollaries 6.4 and 6.5).In Section 6.1 we briefly sketch an extension of the local error bounds to the so called enhanced version of virtual element method [1], which is often the one that can be found in actual implementations.Finally, in Section 7, we present some numerical results.
Throughout the paper, we will write A B to indicate that A ≤ cB, with c independent of the mesh size parameters, and depending on the shape of the elements only through the constants γ 0 and γ 1 in the shape regularity Assumption 2.1.The notation A B will stand for A B A.

Notation and preliminary bounds
In the following we will use the standard notation for Sobolev spaces of both positive and negative index, and for the respective norms (see [36]).Letting Ω ⊂ R 2 denote a bounded polygonal domain, we will consider a family F = {T h } of polygonal tessellations of Ω, depending on a mesh size parameter h.We make the following assumption on the tessellations.Assumption 2.1.There exist constants γ 0 , γ 1 > 0 such that, letting h K denote the diameter of the polygon K, for all tessellation T h ∈ F: (a) All polygons K ∈ T h are star shaped with respect to all points of a ball with center x K and radius ρ K with ρ K ≥ γ 0 h K ; (b) for all K ∈ T h the distance between any two vertices of K is greater than γ 1 h K .
Moreover, for the sake of notational simplicity, we assume that all tessellations are quasi-uniform, that is for all K ∈ T h , we have that h K h.
Trace inequalities.Under Assumption 2.1(a) for all v ∈ H 1 (K), we have We have the following proposition, where the norm • −1,K is defined as Proposition 2.2.Under Assumption 2.1(a), for all u ∈ H 1 (K), we have Proof.We split v as v H + v 0 with ∆v H = 0 and v 0 ∈ H 1 0 (K).The splitting is stable with respect to the H 1 semi norm, that is we have Now, for w ∈ H 1 0 (K) arbitrary, since v H = v on ∂K, integrating by parts twice we can write where n K denotes the outer unit normal to ∂K.Thanks to the arbitrariness of w, dividing both sides by |v H | 1,K and using (2.3) we obtain On the other hand, as ∆v 0 = ∆v, we can write By collecting (2.6) and (2.7) into (2.5)we obtain the desired bound.
For D being a polygonal element K or an edge e of T h , we let P (D) denote the restriction to D of the space of bivariate polynomials of order up to .
Under Assumption 2.1, we have the following polynomial approximation bounds [33]: and, letting h e denote the length of the edge e, for all v ∈ H t (e), 0 ≤ s ≤ t ≤ + 1 (2.9) inf q∈P (e) v − q s,e h t−s e |v| t,e .
Moreover, the following inverse inequalities for polynomial functions hold: for all q ∈ P (K) and q ∈ P (e), and all 0 ≤ s ≤ t (2.10) q t,K h −(t−s) K q s,K , q t,e h −(t−s) e q s,e .

The Virtual Element Method
In order to introduce the notation, let us review the definition of the simple form of the virtual element method that we are going to consider.Letting Ω ⊂ R 2 , we focus on the following model problem: (3.1) which, in weak form, rewrites as: find u ∈ H 1 0 (Ω) such that Let T h ∈ F denote a tessellation of Ω in the family F. For reasons that will be clear in the following, we consider a form of the Virtual Element discretization where we allow different approximation orders on the boundary and in the interior of the elements [17].As usual, for all K ∈ T h we let B k (∂K) = {v ∈ C 0 (∂K) : v| e ∈ P k (e) for all edges e of K}.
The local element space V k m (K) is defined (see [17]) as , and ∆v ∈ P m (K)}.We assume that max{0, k − 2} ≤ m ≤ k.The case m = k − 2 corresponds to the simplest form of the virtual element method, as introduced in [5].

The following inverse inequality holds for all
Remark 3.1.For the sake of simplicity, we do not explicitly include in our analysis the simplest lowest order VEM space and ∆v h = 0}, which can be tackled by the same kind of argument but which would require a separate treatment, in particular when dealing with the terms involving the approximation of the right hand side.
The global discretization space V h is defined as for all K ∈ T h } denote the space of discontinuous piecewise H 1 functions on the tessellation T h , which we endow with the seminorm and norm we also introduce the discontinuous global discretization space , for all K ∈ T h } denote the space of discontinuous piecewise polynomials of order up to k defined on the tessellation T h , we let Π ∇ : H 1 (T h ) → P k (T h ) be defined as The discrete bilinear form a h : where I is the identity of H 1 (T h ) and where, for shortness, here and in the following we use the conventional notation ˆTh The stabilization bilinear form s h : V h × V h is defined as the sum of local contributions where we assume, as usual, that, for all K ∈ T h , the local stabilization bilinear form s K : In particular, for all v h ∈ V k m (K), (3.5) and (3.6) yield (3.7) denote the local counterpart of the bilinear form a h .We recall that, thanks to (3.5) and (3.6), for all v ∈ V h we have that We next let Π 0 m : L 2 (Ω) → P m (T h ) denote the L 2 (Ω) orthogonal projection onto the space P m (T h ) of discontinuous piecewise polynomials of order at most m, and we let The virtual element solution to Problem (3.1) is obtained by solving the following discrete problem: find 3.1.Extension of the discrete operators to H 1 (T h ).To carry out the forthcoming analysis, it will be convenient to extend some of the above operators, which are defined on the discrete space V h , to the whole H 1 (T h ).To this aim, we introduce projectors Π ∇ K : Also the projectors Π ∇ K and Q ∇ K can be assembled, element by element, to global projectors into the space V h of discontinuous virtual element functions.More precisely we define Π ∇ :

With this notation, we can extend the bilinear form a
We remark that, thanks to (3.8), we have As Q ∇ K q = 0 for all q ∈ P k (K), we easily have the following proposition, where the H 1 (K) seminorm bound stems from the best polynomial approximation and the L 2 (K) bound is obtained by a Poincaré inequality, as, by definition, Q ∇ K f is average free.
3.2.Error equations.Letting u ∈ H 1 0 (Ω) and u h ∈ V h respectively denote the solutions of Problem (3.1) and (3.10), we can now write two error equations, satisfied by u − u h .Indeed, for v ∈ H 1 (Ω) arbitrary we have where Then, letting finally yielding the following error equation Combining (3.16) and (3.14) we also have Using Proposition 3.3, we easily see that the following proposition holds.
Moreover, the following proposition provides an a priori bound for the operator δ f appearing at the right hand side of the error equation.
Remark that using the above bounds, in combination with the error equation and an approximation estimate (see e.g.(4.1)), allows to retrieve the following (essentially well known, see [5,17]) bound on the error u − u h : if the solution u and the source term Remark that, as f = −∆u, we have that ρ ≥ r − 1.Moreover, by construction m ≥ k − 2. Then the second term on the right hand side, deriving from the approximation of the source term, is asymptotically dominated by the first term, namely h min{r,k} u 1+r,Ω .

Commutator properties for the VEM space
The local bounds we aim at proving will involve multiplying different quantities by smooth weights.A key role will be played by the error resulting from commutating the action of such weights with different operators appearing in the definition and analysis of the VEM method.To analyze such errors, we start by introducing a local quasi-interpolation operator similar to the one proposed in [11] and defined as follows.Given v ∈ H 1 (K) with v| ∂K ∈ C 0 (∂K) and ∆v ∈ L 2 (K), we let m ∆v, where I K : C 0 (∂K) → B k (∂K) denotes the edge by edge interpolation operator with, as interpolation nodes, the nodes of the k + 1 points Gauss-Lobatto quadrature formula, and where, by abuse of notation, we let Π 0 m : L 2 (K) → P m (K) denotes the L 2 orthogonal projection.As (we can for instance take ωK = ω(x K ), x K being the barycenter of K).
We have the following lemma, which we prove by an approach similar to the one in [19].
with ρ K given by (4.2).Then, using Proposition 2.2 we can write We separately bound the two terms on the right hand side of (4.4), starting from the first one.We remark that, on e edge of K we have that 00 (e), where, we recall, H 1/2 00 (e) can be defined as the space of those functions v ∈ L 2 (e) such that setting v = v in e and v = 0 in ∂K \e it holds that v ∈ H 1/2 (∂K).We endow H 1/2 00 (e) with the norm v H 1/2 00 (e) = | v| 1/2,∂K .We recall that H 1/2 00 (e) is the interpolation space of exponent 1/2 with respect to the interpolation couple (L 2 (e), H 1 0 (e)).Then, using a standard interpolation bound (see [41]) and (4.1), we can write Now, using (2.9) and (2.10), we can write As far as the second term at the right hand side in (4.4) is concerned, we have as well as and then, by triangle inequality, We can bound the three terms by using a standard duality argument, which allows to bound the H −1 (K) norm of any average free function with h K times its L 2 norm, and we obtain ∇v h 0,K , and, using (4.2) as well as the inverse inequality (3.3), Using (4.6) and (4.7) in (4.4) yields Let now v ∈ H s+1 (K).Adding and subtracting Π ∇ K v, using (4.1) and (4.8), then adding and subtracting v and using the polynomial approximation bound (2.8), we have As for all w ∈ H 1 (K) it holds that we immediately have the following corollary.
Remark 4.3.The bound (4.8), valid for all v h ∈ V k m (K), is the so called discrete commutator property of the Virtual Element space.It implies that In particular, we have Another commutativity bound that will play a role later on, is the following lemma.
Lemma 4.4.For all K ∈ T h it holds, for all v ∈ H 1 (K), Proof.Once again we use the splitting (4.2), and, by the linearity of Q ∇ K , we have . With this choice, using Proposition 3.2, we can write , where we also used (2.10).
The third commutativity property that we will need in the forthcoming analysis is stated in the lemma below.Lemma 4.5.Let ω ∈ C ∞ ( Ω) be a fixed weight function.Then, for all K ∈ T h , for all w, v ∈ H 1 (K) it holds that (4.12) the implicit constant in the inequality depending on ω 2,∞,Ω .
Proof.Using once more the splitting (4.2), we have . Now, using Proposition 3.3 and Lemma 4.4, we have The above bound also applies to v, finally yielding (4.12).

Negative norm error estimates
The interior error estimate we aim at proving relies on the validity of different bounds on the error measured in negative norms, both at the global and at the local level.We devote this section to study such bounds.

5.2.
Negative norm error estimates on smooth domains.Before going on, let us consider what happens instead in the case in which Ω is smooth.In such a case, for the solution of (5.2), we have that for all p ≥ 0, ϕ ∈ H p (Ω) implies v ϕ ∈ H p+2 (Ω).We can take advantage of such a fact, provided that, for the discretization, we use a tessellation allowing curvilinear edges at the boundary of Ω, and resort, for boundary adjacent elements, to the VEM with curved edges as introduced in [11], following which, we modify the definition of the space B k (∂K), which, for K with |∂K ∩ ∂Ω| > 0 becomes B k (∂K) = {v ∈ C 0 (∂K) : v| e ∈ P k (e) for all edge e of K interior to Ω, v| ∂K∩∂Ω = 0}.
We also modify the interpolator I h by requiring, for boundary adjacent elements, that I h u = 0 on ∂K ∩ ∂Ω.Using arguments similar to the ones used in Section 4 we can see that, also for boundary adjacent elements, for all u ∈ H 1 0 (Ω) with u| K ∈ H 1+t (K), t ≥ 1 it holds that Then, letting, also for elements with a curved boundary edge, Π ∇ K and Q ∇ K be defined by (3.11) and (3.12), we have, for [11]), under the same assumptions on u, we have that (5.8) |u − Π ∇ K u| 1,K h t |u| 1+t,K .
Both Proposition 3.3 and Proposition 3.4 also hold.Indeed, if u ∈ H 1+s (K), 1 ≤ s ≤ k, with u = 0 on ∂K ∩ ∂Ω we have , where we used a Poincaré inequality.By interpolation we have that Proposition 3.3 follows by a triangle inequality.As its proof essentially relies on Proposition 3.3, Proposition 3.4 also follows.
Remark 5.1.We recall that if Ω is a square, the same smoothness result holds for the solution v ϕ of problem (5.2), as in the case of smooth domains, namely, if ϕ ∈ H p 0 (Ω) it holds that v ϕ ∈ H p+2 (Ω) and v ϕ p+2,Ω ϕ p,Ω .Then, also in such a case we have that (5.9) holds.

Local negative norm error estimates.
For Ω being either a polygonal domain, or a smooth domain discretized by elements where, exclusively for elements adjacent to the boundary, curved edges are allowed, we now prove some bounds on the local error, measured in negative norms.We start by proving the following lemma.Lemma 5.2.Let Ω 0 ⊂⊂ Ω 1 ⊂⊂ Ω be fixed smooth interior subdomains of Ω, and let the solution u to Problem (3.1) satisfy u ∈ H t+1 (Ω 1 ), with 1 ≤ t ≤ k.Let u h ∈ V h denote the solution to Problem (3.10).Then, there exists an h 0 such that, if h < h 0 we have, for p ≥ 0 integer, with γ = min{p + 1, k}, and τ = min{p + 1, m}.
Proof.Let Ω with Ω 0 ⊂⊂ Ω ⊂⊂ Ω 1 , be a fixed intermediate subdomain between Ω 0 and Ω 1 , and let ω ∈ C ∞ 0 (Ω ) with ω = 1 in Ω 0 .We let h 0 be such that for all h < h 0 , all elements K ∈ T h with K ∩ Ω = ∅ satisfy K ⊂ Ω 1 , and we let h < h 0 .Letting e = u − u h , we have where v ϕ is the solution to Observe that, as we assumed that Ω Using the error equation (3.17) for v h ∈ V h with supp v h ⊆ Ω 1 arbitrary, we can write (5.11) We now let v h = I h (ωv ϕ ) and we remark that, as supp(ωv We bound the five terms on the right hand side of (5.11) separately.Using (4.1) we can write, with γ = min{p + 1, k} Adding and subtracting ωv ϕ and using (4.1) and (3.18) we have Moreover, using the fact that f − Π 0 m f is orthogonal to P m (T h ) Π 0 m (I h (ωv ϕ )), we can write, with Ω h ⊆ Ω 1 denoting the union of elements K ∈ T h such with I h (ωv ϕ ) = 0 in K, where we added and subtracted ωv ϕ − Π 0 m (ωv ϕ ).Using (2.8) and (4.1) gives (we recall that f = −∆u so, under our assumptions, we have that f ∈ H t−1 (Ω 1 )).By once again adding and subtracting ωv ϕ and using (4.5) and (3.18) we have Finally, we bound V as in [37] as The thesis follows from the observation that τ ≤ γ and f t−1,Ω 1 u t+1,Ω 1 .
Remark 5.4.In the special case in which the right hand side of (3.2) is computable for all v h ∈ V h with supp v ⊆ Ω1 (this happens if f | Ω 1 is a polynomial of degree at most m) then the term III in the sum at the right hand side of (5.11) vanishes.In such a case, we have a better estimate, namely A recursive application of Lemma 5.2 yields the following lemma.
We can now prove Theorem 6.1.
In order to obtain an explicit a priori estimate on the local error we finally combine (6.1) with the global negative norm error estimates of Sections 5.1 and 5.2.We distinguish two cases: Ω polygon and Ω smooth.If Ω is a polygon, we can use the bound (5.7), and, in (6.1), choose p = k − 1.We immediately obtain the following corollary, where σ 0 is the domain dependent parameter defined in Section 5.1, related to the interior angles of Ω. Corollary 6.4.Let Ω 0 ⊂⊂ Ω 1 ⊂⊂ Ω, with Ω polygonal domain, and let u and u h denote, respectively, the solution to (3.2) and (3.10).Assume that f ∈ H ρ (Ω) and u ∈ H 1+r (Ω) for some ρ and r with 0 ≤ r ≤ k and max{0, r − 1} ≤ ρ ≤ m + 1. Assume also that u| For r = ρ = 0 we obtain the following bound, valid under the minimal global regularity assumptions on u, If, on the other hand, Ω is smooth (or if Ω is a square) using once again (6.1) with p = k − 1 yields the following bound: This time, we have the following corollary.Corollary 6.5.Let Ω 0 ⊂⊂ Ω 1 ⊂⊂ Ω, with Ω smooth domain, and let u and u h denote, respectively, the solution to (3.2), and the solution to (3.10) obtained with the discretization considered in Section 5.2.Assume that u ∈ H 1+r (Ω) for some r with 1 ≤ r ≤ k and that u| Ω 1 ∈ H 1+t (Ω 1 ), with r ≤ t ≤ k.Then we have u − u h 1,Ω 0 h κ , with κ = min{t, m + r}.
Under the minimal global regularity assumption on f , namely f = −∆u ∈ L 2 (Ω), this time we have that u − u h 1,Ω 0 h t (|u| 1+t,Ω 1 + u 1,Ω ) + h min{t,m}+1 f 0,Ω .Observe that we do not have optimality unless m ≥ k − 1. Remark 6.6.While, for the sake of simplicity, we focused our analysis on homogeneous Dirichlet boundary conditions, the result presented in Section 6 extends also to other boundary conditions, such as non homogeneous Dirichlet, or mixed.In such cases, depending on the smoothness of the boundary data, the solution u might lack overall regularity also for very smooth right hand side f .6.1.Extension to the enhanced virtual element method.Let us briefly sketch how the interior estimate (6.1) can be extended to a particularly relevant version of the virtual element method, namely the enhanced VEM [1], characterized by a local discretization space defined as , where the orthogonality is intended with respect to the L 2 (K) scalar product, and where P k (K) ∩ P k−2 (K) ⊥ denotes the space of those polynomials of degree at most k that are orthogonal, in L 2 (K), to all polynomials of degree at most k − 2. Letting V en h ⊂ H 1 0 (Ω) denote the corresponding global virtual element space, we consider the problem: find u h ∈ V en h such that for all v h ∈ V en h (6.7) where a h is defined as before.The space V en h does not fall into the framework which we considered up to now, since for no value of m we have that V k en (K) = V k m (K).As a consequence, the proof of Lemma 4.1 is not valid for such a space.However we know (see [1]) that the functions in V k en (K) and in V k k−2 (K) have the same set of degrees of freedom, and that, letting v h , w h ∈ V k en (K) and v h , w h ∈ V k k−2 (K) denote two couples of functions satisfying (6.8) (which is equivalent to saying that the value of all the degrees of freedom of v h , w h coincide, respectively, with those of v h and w h ) we have (6.9) It is then not difficult to check that u h ∈ V en h is the solution of (6.7) if and only if the corresponding function u h ∈ V h (V h being the "plain" VEM space defined in (3.4) where the enhanced projection Π 0 en : Apart from the definition of the right hand side, equation (6.10) falls in the framework studied in the previous sections.It is not difficult to check that for f ∈ Thanks to this inequality, used for those bounds affected by the altered right hand side, particularly Proposition 3.5 and Lemma 5.2, our analysis carries over, with minor modifications, to Problem (6.10).Moreover, thanks to the higher approximation order of Π 0 en with respect to Π 0 k−2 (Π 0 en satisfies an error bound similar to the one in Proposition 3.2, [1]) Lemma 5.2 now holds with τ = min{p + 1, k}.Then Theorem 6.1, as well as its corollaries, hold, and provide optimal local error bounds for u − u h .More precisely, letting Ω ⊂⊂ Ω 1 ⊂⊂ Ω, for h small enough, under the minimal global regularity assumptions on u and f , we have that u ∈ H t+1 (Ω 1 ), with t ≤ k, implies (6.11) Let now Ω 0 ⊂⊂ Ω and assume that h is sufficiently small so that Ω + 0 ⊂ Ω , where ∅} being the set of all elements that have non empty intersection with Ω 0 .To bound u − u h in Ω 0 we start by observing that, by triangle inequality and (5.8) we have To bound the second term on the right hand side, we add and subtract, element by element, the boundary average, apply a triangle inequality and use a Poincaré inequality to bound the L 2 norm of the boundary-average free terms with their H 1 seminorm, which, in turn, is bound using (3.8),  thus obtaining ,Ω , where we could replace u h with u h thanks to (6.8) and (6.9) and where we used (3.8) once again.Combining (6.12) with (5.8) and (6.13) finally yields the optimal error bound (6.13) u − u h 1,Ω 0 h κ (|u| 1+t,Ω 1 + f 0,Ω ), k = min{t, σ 0 − ε} if Ω is a polygon, t if Ω is smooth.
The solution u has a singularity in (0, 0), and, for both test cases, it verifies u ∈ H s (Ω), for all s < 3/2, but u ∈ H 3/2 (Ω).Consequently, we expect the global H 1 (Ω) error not to converge faster than h 1/2 .On the other hand the solution is smooth in a neighborhood of Ω 0 .According to Corollaries 6.5 and 6.4 we can then expect, for Test 1 and Test 2 respectively, a convergence rate of order k and min{k, 4/3 − ε} (ε arbitrarily small).We solve the problem by the enhanced virtual element method (see Section 6.1) with k = 1, • • • , 4. For both test cases we consider both a sequence of progressively fines structured hexagonal meshes and a sequence of progressively finer regular Voronoi meshes.Examples of the meshes used for the numerical tests are displayed in Figures 1 and 2. For all discretizations the stabilization is chosen to be the simple so called dofi-dofi stabilization, that, under our mesh regularity assumptions, is optimal.In all figures, we display, for reference purpose, dotted straight lines with a slope corresponding to the expected convergence rate for global and local error, namely 1/2 for the global error and, respectively, k and min{4/3, k} for the local error in the squared and in the L-shaped domain.

Figure 1 .
Figure 1.Examples meshes of the unit square.

Figure 2 .
Figure 2. Examples meshes of the L-shaped domain.

Figure 5 ..
Figure 5. Squared domain, Voronoi meshes.On the x axis the meshsize h, on the y axis the global (blue squares) and local (yellow/orange diamond) H 1 errors e 1 Ω , e 1 Ω − 0

Figures 3
Figures 3 through 6 clearly confirm the validity of the a priori estimate.In particular the local error behaves as expected both in the squared domain case and in the L-shaped domain case.