Higher-order adaptive virtual element methods with contraction properties

The realization of a standard Adaptive Finite Element Method (AFEM) preserves the mesh conformity by performing a completion step in the refinement loop: in addition to elements marked for refinement due to their contribution to the global error estimator, other elements are refined. In the new perspective opened by the introduction of Virtual Element Methods (VEM), elements with hanging nodes can be viewed as polygons with aligned edges, carrying virtual functions together with standard polynomial functions. The potential advantage is that all activated degrees of freedom are motivated by error reduction, not just by geometric reasons. This point of view is at the basis of the paper [L. Beirao da Veiga et al., Adaptive VEM: stabilization-free a posteriori error analysis and contraction property, SIAM Journal on Numerical Analysis, vol. 61, 2023], devoted to the convergence analysis of an adaptive VEM generated by the successive newest-vertex bisections of triangular elements without applying completion, in the lowest-order case (polynomial degree k=1). The purpose of this paper is to extend these results to the case of VEMs of order k>1 built on triangular meshes. The problem at hand is a variable-coefficient, second-order self-adjoint elliptic equation with Dirichlet boundary conditions; the data of the problem are assumed to be piecewise polynomials of degree k-1. By extending the concept of global index of a hanging node, under an admissibility assumption of the mesh, we derive a stabilization-free a posteriori error estimator. This is the sum of residual-type terms and certain virtual inconsistency terms (which vanish for k=1). We define an adaptive VEM of order k based on this estimator, and we prove its convergence by establishing a contraction result for a linear combination of (squared) energy norm of the error, residual estimator, and virtual inconsistency estimator.


Introduction
Adaptive Finite Element Methods (AFEM) for self-adjoint coercive problems written in the form ∀v ∈ V iterate the sequence SOLVE → ESTIMATE → MARK → REFINE to produce better and better approximations of u.Their practical efficiency is corroborated by sound theoretical results of convergence, complexity, and optimality, which in various cases (such as, e.g., conforming h-versions) completely explain the behaviour of the adaptive algorithms [15,11,14,16,13].
The standard AFEM realization preserves the conformity of the initial mesh, at the expense of performing a completion step in REFINE: in addition to elements marked for refinement due to their contribution to the global error estimator, other elements are refined.Without this step, one would obtain nonconforming meshes, containing elements with hanging nodes.
In the new perspective opened by the introduction of Virtual Element Methods (VEM) [3,4], elements with hanging nodes can be viewed as polygons with aligned edges, carrying virtual (i.e., non-accessible) functions together with standard polynomial functions.The potential advantage is that all activated degrees of freedom are motivated by error reduction, not just by geometric reasons.On the other hand, in this transformation of an adaptive FEM into an adaptive VEM, one looses the availability of a general convergence theory, which so far is lacking (although results on a posteriori error estimates [8,12] have been obtained, together with efficient practical recipes for refining polytopal meshes [9,10,2]).
Such a shift in perspective inspired the recent papers [5,6], devoted to the analysis of an adaptive VEM generated by the successive newest-vertex bisections of triangular elements without applying completion, in the lowest-order case (polynomial degree k = 1).Despite the simple geometric setup, the investigation faced some VEM-specific obstacles in the analysis, giving answers that could prove useful in the study of more general adaptive VEM discretizations.For instance, a VEM solution u T ∈ V T ⊂ V, defined by the Galerkin projection satisfies an a posteriori error bound of the type T (u T ) + S T (u T , u T ) , where η T (u T ) is a residual-type error estimator, S T (u T , u T ) is the stabilization term that makes the discrete bilinear form B T (u T , v T ) coercive in V, and for simplicity we assume piecewise constant data on the mesh T .Unfortunately, the term S T (u T , u T ) need not reduce under a mesh refinement, as η 2 T (u T ) does: this makes the convergence analysis problematic.However, one of the key results obtained in [5] states that S T (u T , u T ) is dominated by η 2 T (u T ), i.e. S T (u T , u T ) η 2 T (u T ) , provided an assumption of admissibility of the non-conforming meshes generated by successive refinements is fulfilled; such a restriction, which appears to have little practical impact, amounts to requiring the uniform boundedness of the global index of all hanging node, a useful concept introduced in [5] to hierarchically organize the set of hanging nodes.Once the a posteriori error bound is reduced to T (u T ) , the convergence analysis becomes feasible, and a contraction property is proven to hold for a linear combination of the (squared) energy norm of the error and the (squared) residual estimator.
The purpose of this paper is to extend the results in [5] to the case of VEMs of order k ≥ 2 built on triangular meshes.Note that the interest in avoiding the creation of new elements just to satisfy the conformity condition of the mesh becomes more and more evident as the polynomial degree increases.The problem at hand is again a variablecoefficient, second-order self-adjoint elliptic equation with Dirichlet boundary conditions.The geometric concept of hanging node (a vertex for some elements, contained inside an edge of some other elements) is replaced by a functional one, referring to the degrees of freedom associated with the node; once the meaning of hanging node is clarified, the definition of global index of a node, and its role in the analysis, is similar to the one given in [5].
A significant difference with respect to the content of that paper concerns the control of the stabilization term, which does not involve only the residual estimator, but a new term, called the virtual inconsistency estimator and denoted by Ψ T (u T ).It measures the projection error, upon local spaces of polynomials, of certain expressions depending on the operator coefficients and the discrete solution; it vanishes when k = 1 or when the coefficients are constant.The new stabilization bound, which we derive under an admissibility assumption of the mesh, takes the form , which leads to the a posteriori, stabilization-free error control . Correspondingly, we obtain the convergence of the adaptive VEM of order k by proving a contraction result for a linear combination of (squared) energy norm of the error, (squared) residual estimator, and (squared) virtual inconsistency estimator.
Similarly to [5], we assume here that the data D of our boundary-value problem are piecewise polynomials of degrees related to k − 1, on the initial mesh T 0 and consequently on each mesh T derived by newest-vertex bisection.This is not a restriction, since we propose to insert the adaptive VEM procedure just described, which we now consider as a module GALERKIN, into an outer loop AVEM of the form where the module DATA produces, via greedy-type iterations, a piecewise polynomial approximation of the input data with prescribed accuracy, defined on a suitable refinement of the input partition.Manifestly, the target accuracy is matched after a finite number of calls to DATA and GALERKIN.Properties of complexity and quasi-optimality of this two-loop algorithm are investigated in [6] in the linear case k = 1.We plan to do the same for the case k ≥ 2 in a forthcoming paper.
The outline of this paper is as follows.In Sections 2 and 3, we introduce the model boundary-value problem, and its discretization by an enhanced version of the VEM ( [1]).In Section 4 we define the global index of a node, and we formulate the admissibility assumption on the mesh.Two essential properties for bounding the stabilization term are established in Section 5.The a posteriori error estimators are defined in Section 6, whereas stabilization-free a posteriori error estimates are proven in Section 7. In Section 8, we investigate how the a posteriori error estimators are reduced under mesh refinement.These properties are needed to justify the refinement strategy in our adaptive module GALERKIN, which is described in Section 9.The paper ends with the proof of convergence of the loop GALERKIN, reported in Section 10.

VEM spaces of order k ≥ 2
We consider the following Dirichlet boundary value problem in a polygonal domain Ω, where A ∈ (L ∞ (Ω)) 2×2 is symmetric and uniformly positive definite in Ω, c ∈ L ∞ (Ω) and non-negative in Ω, f ∈ L 2 (Ω).Data will be denoted by D = (A, c, f ).The variational formulation of this problem is written as where (•, •) is the scalar product in L 2 (Ω) and B(u, v) := a(u, v) + m(u, v) is the bilinear form associated with Problem (2.1), i.e, We denote the energy norm as In order to find a discrete approximation of the solution of Problem (2.2), we firstly introduce a fixed initial partition T 0 on the domain Ω made of triangular elements E. We will denote by T any refinement of T 0 obtained by a finite number of newest-vertex element bisections.We underline that we are not requiring T to be a conforming mesh, since hanging nodes may arise in the refinement.The classification of nodes, which will play a crucial role in the proofs presented in this paper, is postponed in Section 4.
According to the Virtual Element theory [3], an element E of the triangulation can be viewed as a polygon with more than three edges, if some hanging nodes are sitting on its boundary.We can then denote by E E the set of edges e of element E and E := E∈T E E .We finally define the diameter of an element E as h E = |E| 1/2 and h = max E∈T {h E }.
We introduce the functional spaces needed to apply the Virtual Element Method (VEM).We start by defining the space of functions on the boundary of E, V ∂E,k , which is constituted by the functions that are continuous on the boundary of E and that, when restricted to any edge of ∂E, are polynomials of degree k > 0, i.e, V ∂E,k := {v ∈ C 0 (∂E) : v| e ∈ P k (e), ∀e ⊂ ∂E}.
Then, we define the "enhanced" VEM space in E, as done in [1], such that where is the projector defined by We remark that V E,k contains the polynomial space of degree k on E and its dimension is since in our analysis we consider triangular elements.We notice that in the case k > 1 a function v in V E,k is uniquely defined by • the set of the values at the vertices of E; • the set of the values at the k − 1 equally-spaced internal points on each edge of ∂E; • the set of the moments where the set M p (E), p ≥ 0, is defined as We will denote by µ p (E, v) = 1 |E| E v(x)m(x)dx : m ∈ M p (E) \ M p−1 (E) the vector of the moments of v of order p.By |µ p (E, v)| we will denote the l 2 -norm of this vector.
We can now introduce the global discrete space as On T we need also to give the definition of the space of piecewise polynomial functions on T and its subspace which plays a crucial role in the forthcoming analysis.
We now introduce a series of projectors that will be used in the rest of the paper.For any E ∈ T , we denote by Π 0 p,E : L 2 (E) → P p (E) the L 2 (E)-orthogonal projector onto the space of polynomial of degree p on E. Thanks to the choice of the enhanced space V E,k (2.4), we remark that Π 0 k,E v and Π 0 k−1,E ∇v can be computed for any function v ∈ V E,k , see [1] for the details.To simplify the notation, in the following we will drop the symbol E from Π 0 k,E when no confusion arises.The global L 2 -orthogonal projector is denoted by Π 0 p,T : L 2 (Ω) → W p T .We can also define the Lagrange interpolation operator I E : V E,k → P k (E) on E, which builds a polynomial of degree k using the 3k degrees of freedom on the boundary of E and the moments of order ≤ k − 3, since Moreover, we will denote by I T : V T → W k T the Lagrange interpolation operator that restricts to I E on each E ∈ T .
3 Discretization with data of degree k − 1 In the rest of this paper, we assume that data D = (A, c, f ) are piecewise polynomials of degree k − 1 on the initial partition T 0 , hence on each partition T obtained by newest-vertex refinement.Their values on each element of the triangulation will be denoted by We here define the bilinear forms that we need for the Galerkin discretization problem, starting from a E , m E : We also introduce the symmetric bilinear form s E : where {x i } N E i=1 indicates the set of the degrees of freedom on the boundary of E. Indeed, we remark that in this case the stabilization term can be built without using the internal degrees of freedom, as shown in [7].We assume for s E the existence of two positive constant c s and C s independent on E, such that We define the local stabilizing form as ∀v, w ∈ V E , and the global stabilization form From (3.1), we get where | • | 1,T denotes the broken H 1 -seminorm over T .Thus, we can now define the bilinear form B T (•, •), B T : ) with γ independent of T satisfying γ ≥ γ 0 for some fixed γ 0 > 0. For the loading term we introduce F T : V T → R as since f E has been already approximated with a polynomial of degree k − 1.Note that the equality in (3.3) remains true if f E is an approximation of f of degree k on E.
We have now defined all the forms that appear in the discrete formulation of the Problem (2.2).It reads as find The bilinear form B T is continuous and coercive, hence, there exists a unique and stable solution of the Problem (3.4).Furthermore, the following result extends Lemma 2.6 in [5].Lemma 3.1 (Gakerkin quasi-orthogonality).For any v ∈ V T and w , where u is the solution of (2.2) and u T the solution of (3.4).
Figure 1: Blue squares represent the k + 1 equally-spaced nodes ξ n on the edge S before refinement.Red circles represent the 2k + 1 nodes that arise after refinement.We have denoted by ζ i the new nodes that do not coincide with any ξ n .

The index of a node
A crucial concept, firstly introduced in [5] for the case k = 1, is the global index of a node: it will be used in the proofs of Section 5.In order to extend its definition to the case k > 1, we preliminarily introduce some useful definitions.
Considering the affine function F E : Ê → E mapping the reference element onto an element E ∈ T , we define the physical lattice on E by R E,k := F E ( R Ê,k ), and the set of proper nodes of E as the points of the physical lattice sitting on the boundary of E, i.e., Observe that we implicitly assume that k ≥ 2 is sufficiently small so that interpolation on equally spaced nodes is numerically stable.
Next, we denote by H E the set of hanging nodes of E, i.e., the set of points x ∈ ∂E that are not proper nodes of E, but that are proper nodes of some other element E ′ , i.e., Finally, let N E := P E ∪ H E be the set of all nodes sitting on E.
At the global level, N := E∈T N E will be the set of all nodes of the triangulation T , which we split into the set P := {x ∈ N : x ∈ P E ∀E containing x} of the proper nodes of T , and the set H := N \ P of the hanging nodes of T .
Next, let us clarify what happens when a hanging node is created.Let S be an element edge that is being refined, i.e., split into two contiguous edges S − and S + .Before the refinement, S contains k + 1 equally-spaced nodes ξ n , n = 1, . . .k + 1: the endpoints and the k − 1 internal ones.After the refinement, S contains 2k + 1 nodes, precisely k + 1 equally-spaced nodes on each sub-edge S ± , with the midpoint in common; see Figure 1.The spacing of the 'old' nodes on S was |S| k (where |S| denotes the length of S), whereas the spacing of the 'new' nodes is |S| 2k .Consequently, k + 1 of these nodes coincide with those initially on S, and the new nodes introduced in the refinement are only k.We will denote these latter by ζ i , i = 1, . . ., k.This suggests the following definition.

Definition 4.1 (closest neighbors of a node). With the previous notation, if
We are ready to give the announced definition of global index of a node of the triangulation T .Definition 4.2 (global index of a node).Given a node x ∈ N , we define its global index λ recursively as follows: • if x is a proper node, then λ(x) := 0; (a) Figure 2: Triangulation after the three refinements in the case k = 2 (a) and in the case k = 3 (b).Blue crosses represent the original degrees of freedom.Red squares, green circles and orange triangles are used for the degrees of freedom of the first, second and third refinement, respectively.All nodes are proper, except those on the horizontal line, whose global index is reported.
Figure 2 shows the evolution of the global index after three refinements in the cases k = 2 (a) and k = 3 (b).We remark that, for instance, the midpoint of the horizontal edge is a proper node in case (a), and a hanging node in case (b).
The largest global index in T will be denoted by Λ T := max x∈N {λ(x)}.In this paper, as in [5], we will consider sequences of successively refined triangulations {T } whose global index does not blow up.

Assumption 4.3.
There exists a constant Λ > 0 such that, for any triangulation T generated by successive refinements of T 0 , it holds Any such triangulation will be called Λ-admissible.

Two key properties
In this section we discuss the validity of some results for the degree k > 1 that will be used in the rest of the paper.We will highlight in particular the differences from the case k = 1.

Proposition 5.1 (scaled Poincaré inequality in V T ). There exists a constant
Proof.Let E ∈ T be an element of the triangulation.If E is an element of the original partition T 0 , all its vertices are proper nodes.Otherwise, E has been generated after some refinements by splitting an element E into two elements, E and E ′ .Let L be the common edge shared by E and E ′ .If L is not further refined, then all the nodes on L are proper because they are shared by E and E ′ .If L is refined and k is even, then the midpoint of L is a proper node.
So, let us consider the case k odd and let us assume that L is refined M ≥ 1 times.We focus in particular on the internal node x of L is at distance |L| k from one of the endpoints, Figure 3 shows the case k = 3.This point belongs to one of the M + 1 intervals in which L is refined, having width |L|/2 s , for some 1 ≤ s ≤ M .We remark that s depends on how L has been refined (in the case of uniform refinements of L, one has 2 s = M + 1).We localize the chosen node x in L by defining an m ≥ 0 such that or, equivalently, The interval going from |L| m 2 s to |L|(m+1) is an edge for a smaller element E ′ , thus it contains k − 1 internal nodes.Since they are equi-spaced, their positions are at with n = 0, . . ., k.
By taking n = 2 s − m k, which is compatible with conditions (5.2), we conclude that one of the internal nodes of E ′ coincides with x.
This guarantees that E has at least one proper node x on its boundary.By hypothesis v(x) = 0, and so we can apply the classical Poincaré inequality, 1,E , that concludes the proof.
Remark 5.2.The previous proof exploits the fact that when k > 1, each element of the triangulation contains at least a proper node.This differs from the case k = 1 in which the edges do not contain internal nodes, and then elements with all hanging nodes as vertices are admissible.As a further difference from the case k = 1, we highlight that in Proposition 5.1 the constant C P does not depend on the constant Λ, whose existence has been introduced in Assumption 4.3.
The next result we are going to establish is a hierarchical representation of the interpolation error v − I E v on the boundary ∂E of an element E ∈ T .Assume that v ∈ V E,k , and let L be an edge of E; for simplicity, in the sequel the restriction of v to L, which is a piecewise polynomial of degree k, will be still denoted by v.The subsequent bisections of L which generate the nodes in N E ∩ L allow us to write the difference (v − I E v) |L telescopically as (5.3) here, I 0 = I E |L , I JL is the identity operator, whereas I j v for 1 ≤ j ≤ J L − 1 is the piecewise polynomial of degree k which interpolates v on the partition of L of level j, namely the partition formed by sub-edges of length ≤ |L| 2 j .In order to understand the structure of the detail (I j − I j−1 )v, assume that S is a sub-edge of L of length = |L| 2 j−1 , which is split into two sub-edges S ± of length = |L| 2 j (see again Fig. 1).On S we have two interpolation operators, namely I := I j−1 |S : C 0 (S) → P k (S) and ) , which coincides with the interpolation operator I − : C 0 (S − ) → P k (S − ) when restricted to S − and with the analogous operator I + when restricted to S + .
With the notation introduced just before Definition 4.1, we can quantify the discrepancy between the two interpolation operators by defining the k basis functions See Figure 4 for a graphical representation of these functions in the cases k = 1 (a), k = 2 (b), k = 3 (c).Hence, the difference between the two interpolation operators on S can be written as where d is defined as (5.4) The values of Iv at the k nodes ζ i are a linear combination of the values of Iv at the k + 1 nodes ζ n , where Iv coincides with v. Thus, there exist coefficients α i,n such that The explicit values of these coefficients in the case k = 2 for the two new nodes ζ 1 and ζ 2 are given in these expressions: , Similarly, in the case k = 3, we get where again ξ i ≤ ζ i ≤ ξ i+1 , i = 1, 2, 3. Figure 5 shows both cases.We notice that the coefficients α i,n depend only on the relative positions of the nodes on S, not on the level j of refinement.
Summarizing, at the level j of refinement of the edge L, we get where H L,j is the set of hanging nodes on L created at the level j of refinement, whereas Summing-up over the levels and recalling (5.3), we obtain where We now introduce the subspace of V E,k On X E , we have two norms, namely the seminorm |w| 1,E (which is a norm on X E due to the vanishing of w at the three vertices of E) and the norm Note that, due to Assumption 4.3, the dimension of X E is uniformly bounded by a constant depending on Λ; furthermore, the number of possible patterns of hanging nodes on ∂E, which determine the details d(w, x), is also bounded in terms of Λ.As a consequence, the two norms are equivalent, with equivalence constants depending on Λ.Therefore, Summing-up over all the elements of the triangulation, we arrive at the following result.

Lemma 5.3 (global interpolation error vs hierarchical errors).
There exists a constant C D > 0 depending on Λ but independent of the triangulation T such that Next, we introduce the interpolation operator ) where V 0 T is defined in (2.8), by the following conditions: ) for all 0 ≤ p ≤ k − 3 and for all E ∈ T .
These conditions uniquely identify I 0 T v. Indeed, if x ∈ H is generated by a refinement of level j of an edge L (say, x = ζ i with the notation introduced before Definition 4.1), then (I 0 T v)(x) can be expressed in terms of the values of I 0 T v at the k + 1 nodes (say, ξ n ) created at the previous levels of refinement of L, using the same coefficients as in formula (5.5), i.e., and so on recursively.
The following result provides a representation of the error Proof.Consider an element E ∈ T .Recall that by construction it holds µ p (E, Summing on all the elements of the partition, we get . This concludes the proof. Concatenating Lemma 5.3 and Lemma 5.4, we can prove the second key property of this section.Proposition 5.5 (comparison between interpolation operators).Let T be Λ-admissible.Then, there exists a constant C I > 0, depending on Λ, but independent of T , such that Proof.Given a function v ∈ V T , by the triangle inequality T , so it is enough to bound the last norm on the right-hand side.To this end, considering the vectors and recalling the two Lemmas, the proof can be concluded if we show that Given x ∈ H, assume that it is generated by a refinement of level j of an edge L (say, x = ζ i with the notation introduced before Definition 4.1).Writing v * := I 0 T v for short, and exploiting formulas (5.4) and (5.5), we get (5.9) Thus, we can build a matrix W : l 2 (H) → l 2 (H) such that δ = W d, and we just need to prove that We now organize the hanging nodes with respect to the global index λ ∈ [1, Λ T ].Calling H λ = {x ∈ H : λ(x) = λ}, and H = 1≤λ≤ΛT H λ , the matrix W can be factorized in lower triangular matrices W λ , that change the nodes of level λ, leaving the others unchanged.In particular, where W 1 is just the identity matrix I, whereas each other matrix W λ differs from the identity only in the rows of block λ.In each of these rows, all entries are zero, but the entries α i,n in the off-diagonal part and 1 on the diagonal.In order to estimate W λ , we use the Hölder inequality where in the last inequality it has been used the fact that a hanging node of global index < λ may appear at most 5 times on the right-hand side of (5.9), since at most five edges meet at a node [5, Proposition 3.2].These bring us to the following bound . and the proof is concluded.

A posteriori error estimator
With the aim of discussing the a posteriori error analysis, and following [12], we define the a posteriori error estimators, starting from the internal residual over an element E, i.e., for any v ∈ V E,k .We highlight that in the case k = 1, with piecewise constant data, the diffusion term in the residual vanishes.Furthermore, we define the jump residual over e, where e is an edge shared by two elements E 1 and E 2 of the partition T , as where n i denotes the unit normal vector to e pointing outward with respect to E i ; we set j T (e; v) = 0 of e ∈ ∂Ω.Then, let the local residual estimator associated with E be and the global residual estimator as the sum of the local residuals In contrast to what has been done for the case k = 1, we also need to introduce the virtual inconsistency terms, defined by as well as their sum

A posteriori error estimates
In this section we present one of the main results of this paper, a stabilization-free a posteriori error bound.In this view, we firstly start by introducing the classical Clément operator upon the space V 0 T , Ĩ0 T : V → V 0 T ; it is defined at the proper nodes on the skeleton of T as the average of the target function on the support of the associated basis functions, whereas the internal moments (if any) coincide with those of the target function.
The scaled Poincaré inequality (Proposition 5.1) and Proposition 5.5 guarantee the validity of the error estimate for Ĩ0 T .Given these propositions, its proof does not involve the polynomial degree k, hence, it does not change with respect to the one presented in [5].
where the hidden constant depends on Λ but not on T .
We can now prove the following results, which is similar to Theorem 13 in [12], but with a slightly modified proof.

Proposition 7.2 (upper bound).
There exists a constant C apost > 0, independent of u, T , u T and γ, such that Proof.For any v ∈ V, using the definition of Problem (2.2), we have that The first term can be written as The addend I 1 can be expressed as which can be bounded by using Lemma 7.1, On the other hand, noting that and applying again Lemma 7.1, the addend I 2 can be bounded as follows: Looking now at the term II, we have by Lemma 3.1 Finally, by taking v := u − u T ∈ V, we get which, using the coercivity of B, concludes the proof.
We now report a bound for the local residual estimator, proved in [12](Theorem 16).

Proposition 7.3 (local lower bound).
There holds where The hidden constant is independent of γ, h, u and u T .
Summing on all the elements of the partition, we get the following corollary.

Corollary 7.4 (global lower bound).
There exists a constant c apost > 0, independent of u, T , u T and γ, such that In the following proposition we present a bound of the stabilization term.We remark that in the case k = 1 the inconsistency term does not appear.Proposition 7.5 (bound of the stabilization term).There exists a constant C B > 0 independent of T , u T and γ, such that Proof.From the definition (3.2) of the form B T and from (3.4), Defining e T := u T − w, we get We notice that By substituting (7.5) and (7.6) into (7.4), it results With the same strategy used in [5], for any δ > 0, we get where Posing now w = I 0 T u T and applying Proposition 5.1, we get Φ T (u ,Ω , whereas Proposition 5.5 yields Combining Propositions 7.2 and 7.5, we arrive at the following key result.Corollary 7.6 (stabilization-free a posteriori error upper bound).It holds where 8 The effect of a mesh refinement In view of the convergence analysis of the adaptive algorithm GALERKIN, in this section we analyse the effect of refining the partition T by applying one or more newest-vertex bisections to some of its elements.Specifically, in Sect.8.1 we prove that the residual estimator (6.2) is reduced by a fixed fraction (up to an addend proportional to the stabilization term) when the element E is split into two elements by one bisection.We prove a similar result for the inconsistency term estimator (6.4), provided a suitable number of bisections is applied to E. Next, in Sect.8.2 we establish a quasi-orthogonality property in the energy norm between the solutions on two partitions, one being a refinement of the other.

Reduction of estimators under refinement
Let us consider an element E in T which is bisected into elements E 1 and E 2 ; the refined partition containing these two elements will be denoted by T * .Given v ∈ V T , we notice that v is known on ∂E, and in particular at the new vertex of E 1 and E 2 produced by the bisection.Denoting by e = E 1 ∩ E 2 the new edge, we associate a function In the following we will write v instead of v * when no confusion arises.

The residual estimator
Let η T (E; v, D) be defined in (6.2) and η T * (E; v, D) be the sum of the local residual estimators on the two newly formed elements , defined as follows: where we recall that We notice that, since D does not change under refinement, the functions will be denoted again by f E , c E and A E , respectively.Lemma 8.1 (local residual estimator reduction).There exist constants µ r ∈ (0, 1) and c er,1 > 0 such that for any Proof.Recalling the definition (6.1), we have the following residuals The second term can be bounded by using the inverse inequality and the minimality of Π 0 k−1,Ei as follows: while, for the last term, using the Poincaré inequality we have Finally, taking an appropriate value of ǫ and setting µ := 1+ǫ 2 ∈ (0, 1) (for instance, if ǫ = 1 2 , µ = 3 4 ) we get where C > 0 is a constant.
For the jump condition, we will essentially use the proof given in [5, Lemma 5.2].In particular, we write j T * (e; v) = j T (e; v) + (j T * (e; v) − j T (e, v)) and for any ǫ > 0 with On the new edge we notice that j T (e; v) = 0, then, We now define for any edge e ∈ E Ei , we denote by E i,e ∈ T * (E i ) the element such that e = ∂E i ∩ ∂E i,e .Then, )∇v 0,e , where Êi,e indicates the parent of E i,e .Using the trace inequality we have Using now the minimality property of Π 0 k−1,E ′ and Π 0 k−1, Ê′ , we easily get as above which, for a sufficiently small ǫ, concludes the proof.
From this Lemma and the Lipschitz continuity of the residual estimator with respect to the argument v (whose proof is independent of the used polynomial degree, so we refer to [5, Lemma 5.3]), we immediately deduce the following result.
Proposition 8.2 (residual estimator reduction on refined elements).There exist constants µ r ∈ (0, 1), c er,1 > 0 and c er,2 > 0 independent of T such that for any v ∈ V T and w ∈ V Given v ∈ V T and E ∈ T , consider the two virtual inconsistency terms Ψ T ,A (E, v, D) and Ψ T ,c (E, v, D) introduced in (6.3).When E is bisected into E 1 and E 2 , the term Ψ T ,c (E, v, D) is reduced by a factor µ c < 1 up to an addend proportional to the stabilization term, i.e., there exists c vi,c > 0 such that This stems from the presence of the factor h E in front of the norm , with an argument similar to the one used in the proof of Lemma 8.1.
Due to the lack of the factor h E , a reduction result similar to (8.2) does not hold for Ψ T ,c (E, v, D).Indeed, since A E Π 0 k−1,E ∇v ∈ P 2k−2 (E), one may ask whether a constant µ < 1 esists such that Unfortunately, the answer is no, as it can be seen numerically, working on the reference element Ê by affinity and identifying µ 2 as the largest eigenvalue of a generalized eigenvalue problem.However, the same numerics indicates that if Ê is split into 2 m triangles of equal area by m successive levels of uniform bisections, then µ 2 becomes < 1 for m large enough, as seen in Table 1.This is indeed predicted by the following result.
Lemma 8.3.Let E ∈ T .For any polynomial degree k ≥ 1 there exists a minimal m ∈ N and a constant µ = µ m < 1 independent of E such that, if E is partitioned into 2 m elements E i of equal area by m levels of uniform newest vertex bisection, it holds Proof.Since by construction h Ei = 2 −m/2 h E , classical approximation results give for some constant C k depending on k.Replacing q by q − Π 0 k−1,E q leaves the left-hand side unchanged, whereas on the right-hand side an inverse inequality yields One concludes taking as m the smallest integer such that Based on these results, let T * m be a refinement of T in which the element E has undergone m levels of uniform refinements by newest vertex bisection, and has been replaced by 2 m subelements E i .Given v ∈ V T , let us set Lemma 8.4.There exist constants ρ A < 1 and c vi,A > 0 such that for any v ∈ V E,k Proof.Write sum over i, and conclude using (8.4) and the usual arguments based on the minimality of the L 2 -orthogonal projections.
Let us set Applying a bound similar to (8.2) to the successive level of refinements, we arrive at the following result.Lemma 8.5.There exist constants µ vi < 1 and c vi,1 > 0 such that for any v ∈ V E,k Combining this estimate with the Lipschitz continuity property of the virtual inconsistency estimator, we obtain the following result.Proposition 8.6 (virtual inconsistency estimator reduction on refined elements).There exist constants µ vi ∈ (0, 1), c vi,1 > 0 and c vi,2 > 0 independent of T such that for any v ∈ V T and w ∈ V T * m , and any element Next result extends Corollary 5.8 in [5].Proposition 8.8 (quasi-orthogonality of energy errors without stabilization).Given any δ ∈ 0, 1  4 , there exists γ δ > 0 such that for any γ > γ δ , it holds Combining them, we have Doing the same on T * and defining Employing Proposition 8.7, we obtain By choosing γ such that we get , which concludes the proof by observing that 1+2δ 1−δ ≤ 1 + 4δ and δ 1−δ ≤ 2δ, when δ ≤ 1 4 .

The module GALERKIN
Let us consider a Λ-admissible input mesh T 0 , a set of approximated data D which consist of piecewise polynomials of degree k − 1 on T 0 , and a tolerance ǫ > 0. The call produces a Λ-admissible refined mesh T and the Galerkin approximation u T ∈ V T , such as where u is the solution of Problem (2.2) and C G = c B max {C U1 , C U2 }, with c B is defined in (2.3) and C U1 , C U2 in Corollary 7.6.We obtain it by iterating the sequence At each step, a Λ−admissible mesh T j and the associated solution u j of the discrete Problem (3.4) are produced.The process stops when the condition η 2 Tj (u j , D) + Ψ 2 Tj (u j , D) ≤ ǫ 2 is reached.In particular, the modules are defined as follows: [15] and finds an almost minimal set M of elements in T such that for a given parameter θ ∈ (0, 1); • [T * ] = REFINE(T , M, Λ) returns a Λ-admissible refined mesh obtained from T by suitable newest-vertex bisections of the elements in M, and possibly of other elements to fullfil the Λ-admissibility condition.
It is worth adding some details about the procedure REFINE.Let E ∈ M be an element marked for refinement.For simplicity, hereafter let us set η := η T (E; u T , D) and Ψ := Ψ T (E; u T , D).The refinement of E is performed as follows: • if η ≥ Ψ, then E is bisected once; • if η < Ψ, then E is bisected m-times, where m has been introduced in Sect.
In the case η < Ψ, η * + Ψ * ≤ max(µ m r , µ vi )(η + Ψ) + c S 1/2 T (E) (u T , u T ) .In all cases, it holds which shows that in each marked element the sum of the two estimators is reduced under refinement, up to the stabilization term.Note that for values k = 2 or 3 of the polynomial degree of practical use, two bisections (m = 2) are enough when η < Ψ.
This refinement may create non-admissible hanging nodes, i.e., hanging nodes with global index larger than Λ.To remove them and guaranteee Λ-admissibility of T * , further refinements should be applied.For the realization of this technical part, we refer to Sect.11.1 in [6].
The following section proves the convergence of the GALERKIN algorithm.
10 Convergence property of GALERKIN  The conditions on γ and δ which lead to the desired estimate are given in (8.6), (10.3) and (10.4).

Conclusions
In this paper, we presented an adaptive VEM of order k ≥ 2 on nonconforming triangular meshes.In the analysis, the space V 0 T of continuous, piecewise polynomials functions of degree k on the triangulation T plays a fundamental role.Indeed, it is contained in the global VEM space, V 0 T ⊆ V T , and guarantees a quasi-orthogonality property for any refinement T * of T , since V 0 T ⊆ V 0 T * .By pivoting on this space, we proved an a posteriori error estimate which does not contain the stabilization term appearing in the VEM discrete formulation.Consequently, we established the convergence of the adaptive VEM algorithm, by a contraction argument.

Figure 3 :
Figure 3: The case k = 3 with 3 refinements of the edge L (in blue) is shown.Red, green and orange lines are the lines needed to refine L the first, the second and the third time respectively.Blue crosses are the degrees of freedom on L of the function living on E. Red squares, green circles, orange diamonds are the degrees of freedom on L generated after the first, the second and the third refinement of L.

Figure 4 :
Figure 4: Blue square are the k + 1 equi-spaced original nodes on the blue edge.Red points represent the nodes added after the refinement of the interval.Black lines show the shapes of the basis ψ i , i = 1, . . ., k, in the case k = 1 (a), k = 2 (b), k = 3 (c).

Figure 5 :
Figure 5: Black points are the proper nodes.Red points represent the hanging nodes generated after a refinement.In (a) the case k = 2 is showed, ζ 1 is the hanging node obtained after the refinement of ξ 1 and ξ 3 and it is the midpoint of ξ 1 and ξ 2 .We notice that if we have called the other red point ζ 2 , ξ 1 and ξ 3 would have been switched.Analogusly, (b) represents the case k = 3.

) 8 . 2
Quasi-orthogonality property Let u T * ∈ V T * be the solution of Problem (3.4) on the refined mesh T * .Hereafter we establish relations between the two energy errors |||u − u T ||| and |||u − u T * |||.The first result follows from Proposition 5.5 and Lemma 3.1; the proof is independent of the used polynomial degree, so we refer to [5, Proposition 5.7].Proposition 8.7 (comparison of the energy error under refinement).For any δ ∈ (0, 1] there exists a constant C E > 0 independent of T and δ such that