A posteriori error estimation and adaptivity in hp virtual elements

An explicit and computable error estimator for the hp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$hp$$\end{document} version of the virtual element method (VEM), together with lower and upper bounds with respect to the exact energy error, is presented. Such error estimator is employed to provide, following the approach of Melenk and Wohlmuth (Adv Comput Math 15(1–4):311–331, 2001), hp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$hp$$\end{document} adaptive mesh refinements for very general polygonal meshes. In addition, a novel VEM hp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$hp$$\end{document} Clément quasi-interpolant, instrumental for the a posteriori error analysis, is introduced. The performances of the adaptive method are validated by a number of numerical experiments.


Introduction
Among the novel technologies developed in the last two decades, polygonal methods [8,14,28,30,34,40,45,47] received an increasing attention, thanks to the advantages that they entail compared to methods based on standard simplicial/quadrilateral meshes.One of the appealing features of polygonal methods is that they dovetail particularly well with adaptive mesh refinements, since polygonal meshes do not require any post-processing within the adaptive procedure.
In the mare magnum of the various theoretical and practical applications of mesh adaptivity, a particular place is occupied by the hp a posteriori analysis in Galerkin methods and in particular in finite element method (FEM) [3,26,31,33,39].Here, one combines mesh refinements with the increase of the dimension of the local approximation spaces.
This results in particularly competitive methods able to capture and to solve in an efficient way the singular behaviour of solutions to partial differential equations (PDEs), but comes at the price of a more involved framework.For instance, one has to face with the problem of deciding whether to refine in h (mesh refinement) or in p (increasing the dimension of local space), has to construct an error estimator with explicit bounds in terms of the distribution of local polynomial degrees, and has to construct spaces possibly tailored for high polynomial order (such that, e.g., the condition number of the stiffness matrix is sufficiently robust in terms of p).
Amidst the polygonal methods, we pinpoint the virtual element method (in short VEM, introduced in [8,10]), which has established itself as one of the most successful technologies in this context and enjoyed an increasing community (a very limited list of papers being [4,6,16,18,21,22,24,32,42,44,49,50]).Differently from other polygonal approaches, the VEM is built without an explicit knowledge of the functions in the approximation spaces.The construction of the method is based on two main ingredients, which are computable only via the degrees of freedom with whom the spaces are endowed: projections into polynomial spaces and particular bilinear forms stabilizing the method.Virtual elements are particularly suitable for the hp framework not only due to the mesh refinement flexibility, but also due to the fact that variable (from element to element) polynomial degrees can be included very naturally in the formulation, without resorting to"ad hoc" modifications as in finite elements.
The a posteriori version of VEM was investigated in [17,20,25,43] but the hp adaptivity has never been targeted before and is the topic of the present work.This paper falls in a series of articles covering different aspects of the p and hp version of VEM, namely a priori error analysis on quasi-uniform meshes [11], approximation of corner singularities [12], multigrid algorithms [5], ill-conditioning in two [36] and three dimensional problems [29], Trefftz and non-conforming approaches [27,37], theoretical and numerical analysis of the stabilization typical of VEM [35].
In the present work, we introduce a VEM computable error estimator for a two dimensional model problem and we prove lower and upper bounds with respect to the energy error that are explicit both in the mesh size and in the distribution of local degrees of accuracy (which are the VEM counterpart of polynomial degrees for FEM).Additionally, we recall from [35,39] some hp inverse inequalities on triangles and polygons.Besides, we introduce and prove error estimates of a novel hp Clément type quasi-interpolant, which is constructed starting from the quasi-interpolant of [38] on subtriangulations of the underlying polygonal meshes.Such estimates are an interesting result on their own, independently of their application in the a posteriori error analysis.
We also present many numerical experiments; more precisely, after verifying the lower and upper bounds of the error estimator, we recall from [39] an hp refinement strategy and we apply it to the proposed adaptive scheme.Such refinement strategy heuristically suggests which are the marked elements where the solutions is supposed to be "smooth" or "singular"; in the former case, p refinement are effectuated (since p VEM converges exponentially when approximating analytic solutions, see [11]), in the latter, h refinement are performed (geometric refinements towards the singularities also give exponential convergence in terms of the number of degrees of freedom, see [12]).
Possible further developments of this work include extension to more complex problems and to the three dimensional case.Moreover, one could also take into account in the adaptive strategy an hp derefinement process including the agglomeration of elements.
The outline of the paper is the following.We begin in Section 2 by presenting the model problem and the hp version of the virtual element method and next, in Section 3, we introduce a set of technical tools needed for the a posteriori error analysis; in particular, we discuss about approximation by functions in the virtual element space, hp polynomial inverse estimates and extension operators from an edge into the interior of a triangle.Successively, in Section 4 we deal with the a posteriori error analysis: we introduce an error estimator and we assert the "reliability" and the "efficiency" of such an estimator.Finally, a number of numerical experiments, including a comparison with the pure h refinement strategy, are presented in Section 5.

Notation
In the remainder of the paper, we will employ the following notation.Given D a measurable open set in R 2 , we denote by L 2 (D) and H s (D), s ∈ N, the standard Lebesgue and Sobolev spaces endowed with the standard inner products and seminorms: The Sobolev norms are denoted by: Sobolev spaces with fractional order are defined via interpolation theory [48].
It is worth to highlight with a separate notation the H 1 inner product over domain D: We refer to [1] for the definition of spaces, inner products and (semi)norms.Given ∈ N, we set P (D) to be the space of polynomials of degree over D. Finally, we write a b and a ≈ b in lieu of: there exist positive constants c 1 , c 2 and c 3 independent of the polynomial degree and of the mesh-size such that a ≤ c 1 b and c 2 a ≤ b ≤ c 3 a, respectively.

Model problem and the hp version of the virtual element method
In this section, we briefly review the model problem and the hp version of the virtual element method.
Given Ω ⊂ R 2 a polygonal domain and f ∈ L 2 (Ω), we consider the Poisson problem with homogeneous Dirichlet boundary conditions and its weak formulation where we have set Next, we introduce a virtual element method based on polygonal meshes and with nonuniform local degree of accuracy for the approximation of problem (2).We note that virtual element methods with nonuniform degree of accuracy were firstly introduced in [12].Given {T n } n∈N a sequence of decomposition of Ω into polygons with straight edges, we associate to each T n , n ∈ N, its set of vertices V n and boundary vertices V b n , its set of edges E n and boundary edges E b n .To each K in T n , we associate its diameter h K , its barycenter x K , its set of vertices V K and its set of edges E K .To each edge s ∈ E b n , we associate once and for all a unit normal n s = n.We require that T n is a conforming polygonal decomposition of Ω for every n ∈ N, i.e. every internal edge belongs to the intersection of the boundary of two polygons.
Besides, we demand the two following geometric assumptions on sequence T n .
(D1) Every K ∈ T n is star-shaped with respect to a ball (see [23]) with radius greater than or equal to γh K , where γ is a universal positive constant.
(D2) Given K ∈ T n , each of his edges s ∈ E K has length greater than or equal to γh K , where γ is a universal positive constant.
A consequence of assumptions (D1) and (D2) which will be extensively used in the following is highlighted in Remark 1.
Remark 1. Owing to assumptions (D1) and (D2), the following fact holds true.The subtriangulation T n = T n (K) of K, obtained by joining the vertices of K to the center of the ball with respect to which K is star-shaped, is made of triangles that are star-shaped with respect to balls with radius greater than or equal to γh K , where γ is a universal positive constant depending only on γ and γ.For a proof of this fact, see [35,Chapter 2].In the forthcoming analysis, we will often make usage of standard functional inequalities (Poincaré, trace, . . .).If not explicitly mentioned, the constants in such inequalities, and therefore also in all the results we will present, depend solely on the shape regularity of T n and T n .
Next, we associate to each element K ∈ T n a local degree of accuracy p K .To each boundary edge s ∈ E b n , we associate p s equal to p K , where K is the unique polygon having s as an edge.On the other hand, to each internal edge s ∈ E n ⊂ E b n we associate p s equal to max(p K1 , p K2 ), where Given now T n the regular subtriangulation of Ω obtained by gluing the local triangulations introduced in Remark 1, we let V n and E n be the sets of vertices and edges of T n , respectively.To each T ∈ T n with T ⊂ K, we associate p T equal to p K .We denote by p the vector of degrees of accuracy associated with polygonal mesh T n , whereas we denote by p the vector of degrees associated with triangular subdecomposition T n .To each edge s of T ⊂ K we associate p s equal to p s , where s = T ∩ K = s.To the other edges s of T , we associate p s equal to p K .We also associate with each vertex (3) See Figure 1 for a graphical idea regarding the distribution of the degrees over T n and T n .We also ask for the following assumption on the local degree of accuracy distribution associated with T n : (P1) for every K and K in T n , there exists a positive universal constant c such that: As a consequence of the assumption (P1) and the above construction, it also holds p s ≈ p K whenever s is an edge of K.
We introduce now the virtual element method associated with the Poisson problem (2).In particular, following [8], we introduce a finite dimensional subspace V n of V , a discrete bilinear form a n : V n × V n → R and a discrete right-hand side f n , such that the method is well-posed and its output u n approximates the solution u of (2).In Section 2.1, we describe the virtual element space V n .Next, in Section 2.2, we present a computable choice for the discrete bilinear form a n .Finally, in Section 2.3, we discuss the construction of the discrete right-hand side f n , • n .

The virtual element space
The present section is devoted to introduce the virtual element space V n .
We begin by defining the local virtual element spaces.Given K ∈ T n , we set first The local virtual element space over polygon K reads Let us consider the following set of linear functionals defined on V (K): given v n ∈ V n (K), • the point-values of v n at the vertices of K; • for every edge s ∈ E K , the point-values at the p s − 1 internal Gauß-Lobatto nodes on s; • the internal moments where {m α , α ∈ N 2 , |α| ≤ p K − 2} is any basis of P p K −2 (K), provided that m α ∞,K ≈ 1 and that it is invariant with respect to homotetic transformations.
These linear functionals are a set of unisolvent degrees of freedom (dofs), see [8].
Remark 2. The choice of the polynomial basis dual to internal moments play a fundamental role in the prospective ill-conditioning of the high-order virtual element method.We refer to [29,36] for a deep insight on this topic and to Section 5 for an explicit choice of such polynomial basis.
Although the space V n (K) is not defined explicitly inside the element K, by means of the set of degrees of freedom above, it is possible to compute a couple of local projectors, see [10].The first operator is the The second one is an When no confusion occurs we will write Π 0 p K −2 and Π ∇ p K in lieu of Π 0,K p K −2 and Π ∇,K p K .The global virtual element space is obtained as in FEM by a standard conforming coupling of the local sets of degrees of freedom and by enforcing homogeneous boundary conditions on ∂Ω: We also introduce the set of global degrees of freedom and the global canonical basis where δ i,j is the Kronecker delta.
Remark 3. We point out that if we aim to approximate a Poisson problem (1) with a nonhomogeneous Dirichlet boundary condition g ∈ H 1 2 (∂Ω), then the global space ( 10) is modified into , where g GL|s denotes the Gauß-Lobatto interpolant of g over the boundary edge s.In fact, Gauß-Lobatto interpolation on the boundary guarantees optimal approximation properties for the p and hp versions of VEM, see [27].However, in order to avoid cumbersome technicalities in the following, we assume henceforth to deal with homogeneous Dirichlet boundary conditions.

The discrete bilinear form
Here, we introduce the discrete bilinear form a n associated with method (5).To this purpose, we follow the guidelines of [8,11].Given a K n : V n (K) × V n (K) → R local bilinear forms, the (global) discrete bilinear form is a sum of such local contributions: where the local discrete bilinear forms a K n are assumed to satisfy (A1) polynomial consistency: for every K ∈ T n , (A2) local stability: for every K ∈ T n , there exist two positive constants α * (p K ) ≤ α * (p K ) independent of h K but possibly depending on the local degree of accuracy p K , such that The assumption (A1) implies that the discrete bilinear form is locally exact for piecewise polynomial functions having a proper degree.On the other hand, the assumption (A2) guarantees the coercivity and the continuity of the global discrete bilinear form and hence the well-posedness of method (5).
A number of computable discrete bilinear forms can be found in the literature.We refer to [12] for what concerns the first "p-explicit" stabilization available in the VEM literature and to [35,Chapter 6] and to [15] for a survey of the topic.
Following [8], the local bilinear forms a K n can be built with the aid of an auxiliary local stabilizing bilinear form S K : ker(Π ∇ p K ) × ker(Π ∇ p K ) → R. For every K ∈ T n , the bilinear form S K must satisfy where c * (p K ) and c * (p K ) are two positive constants independent of h K but possibly depending on the local degree of accuracy p K .The local discrete bilinear form where the energy projector Π ∇ p K is defined in ( 9), satisfies assumptions (A1) and (A2) with [8].
An explicit choice of the stabilization, and consequently also of the discrete bilinear form, will be discussed in Section 5.

The discrete right-hand side
Here, we present the discrete right-hand side of method (5).In particular, we split the global discrete right-hand side f n , • n as a sum of local contributions: where the local discrete duality pairing are defined as For a deeper analysis concerning the discrete right-hand side, we refer to [2,9,11].

Technical tools
In this section, we introduce and discuss some technical tools that we will employ in the a posteriori error analysis.In Section 3.1, we prove local L 2 and H 1 hp approximation properties in terms of functions in the virtual element space.In Section 3.2, we recall a couple of hp polynomial inverse estimates, whereas, in Section 3.3, we recall from [39] the existence of a lifting operator of a polynomial from any edge of a triangle into its interior with good stability properties.Henceforth, we adopt the following notation concerning spaces of piecewise continuous polynomials over polygonal and triangular meshes.Given T n either T n or T n , we write where the choice of the polynomial degrees p s on the edges is picked as in Section 2.

Local hp approximation estimates
In this section, we discuss about approximation properties of functions in the local virtual element spaces V n (K) defined in (6).In particular, we prove local approximation estimates in the L 2 and H 1 (semi)norms and in the L 2 norm on the boundary.Previous results on hp-VEM approximation can be found in [11,12].We first need to recall a technical tool from [38,39].Given T n the (regular) triangular subdecomposition of Ω, defined locally by the subtriangulations of the polygons of T n introduced in Remark 1, and given a vertex V ∈ V n of subtriangulation T n , we set the triangular patch around vertex V as Given D either a polygon of T n or a triangle of T n , we define the triangular patch around D as Besides, given an edge s ∈ E n , we also define the triangular patch around s as In Figures 2, 3     It can be proven that, using assumptions (D1)-(D2), diam(ω D ) ≈ diam(D), D being either a polygon or a triangle and ω(D) being the associated triangular patch.
We now recall the following polynomial Clément-type approximation result over triangles.
Lemma 3.1.Let the assumptions (D1)-(D2) be valid.Then, given u ∈ H 1 0 (Ω), there exists 14), such that, for each triangle T ∈ T n and for each s edge of E n (i.e. of the skeleton of T n ), the following estimates hold true: where c M is a positive constant independent on u, h T and p T but depending on the shape-regularity of T n and hence of T n .Above, the symbols ω T and ω s represent the triangular patch around triangle T and the triangular patch around edge s defined in ( 16) and (17), respectively.
Proof.The proof, based on the partition of unity method [7], can be found in [38,Theorem 3.3].
Lemma 3.1 is a key ingredient for the following result which asserts the existence of an hp VEM Clément quasi-interpolant.Proposition 3.2.Let the assumptions (D1)-(D2) be valid.Given u the solution to (5), a convex polygon K ∈ T n and s any of its edges, there exists a function u I ∈ V n , such that its restriction to K satisfies the following estimates: where c is a positive constant independent on u, h K and p K but depending on the shape-regularity of T n and hence of T n , and where ω K and ω s follow the notation introduced in ( 16) and (17).
Proof.Given K ∈ T n convex, we define u I ∈ V n (K) as the solution to the following Poisson problem with nonhomogeneous Dirichlet boundary conditions where u M is the quasi-interpolant from [38] introduced in Lemma 3.1.In practice, we fix u I on the boundary to be the polynomial quasi-interpolant guaranteeing hp approximation properties on triangles and we fix the internal moments of u I by demanding Roughly speaking, this last condition enforces the (distributional) identity ∆u I = ∆u M in an approximated sense.
Next, we investigate the bound on the H 1 seminorm in (19a).Given any bubble function w in V n (K) ∩ H 1 0 (K), we observe that, owing to (20), Therefore, since any function w in the virtual element space with w |∂K = u M |∂K = u I|∂K can be written as u I + w, one easily gets As a consequence, where we choose u I ∈ V n (K) satisfying the following local problem: Plugging (23) in (22), we get whence, applying (18a) and [11,Lemma 4.2] to the first and second terms on the right-hand side of (24), respectively, We underline that the constant c depends only on the shape-regularity of T n and hence of T n .
Finally, we investigate the bound on the L 2 norm in (19a).To this purpose, we employ the hypothesis that K is convex and we consider the auxiliary problem The following standard a priori bound on the solution to the auxiliary problem (25) for convex K is valid: Owing to the fact that u M − u I | ∂K = 0 and to the orthogonality of u M − u I with respect to virtual bubble functions (21), we have where ψ II is a function in V n (K) ∩ H 1 0 (K) providing optimal hp approximation estimates of ψ, see e.g.[11,Lemma 4.3], being c a positive constant depending only on the shape-regularity of K.
Recalling (26), we obtain The assertion follows from a triangle inequality and the bound on the H 1 seminorm in (19a).
for all ε > 0 arbitrarily small, where α K denotes the largest interior angle of K.

Inverse estimates on triangles and polygons
In this section, we recall two hp polynomial inverse inequalities that will be instrumental in the forthcoming analysis.We begin with the following 1D result.
given ψ the quadratic bubble function over I and given 0 ≤ α 1 ≤ α 2 , there exists a positive constant such that, for all q ∈ P p ( I), p ∈ N, Proof.See [19, Lemma 4].
Given K ∈ T n , we define a piecewise bubble function ψ K over T n (K) as follows: , where ψ T is the cubic bubble function on triangle T for all T ∈ T n (K), (29) where we recall that T n (K) is the subtriangulation of K introduced in Remark 1.
The following hp polynomial inverse inequality over a polygon holds true.
Lemma 3.4.Given K a polygon in T n , let ψ K be the piecewise bubble function associated with the subtriangulation T n (K) of K defined in (29).For all −1 < α ≤ β, there exist a positive constants c depending only on α, β, and the shape-regularity of the subtriangulation T n (K), such that, for all q ∈ P p K (K),

A lifting operator
In this section, we recall the existence of a lifting operator from an edge of a triangle to its interior.
Lemma 3.5.Given T ∈ T n a triangle and s any of its edges and given 1/2 < α ≤ 2, let ψ s be the quadratic bubble function associated with edge s.Then, for every q ∈ P ps (s) and for every ε > 0, there exist a lifting E(q) ∈ H 1 (T ), such that where c α is a positive constant depending only on α and on the shape of T .
Proof.The proof follows from [39, Lemma 2.6] and a scaling argument.
Note that, since E(q) vanishes on ∂T \ s, then it can be extended to 0 on the remaining part of the polygon containing T as part of its subtriangulation.As a consequence, ( 32) and (33) can be "generalized" in a straightforward way, employing on the left-hand side (semi)norms on the complete polygon in lieu of their counterparts on triangle T .

A posteriori error analysis
In this section, we build an error estimator and we prove that it can be upper and lower bounded by the H 1 error of the method, with constants explicit in terms of the discretization parameters.
The remainder of the the section is structured as follows.In Section 4.1, we write the residual equation, whereas, in Section 4.2, we construct an error estimator and we show that the energy error can be bounded in terms of such an estimator with an explicit dependence in terms of the discretization parameters.Next, in Section 4.3, we show that the local estimator can be bounded by the H 1 seminorm of the quasi-local error plus a couple of terms involving the oscillation of the right-hand side of the problem (1) and the stabilization of the method.Finally, in Section 4.4, we summarize the foregoing results.
Henceforth, we assume for the sake of clarity that all the polygons in T n are convex.The nonconvex case is discussed in Remark 5. Furthermore, we also assume in the forthcoming analysis that p K ≥ 2 for all K ∈ T n ; this allows us to avoid a cumbersome notation when treating the discrete right-hand side, cf.Section 2.3.Thanks to this assumption, we will write f n , v n n as (f n , v n ) 0,Ω , where f n| K = Π 0 p K −2 f for all K ∈ T n .Finally, we set the jump across an edge s ∈ E n as where v + and v − are the restriction (assumed sufficiently regular) of v over ∂K + and ∂K − , respectively, and where s ⊆ ∂K + ∩ ∂K − .Note that we are assuming that the unit normal n is pointing outside K + and inside K − .

The residual equation
We begin by introducing the residual equation.Let e := u − u n be the difference between the solution to the continuous (2) and discrete (5) problems, respectively.One has, for any v ∈ H 1 0 (Ω) and We rewrite the last term on the right-hand side of (35) by observing that, for all w ∈ H 1 0 (Ω), the following holds true: where we recall that the jump • s is defined in (34) and where we recall that we are assuming that the unit normal n of each edge is fixed once and for all.We highlight the presence of a polynomial projector in (36); without such a projector we would not be able to compute exactly the jump across the edges of normal derivatives of functions in the virtual element space, which we anticipate will be part of the error residual of the method.
Plugging (36) in (35) The first two and the last term on the right-hand side of (37) are analogous to the terms appearing in the FEM counterpart of the residual equation, see [39,Lemma 3.1].The only difference here, is that we consider the H 1 projection of the discrete solution.Note that in the FEM framework such a projection applied to the solution coincides with the solution itself (since the solution is a piecewise polynomial).
The remaining terms (virtual element consistency terms) on the right-hand side of (37) are instead typical of the virtual element setting and they take into account the fact that the discrete bilinear form is only an approximation to the exact one.

Error estimator and upper bounds
In this section, we introduce a computable error estimator and we prove the lower and upper bounds of the energy error in terms of such estimator.
To this purpose, we use the residual equation ( 37) with v = e := u−u n and χ n = e I = (u−u n ) I , where we recall that the local hp approximation properties of e I are described in Proposition 3.2.We obtain We estimate the five local terms separately.We start with the term I. Applying Proposition 3.2, which is valid since we are assuming K convex, we get where we recall that the hidden constant depends solely on the shape-regularity of T n and hence of T n , see Remark 1.In the following, when no confusion occurs and when unnecessary, we will omit to highlight such dependence.Secondly, we investigate II, the term involving the oscillation of the right-hand side.Owing to L 2 orthogonality and hp approximation properties of such projector, see e.g.[11,Lemma 4.2], we write where we recall that the L 2 projector Π 0 p K −2 is defined in (8).
Next, we deal with V , the projected jump edge residual term: where in the last but one inequality we employed (19b).The first virtual element consistency term IV can be bounded instead using (19a) (in the last but one inequality) and the stability property (12) (in the last inequality): We emphasize that we make appear the stabilization term since in the definition of the error residual we want to have computable quantities only.
Finally, we study III, the second virtual element consistency term : whence, recalling Proposition 3.2 and the stability property (12), In order to simplify the notation, we write henceforth Collecting the estimates on the five terms on the right-hand side of (38), we get, after some trivial algebra, where we have set the local error estimators as For ease of notation, we define next η p K , which takes into account both the bulk and the edge error estimators η K and η s , s ∈ E K , defined in (41), as We highlight two facts.The first one is that, apart from the term involving the stabilization, the upper bound is analogous to the one of hp FEM, see [39,Lemma 3.1].The second observation is that for nonconvex K, we would have an additional factor p for all ε > 0 arbitrarily small in front of η 2 K in (40), where α K denotes the largest interior angle of K, see Remark 4.

Lower bound
In this section, we bound the local error estimator η p K introduced in ( 42) with the quasi-local H 1 seminorm of the error of the method plus a term related to the oscillation of the right-hand side and a term related to the stabilization of the method.After recalling a technical tool in Section 4.3.1,we prove in Sections 4.3.2 and 4.3.3 the bounds on the internal and edge residuals, respectively, in terms of the local energy errors.

An inverse inequality
Here, we recall an hp polynomial inverse inequality, which will be used in Section 4.3.2.We stress that the triangular counterpart of this result was shown in [39,Lemma 3.4].
Lemma 4.1.Given a polygon K ∈ T n , we set ψ K the piecewise bubble function associated with polygon K ∈ T n as in (29).Then, for all q ∈ P p K (K) and 1/2 < α ≤ 2, the following holds true: where the hidden constant depends on the shape-regularity of T n and hence of T n .

Bounding the internal residual
In this section, we bound the local internal residual R K appearing on the right-hand side of (42).
To this purpose, we note that plugging χ n = 0 in (37), we have Let us focus our attention on a single element K. Given ψ K the piecewise bubble function associated with element K defined as in ( 29), we extend it to 0 outside K.We then choose From ( 45), we get As a consequence, applying the stability bounds (12), recalling that R K ∈ P p K (K) and applying the hp polynomial inverse estimate on polygons (43), one has, for all 1/2 < α ≤ 2, or, equivalently, We are now ready to prove the bound on the internal residual.Recalling that 1/2 < α ≤ 2, the hp polynomial inverse estimate on polygons (30) implies We pick α = 1/2 + ε, with ε > 0 arbitrarily small, and get

Bounding the edge residual
Next, we bound the edge residual R s appearing on the right-hand side of (42).Without loss of generality, s is an internal edge; the general case follows analogously.We henceforth consider a function R s given by E(R s ), where the lifiting operator E is defined in Lemma 3.5.We recall that the restriction of R s on s is equal to ψ α s R s , being ψ s defined as the quadratic edge bubble function on s and where 1/2 < α ≤ 2.
In the following, we assume that R s can be extended to 0 outside T 1 ∪ T 2 , see (17).Let us denote by K 1 and K 2 the two polygons containing the triangles T 1 and T 2 .
We substitute v = R s and χ n = 0 in (37), obtaining We observe that the following bound ot the third term on the right-hand side of (47) holds true: We deduce from (47) that and, combining this with (48), we arrive at Recalling that p s ≈ p Ki for all edges s ⊂ ∂K i , i = 1, 2, see (4), and applying Lemma 3.5 with 1/2 < α ≤ 2 on |R s | 1,Ki and R s 0,Ki , we obtain, for every ε > 0, Therefore, we can write Hence, squaring both sides, one deduces Using that h s ≤ h Ki for i = 1, 2, also recalling ( 4) and ( 46), we get , whence, using again (4), Selecting ε = p −2 s and using once more (4), one obtains So far, we have assumed that 1/2 < α ≤ 2. In order to get the desired bound on the edge residual, i.e. the one with α = 0, we apply Lemma 3.3 with α 1 = 0 and α 2 = 1/2 + ε, getting Finally, note that for any value of p Ki ∈ N it holds that p 2 p 2 K i Ki ≤ 3 and thus such term can be discarded.

Conclusions
We collect here the lower and upper bounds discussed in the foregoing Sections 4.2 and 4.3.Note that such bounds are explicit in h and p. Theorem 4.2.Assume that the assumptions (D1)-(D2)-(P1) hold true.Let u and u n be the solutions to (2) and (5), respectively, and let e = u − u n .For all K ∈ T n , let η p K be the error residual defined in (42) and let ρ p K and ζ p K be defined in (41).Then, assuming that all the polygons K ∈ T n are convex, the following global upper bound holds true: Further, for every K ∈ T n and for all ε > 0, the following local lower bounds hold true: where The hidden constants in (49) and (50) depend solely on ε and on the shape-regularity of T n and hence of T n , see Remark 1.
On the light of Theorem 4.2, the global (computable) error estimator that we propose is Note that in the bounds ( 49) and ( 50) there appear additional multiplicative terms depending on max(α −1 * (p K ), α * (p K )).On the other hand, such terms were numerically shown in [12] to have a very mild dependence in terms of p K .We refer to Remark 8 for additional comments regarding the effects of the "stability" terms α −1 * (p K ) and α * (p K ).We also underline that on the right-hand side of (50), in addition to the energy error, other two terms appear.The term ρ p K is related to the oscillation of f , the datum of the problem (1), and is typical also in the finite element framework.On the other hand, the term ζ p K deals with the nonexactness of the discrete bilinear form and is standard in a posteriori error analysis of VEM, see [17,25].
Remark 5.In presence of nonconvex polygons, the bound (49) needs to be modified; in particular, it appears in front of η 2 p K a suboptimal factor p −ε for all ε > 0 arbitrarily small, α K being the largest interior angle of K.

Numerical results
In this section, we present a set of numerical experiments investigating on the performances of the error estimator introduced in (51).More precisely, we validate in Section 5.1 the estimates presented in Theorem 4.2, whereas, in Section 5.2, we recall from [39] an hp refinement algorithm which permits us to apply the adaptive hp virtual element method on a number of test cases.Before presenting the results, we have to settle some features of the method, which so far were kept at a very general level.First of all, we underline that in the virtual element setting, it is not possible to compute explicitly the exact H 1 error (since functions in virtual element spaces are not known in closed-form) and therefore we compute instead the following quantity: where (9).It is possible to show that the two quantities converge to zero with the same h and p rate.
Another important aspect of the method is the choice of the stabilization in (13).We adopt here the so called "D-recipe", which was firstly introduced in [13] and whose performances were investigated in deep in [29,36].If we denote by S K the matrix representing the stabilization S K on element K with respect to the canonical basis (11), then we set where δ i,j denotes the Kronecker delta.Among the possible stabilizations available in the virtual element literature, the one defined in (53) is one of the most appealing in terms of robustness of the method, both when considering high degrees of accuracy and in presence of badly-shaped elements.
Finally, we address the issue of picking a "clever" polynomial basis dual to the internal moments (7).In particular, following again [29,36], we employ a basis {m α } p−2 |α|=0 which is L 2 (K) orthonormal for all K ∈ T n .Such a basis can be built, for instance, by orthonormalizing a basis of scaled and shifted monomials.Picking an orthonormal basis dual to internal moments is a choice particularly suited for the p version of the method since it allows to effectively damp the condition number of the stiffness matrix for high values of the polynomial degree.
In the numerical experiments, both the error estimator and the computable error (52) are normalized by |u| 1,Ω .

Performances in terms of p of the error estimator
In this section, we investigate the performances of the error estimator for uniform p refinements; that is, we fix a coarse mesh and we achieve convergence by raising the polynomial degree in all mesh elements (note that the performances of h refinements where the topic of [25, Section 6.1] and therefore are not explored here).To this end, we consider three test cases with known exact solution As usual, the loading term and the (Dirichlet) boundary conditions are set in accordance with the exact solution.Function u 1 is analytic, function u 2 belongs to H 3−ε (Ω 1 ) for all ε > 0 and is the "natural" singular solution to an elliptic problem over square Ω 1 , function u 3 belongs to H 5 3 −ε (Ω 2 ) for all ε > 0 and is the "natural" singular solution to an elliptic problem over the L-shaped domain Ω 2 .
We test the performances of the error indicator employing two types of mesh, namely a cartesian and a Voronoi mesh.In Figure 6, we depict the two meshes (on domain Ω 1 ), their counterparts in the L-shaped domain Ω 2 being analogous.In Figures 7 and 8, we plot (for different values of the degree p) the computable error (52) versus the computed error estimator (51) for the three exact solutions in (54) on the two meshes of Figure 6.
We observe that, when approximating analytic solutions, the error estimator and the computed error behave practically in the same way.For singular functions, the behaviour is instead slightly different and in particular the error estimator is slightly smaller than the computed error.This is due to the small pollution effect of the stabilization, see ( 49)- (50); however, if the solution is analytic, the pollution factors are negligible, since all the local error estimators, as well as the exact error, decrease exponentially in terms of p, see [11].

The hp adaptive refinement strategy
In this section, we present an hp refinement strategy, based on the standard procedure Before discussing the hp refinement strategy, we describe how to perform h refinements.The standard strategy that one can follow consists in subdividing a polygon K with N K edges into a number of quadrilaterals smaller than or equal to N K obtained by connecting the barycenter of K to the midpoints of its edges.Note that, whenever two edges belong to the same straight line, the procedure is slighly modified: in such case we connect the barycenter with the resulting hanging node, see Figure 9.
We stress that this procedure can be performed only on elements containing their barycenter (e.g.convex elements), but, importantly, starting from a convex element this construction generates only convex elements.Besides, the hanging nodes popping-up with this strategy fit naturally in the polygonal framework of VEM.Note that, when dealing with cartesian meshes, the refinement strategy of Figure 9 reduces to standard square refinement.
We now briefly review the hp refinement strategy of [39,Section 4.2], here adapted to the present method.The basic idea behind this approach is that, on the elements on which the exact solution to the continuous problem is "regular", one performs p refinement (since the p version of the method leads to exponential convergence of the error in terms of p when approximating analytic solutions, see [11]), whereas, on the elements where the solution has a "singular" behaviour, one performs h refinement (since geometric mesh refinements lead to exponential convergence of the error in terms of the number of degrees of freedom when approximating singular functions, cf.[12]).
The hp refinement strategy reads as follows.At the n-th level of the algorithm, one marks the elements on which the error estimator is sufficiently large.To this end, one marks all the elements such that η 2 comp,K,n ≥ ση 2 a for some σ ∈ (0, 1), where being η 2 comp,K,n and η comp,n computed as in (51) for all K ∈ T n and for all n ∈ N. Let us now be given, at the n-th adaptive refinement step, a "prediction-estimator" η pred,K,n on each element K of the mesh T n .Such estimator is instrumental in order to heuristically decide whether (elementwise) the exact solution is regular or not.
If one performs on element K an h refinement, then the expectations are, under the hypothesis of analyticity of u, that the error reduces by a factor decreasing exponentially in terms of p, see [11,12].Therefore, given N K the number of son elements of K, one sets the prediction-estimator η pred,Ki,n on K i to: for some γ h fixed positive parameter.Note that the term 0.5 in the above equation stems from the fact that we are roughly dividing the element diameter h K by two when performing the procedure in Figure 9. Instead, if one performs on element K a p refinement, one expects, under the hypothesis of analyticity of u, a reduction of the error by a fixed factor γ p .Thus, one sets the prediction-estimator η pred,p K +1 on K to: for some 0 < γ p < 1 fixed parameter.
On the elements that are not marked for refinement, one updates in a trivial fashion as for some γ n fixed positive parameter.
So far, we have constructed, at the n-th level of the hp refinement process, a set of predictionestimators for the n + 1-th level; as already underlined, such prediction-estimators have the scope of "guessing" elementwise the behaviour of the exact solution.
What one does at this point is that on all marked elements he checks whether If this is the case, the actual error estimator is "larger" than the predicted one, where we have assumed that the solution was analytic; this means that the solution is not sufficiently "regular" on the element and therefore an h refinement has to be performed.The p refinement is effectuated otherwise.
Since in the first refinement step the prediction-estimators are not available (since there is not a 0 level in the adaptive construction of the spaces), the hp refinement boils down to a simple h refinement.To this purpose, in the first iteration of the adaptive algorithm, it suffices to set on all elements K.
In Algorithm 1, we present the hp refinement process discussed so far.
given fixed positive parameters σ, γ h , γ p , and γ n : end if end for Remark 6.It may occur within the adaptive algorithm that the assumptions (D2), guaranteeing the presence of edges with size comparable to that of the element to whom they belong, and (P1), guaranteeing comparable degrees of accuracy on neighbouring elements, are not valid.However, from the forthcoming numerical experiments, it is evident that the method is robust in this respect.
Remark 7. The approach of Algorithm 1 is not the only one available in the framework of hp Galerkin methods.We refer to [41] for an overview of different hp refinement approaches.
In our experiments, we set the parameters in Algorithm 1 to [39] σ = 0.75, γ h = N K , γ p = 0.4, γ n = 1, where we recall that N K is the number of elements into which polygon K is split in case of h refinement.
Moreover, we test the method on the test case with known solution u 3 , see (54), and where we recall that Ω 1 is the square domain introduced in (54).We underline that the solution u 4 is analytic but has a steep derivatives around (0.5, 0.5).
Remark 8. Algorithm 1 is based on h refinements where the solution is assumed to have a "singular" behaviour, whereas it is based on p refinements where the solution is assumed to be "smooth".Therefore, the pollution factor max(α −1 * (p), α * (p)) appearing in the estimates of Theorem 4.2, behaves in the following heuristic fashion.It does not blow up when doing h refinements since it is independent of h.Instead, when doing p refinements, it multiplies the factor ζ p K defined in (41) which, since we are assuming that the solution is analytic in case of p refinements, is converging exponentially, in terms of p. Thus, the pollution factor, which typically grows algebraically in terms of p, see [5, Lemma 2.5] and [12,Theorem 2], is not in principle spoiling the behaviour of the adaptive scheme.This observation is confirmed by the numerical experiments below.
In Figures 10 and 11, we show the hp mesh refinement after 3 and 10 steps of Algorithm 1 applied to the two test cases with known solution u 3 and u 4 , starting from a coarse cartesian mesh.
Furthermore, we depict in Figures 12 and 13 the same set of experiments starting from a coarse Voronoi mesh.In both cases, the h refinements are performed as in Figure 9.
From Figures 10 and 12, it is possible to appreciate the fact that Algorithm 1 is actually catching the singular behaviours of the solution u 3 .In particular, the mesh is geometrically refined towards such singularity.Note that, taking into account that the degrees of accuracy increase where the solution is smooth, it may seem that p refinement are effectuated also towards the re-entrant corner.This is not what actually happens, since lots of geometric refinements are performed and high degree of accuracy is picked only in outer layers.
Regarding Figures 11 and 13, we can observe that the algorithm indeed detects that the target function is regular and therefore mainly applies p-refinements; such refinements are slightly more concentrated in the area with higher derivatives.In order to have a more quantitative evaluation on the algorithm performance, we plot in Figure 14 the error (52) and the error estimator (51) in terms of the square root of the total number of degrees of freedom, for the case with exact solution u 4 .We can appreciate the exponential convergence of the method, which resembles that of the p version of VEM when approximating analytic solutions [11].
Finally, we concentrate on the singular solution case u 3 .We aim to investigate whether Algorithm 1 leads to a decay of the error which is exponential in terms of the cubic root of the overall degrees of freedom also in this more complex situation.This is indeed what one expects from the a priori error analysis of hp FEM [46] and of hp VEM [12], employing graded meshes towards the singularity and increasing elementwise the degrees of accuracy.In Figure 15, we plot the computable error (52) and the error estimator (51).We observe that, after some adaptive refinement steps, the decay gets the desired exponential slope.
As a final experiment, we compare in Figure 16 the performances of the hp adaptive refinement strategy of Algorithm 1, with a pure h adaptive refinement strategy, which is equivalent to the scheme presented in Algorithm 1 setting η 2 pred,K,n = 0 for all K ∈ T n and for all n ∈ N. From Figure 16, it is evident that the hp adaptive strategy overperformes the h refinement strategy, since, on the portions of the domain where the solution is sufficiently smooth, the p version is able to capture with very few degrees of freedom the same accuracy of an h version with many more degrees of freedom.It is however worth to underline that this comes at the price of having, for high p, a more densily populated stiffness matrix.

Figure 1 :
Figure 1: Local degrees of accuracy distribution.Left: distribution over Tn, the polygonal decomposition, and Vn of Tn.Center-left: distribution over En, the set of edges of Tn.Center-right: distribution over Tn, the subtriangulation of Tn, and Vn, the set of vertices of Tn.Right: distribution over En, the set of edges of Tn.
, 4 and 5 we provide graphical examples of patches around a vertex, a polygon, a triangle and an edge.

Figure 2 :
Figure2: Left: a polygonal mesh Tn and its subtriangulation Tn, in a red circle a vertex V. Right: in light blue, ω V , the triangular patch around vertex V, defined in(15).

Figure 3 :
Figure3: Left: a polygonal mesh Tn, in light red a polygon K. Center: the regular subtriangulation Tn.Right: in light blue, ω K , the triangular patch around K, defined in(16).

Figure 4 :
Figure4: Left: a polygonal mesh Tn.Center: the regular subtriangulation Tn, in light red a triangle T .Right: in light blue, ω T , the triangular patch around T , defined in(16).

Figure 5 :
Figure5: Left: a polygonal mesh Tn, in red an edge s.Center: the regular subtriangulation Tn.Right: in light blue, ωs, the triangular patch around edge s, defined in(17).

Figure 9 :
Figure9: Standard h refinement strategy.One connects the barycenter of the element with the midpoints of every edge.As an exception, whenever two edges belong to the same straight line, one connects the barycenter with the resulting hanging node to the straight line.The red dots denote the vertices; the original polygon is denoted by solid lines.

Figure 10 :
Figure 10: Steps 3 and 10 of the hp refinement algorithm applied to the test case with known solution u 3 introduced in (54).Starting mesh: a coarse Cartesian mesh.Refinement strategy as in Figure 9.

Figure 11 :Figure 13 :
Figure 11: Steps 3 and 10 of the hp refinement algorithm applied to the test case with known solution u 4 introduced in (56).Starting mesh: a coarse Cartesian mesh.Refinement strategy as in Figure 9.

Figure 14 :
Figure 14: Error decay of the computable error (52) and of the error estimator (51) in terms of the square root of the degrees of freedom employing the hp adaptive refinement strategy in Algorithm 1. Exact solution: u 4 introduced in (56) on cartesian (left) and Voronoi (right) meshes.

Figure 15 :
Figure 15: Error decay of the computable error (52) and of the error estimator (51) in terms of the cubic root of the degrees of freedom employing the hp adaptive refinement strategy in Algorithm 1. Exact solution: u 3 introduced in (54) on cartesian (left) and Voronoi (right) meshes.

Figure 16 :
Figure 16: Comparison between the hp and h adaptive refinement in Algorithm 1. Exact solution: u 3 introduced in (54) on cartesian (left) and Voronoi (right) meshes.On the x−axis, we depict the cubic root of the number of degrees of freedom.