A priori and a posteriori error analysis of the lowest-order NCVEM for second-order linear indefinite elliptic problems

The nonconforming virtual element method (NCVEM) for the approximation of the weak solution to a general linear second-order non-selfadjoint indefinite elliptic PDE in a polygonal domain Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} is analyzed under reduced elliptic regularity. The main tool in the a priori error analysis is the connection between the nonconforming virtual element space and the Sobolev space H01(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1_0(\Omega )$$\end{document} by a right-inverse J of the interpolation operator Ih\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_h$$\end{document}. The stability of the discrete solution allows for the proof of existence of a unique discrete solution, of a discrete inf-sup estimate and, consequently, for optimal error estimates in the H1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1$$\end{document} and L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} norms. The explicit residual-based a posteriori error estimate for the NCVEM is reliable and efficient up to the oscillation terms. Numerical experiments on different types of polygonal meshes illustrate the robustness of an error estimator and support the improved convergence rate of an adaptive mesh-refinement in comparison to the uniform mesh-refinement.


Introduction
The nonconforming virtual element method approximates the weak solution u ∈ H 1 0 (Ω) to the second-order linear elliptic boundary value problem for a given f ∈ L 2 (Ω) in a bounded polygonal Lipschitz domain Ω ⊂ R 2 subject to homogeneous Dirichlet boundary conditions.

General introduction
The virtual element method (VEM) introduced in [4] is one of the well-received polygonal methods for approximating the solutions to partial differential equations (PDEs) in the continuation of the mimetic finite difference method [7].This method is becoming increasingly popular [1,3,5,6,16,17] for its ability to deal with fairly general polygonal/polyhedral meshes.On the account of its versatility in shape of polygonal domains, the local finite-dimensional space (the space of shape functions) comprises non-polynomial functions.The novelty of this approach lies in the fact that it does not demand for the explicit construction of non-polynomial functions and the knowledge of degrees of freedom along with suitable projections onto polynomials is sufficient to implement the method.
Recently, Beirão da Veiga et al. discuss a conforming VEM for the indefinite problem (1.1) in [6].Cangiani et al. [17] develop a nonconforming VEM under the additional condition which makes the bilinear form coercive and significantly simplifies the analysis.The two papers [6,17] prove a priori error estimates for a solution u ∈ H 2 (Ω) ∩ H 1 0 (Ω) in a convex domain Ω.The a priori error analysis for the nonconforming VEM in [17] can be extended to the case when the exact solution u ∈ H 1+σ (Ω) ∩ H 1 0 (Ω) with σ > 1/2 as it is based on traces.This paper shows it for all σ > 0 and circumvents any trace inequality.Huang et al. [31] discuss a priori error analysis of the nonconforming VEM applied to Poisson and Biharmonic problems for σ > 0. An a posteriori error estimate in [16] explores the conforming VEM for (1.1) under the assumption (1.2).There are a few contributions [9,16,34] on residual-based a posteriori error control for the conforming VEM.This paper presents a priori and a posteriori error estimates for the nonconforming VEM without (1.2),but under the assumption that the Fredholm operator L is injective.
Since the bounded linear operator L is a Fredholm operator [30, p. 321], (A3) implies that L is bijective with bounded inverse L −1 : H −1 (Ω) → H 1 0 (Ω).The Fredholm theory also entails the existence of a unique solution to the adjoint problem, that is, for every g ∈ L 2 (Ω), there exists a unique solution Φ ∈ H 1 0 (Ω) to The bounded polygonal Lipschitz domain Ω, the homogeneous Dirichlet boundary conditions, and (A1)-(A2) lead to some 0 < σ ≤ 1 and positive constants C reg and C * reg (depending only on σ, Ω and coefficients of L) such that, for any f, g ∈ L 2 (Ω), the unique solution u to (1.1) and the unique solution Φ to (1.4) belong to H 1+σ (Ω) ∩ H 1 0 (Ω) and satisfy (The restriction σ ≤ 1 is for convenience owing to the limitation to first-order convergence of the scheme.)

Weak formulation
(with piecewise versions a pw , b pw , c pw and B pw for ∇ replaced by the piecewise gradient ∇ pw and local contributions a P , b P , c P defined in Subsection 3.1 throughout this paper).The weak formulation of the problem (1.1) seeks u ∈ V such that B(u, v) = (f, v) for all v ∈ V. (1.8) Assumptions (A1)-(A3) imply that the bilinear form B(•, •) is continuous and satisfies an inf-sup condition [11] 0 < β 0 := inf . (1.9)

Main results and outline
Section 2 introduces the VEM and guides the reader to the first-order nonconforming VEM on polygonal meshes.It explains the continuity of the interpolation operator and related error estimates in detail.Section 3 starts with the discrete bilinear forms and their properties, followed by some preliminary estimates for the consistency error and the nonconformity error.The nonconformity error uses a new conforming companion operator resulting in the well-posedness of the discrete problem for sufficiently fine meshes.Section 4 proves the discrete inf-sup estimate and optimal a priori error estimates.Section 5 discusses both reliability and efficiency of an explicit residual-based a posteriori error estimator.Numerical experiments in Section 6 for three computational benchmarks illustrate the performance of an error estimator and show the improved convergence rate in adaptive mesh-refinement.
(M1) Admissibility.Any two distinct polygonal domains P and P in T ∈ T are disjoint or share a finite number of edges or vertices.to every point of a ball of radius greater than equal to ρh P and every edge E of P has a length |E| greater than equal to ρh P .
Here and throughout this paper, h T | P := h P denotes the piecewise constant mesh-size and T(δ) := {T ∈ T : h max ≤ δ ≤ 1} with the maximum diameter h max of the polygonal domains in T denotes the subclass of partitions of Ω into polygonal domains of maximal mesh-size ≤ δ.
Let |P | denote the area of polygonal domain P and |E| denote the length of an edge E. With a fixed orientation to a polygonal domain P , assign the outer unit normal n P along the boundary ∂P and n E := n P | E for an edge E of P .Let E (resp.E) denote the set of edges E of T (resp. of T ) and E(P ) denote the set of edges of polygonal domain P ∈ T .For a polygonal domain P , define mid(P ) := 1 |P | P x dx and mid(∂P ) := 1 |∂P | ∂P x ds.
Remark 1 (consequence of mesh regularity assumption).There exists an interior node c in the sub-triangulation T (P ) := {T (E) = conv(c, E) : E ∈ E(P )} of a polygonal domain P with h T (E) ≤ h P ≤ C sr h T (E) as illustrated in Figure 2.2.Each polygonal domain P can be divided into triangles so that the resulting sub-triangulation T | P := T (P ) of T is shape-regular.The minimum angle in the sub-triangulation solely depends on ρ [13, Sec.2.1].

Lemma 2.1 (Poincaré-Friedrichs inequality).
There exists a positive constant C PF , that depends solely on ρ, such that holds for any f ∈ H 1 (P ) with j∈J E(j) f ds = 0 for a nonempty subset J ⊆ {1, . . ., m} of indices in the notation Some comments on C PF for anisotropic meshes are in order before the proof gives an explicit expression for C PF .
Example 2.1.Consider a rectangle P with a large aspect ratio divided into four congruent sub-triangles all with vertex c = mid(P ).Then, m = 4 and the quotient of the maximal area divided by the minimal area of a triangle in the criss-cross triangulation T (P ) is one.Hence C PF ≤ 1.4231 (from the proof below) is independent of the aspect ratio of P .
Proof of Lemma 2.1.The case J = {1, . . ., m} with f ∈ H 1 (P ) and ∂P f ds = 0 is well-known cf.e.g.[13, Sec.2.1.5],and follows from the Bramble-Hilbert lemma [14,Lemma 4.3.8]and the trace inequality [13,Sec. 2.1.1].The remaining part of the proof shows the inequality (2.1) for the case J ⊆ {1, . . ., m}.The polygonal domain P and its triangulation T (P ) from Figure 2.2 has the center c and the nodes z 1 , . . ., z m for the m := |E(P )| = | T (P )| edges E(1), . . ., E(m) and the triangles T (1), . . ., T (m) with T (j) = T (E(j)) = conv{c, E(j)} = conv{c, z j , z j+1 } for j = 1, . . ., m.Here and throughout this proof, all indices are understood modulo m, e.g., z 0 = z m .The proof uses the trace identity for f ∈ H 1 (P ) as in the lemma.This follows from an integration by parts and the observation that (x−c)•n F = 0 on F ∈ E(T (j))\E(j) and the height (x−c)•n E(j) = 2|T (j)| |E(j) of the edge E(j) in the triangle T (j), for x ∈ E(j); cf.[24,Lemma 2.1] or [25,Lemma 2.6] for the remaining details.Another version of the trace identity (2.2) concerns conv{z j , c} =: F (j) = ∂T (j − 1) ∩ ∂T (j) and reads in T (j − 1) and T (j).The three trace identities in (2.2)-(2.3)are rewritten with the following abbreviations, for j = 1, . . .m, Let t min = min T ∈ T (P ) |T | and t max = max T ∈ T (P ) |T | abbreviate the minimal and maximal area of a triangle in T (P ) and let Π 0 f ∈ P 0 ( T (P )) denote the piecewise integral means of f with respect to the triangulation T (P ).The Poincaré inequality in a triangle with the constant C P := 1/j 1,1 and the first positive root j 1,1 ≈ 3.8317 of the Bessel function J 1 from [24, Thm.2.1] allows for This and the Pythagoras theorem (with f − Π 0 f ⊥ P 0 ( T (P )) in L 2 (P )) show It remains to bound the term Π 0 f 2 L 2 (P )) .The assumption on f reads j∈J E(j) f ds = j∈J |E(j)|x j = 0 for a subset J ⊂ {1, . . ., m} so that 0 ∈ conv{|E(1)|x 1 , . . ., |E(m)|x m }.It follows 0 ∈ conv{x 1 , . . ., x m } and it is known that this implies for a constant M = 2) in the form x j = f j + a j to deduce from a triangle inequality and (2.5) that This shows that This and the Cauchy-Schwarz inequality imply the first two estimates in with the definition of h P and t min in the end.The inequality T (j Lemma 2.7] and the Cauchy-Schwarz inequality show, for j = 1, . . ., m, that The combination of the previous three displayed estimates result in This and (2.4) conclude the proof with the constant C 2 PF = (M + 1/4)(t max /t min ) + C 2 P .In the nonconforming VEM, the finite-dimensional space V h is a subset of the piecewise Sobolev space The piecewise H 1 seminorm (piecewise with respect to T hidden in the notation for brevity) reads

Local virtual element space
The first nonconforming virtual element space [3] is a subspace of harmonic functions with edgewise constant Neumann boundary values on each polygon.The extended nonconforming virtual element space [1,17] reads V h (P ) := v h ∈ H 1 (P ) : ∆v h ∈ P 1 (P ) and ∀E ∈ E(P ) ∂v h ∂n P E ∈ P 0 (E) .
(2.6) Definition 2.2 (Ritz projection).Let Π ∇ 1 be the Ritz projection from H 1 (P ) onto the affine functions P 1 (P ) in the H 1 seminorm defined, for v h ∈ H 1 (P ), by (∇Π ∇ 1 v h − ∇v h , ∇χ) L 2 (P ) = 0 for all χ ∈ P 1 (P ) and Remark 2 (integral mean).For P ∈ T and f ∈ H 1 (P ), ∇Π ∇ 1 f = Π 0 ∇f .(This follows from (2.7.a) and the definition of the L 2 projection operator Π 0 (acting componentwise) onto the piecewise constants P 0 (P ; R 2 ).) Remark 3 (representation of Π ∇ 1 ).For P ∈ T and f ∈ H 1 (P ), the Ritz projection Π ∇ 1 f reads The enhanced virtual element spaces [1,17] are designed with a computable L 2 projection Π 1 onto P 1 (T ).The resulting local discrete space under consideration throughout this paper reads The point in the selection of V h (P ) is that the Ritz projection Π ∇ 1 v h coincides with the L 2 projection Π 1 v h for all v h ∈ V h (P ).The degrees of freedom on P are given by Proof.Let E(1), . . ., E(m) be an enumeration of the edges E(P ) of the polygonal domain P in a consecutive way as depicted in Figure 2.2.a and define W (P ) := P 1 (P ) × P 0 (E(1)) × • • • × P 0 (E(m)).Recall V h (P ) from (2.6) and identify the quotient space V h (P )/R ≡ f ∈ V h (P ) : ∂P f ds = 0 with all functions in V h (P ) having zero integral over the boundary ∂P of P .Since the space V h (P ) consists of functions with an affine Laplacian and edgewise constant Neumann data, the map is well-defined and linear.The compatibility conditions for the existence of a solution of a Laplacian problem with Neumann data show that the image of S is equal to . ., g m ) ∈ W (P ) : (The proof of this identity assumes the compatible data (f 1 , g 1 , . . ., g m ) from the set on the right-hand side and solves the Neumann problem with a unique solution u in V h (P )/R and S u = (f 1 , g 1 , . . ., g m ).)It is known that the Neumann problem has a unique solution up to an additive constant and so S is a bijection and the dimension m + 2 of V h (P )/R is that of R(S).
The linearly independent functions ψ 3 , . . ., ψ m+2 belong to V h (P ) and so dim(V h (P )) ≥ m.Since V h (P ) ⊂ V h (P ) and three linearly independent conditions (1 − Π ∇ 1 )v h ⊥ P 1 (P ) in L 2 (P ) are imposed on V h (P ) to define V h (P ), dim(V h (P )) ≤ m.This shows that dim(V h (P )) = m and hence, the linear functionals dof E = − E • ds for E ∈ E(P ) form a dual basis of V h (P ).This concludes the proof of (b).
Remark 4 (stability of L 2 projection).The L 2 projection Π k for k = 0, 1 is H 1 and L 2 stable in V h (P ), in the sense that any v h in V h (P ) satisfies (The first inequality follows from the definition of Π k .The orthogonality in (2.9) and the definition of Π 1 imply that the Ritz projection Π ∇ 1 and the L 2 projection Π 1 coincide on the space V h (P ) for P ∈ T .This with the definition of the Ritz projection Π ∇ 1 verifies the second inequality.) Definition 2.4 (Fractional order Sobolev space [14]).Let α := (α 1 , α 2 ) denote a multi-index with α j ∈ N 0 for j = 1, 2 and |α| := α 1 + α 2 .For a real number m with 0 < m < 1, define with v α as the partial derivative of v of order α.Define the seminorm | • | 1+m and Sobolev-Slobodeckij norm • 1+m by Proposition 2.5 (approximation by polynomials [29, Thm.6.1]).Under the assumption (M2), there exists a positive constant C apx (depending on ρ and on the polynomial degree k) such that, for every v ∈ H m (P ), the L 2 projection Π k (P ) on (2.12)

Global virtual element space
Define the global nonconforming virtual element space, for any T ∈ T, by Example 2.2.If each polygonal domain P is a triangle, then the finite-dimensional space V h coincides with CR-FEM space.(Since the dimension of the vector space V h (P ) is three and P 1 (P ) ⊂ V h (P ), V h (P ) = P 1 (P ) for P ∈ T .) Lemma 2.6.There exists a universal constant C F (that depends only on ρ from (M2)) such that, for all T ∈ T, any v h ∈ V h from (2.13) satisfies Proof.Recall from Remark 1 that T is a shape regular sub-triangulation of T into triangles.Since V h ⊂ H 1 ( T ) and the Friedrichs' inequality holds for all functions in H 1 ( T ) [14,Thm. 10.6.16],there exists a positive constant C F such that the (first) inequality holds in The (second) equality follows for v h ∈ H 1 (P ) with P ∈ T .
Lemma 2.6 implies that the seminorm 1,pw in V h with mesh-size independent equivalence constants.

Interpolation
Definition 2.7 (interpolation operator).Let (ψ E : E ∈ E) be the nodal basis of V h defined by dof E (ψ E ) = 1 and dof F (ψ E ) = 0 for all other edges F ∈ E \ {E}.The global interpolation operator Since a Sobolev function v ∈ V has traces and the jumps [v] E vanish across any edge E ∈ E, the interpolation operator I h is well-defined.Recall ρ from (M2), C PF from Lemma 2.1, and C apx from Proposition 2.5.

Theorem 2.8 (interpolation error). (a)
There exists a positive constant C Itn (depending on ρ) such that any v ∈ H 1 (P ) and its interpolation I h v ∈ V h (P ) satisfy (c) The positive constant Proof of (a).The boundedness of the interpolation operator in V h (P ) is mentioned in [17] with a soft proof in its appendix.The subsequent analysis aims at a clarification that C I depends exclusively on the parameter ρ in (M2).The elementary arguments apply to more general situations in particular to 3D.Given I h v ∈ V h (P ), q 1 := −∆I h v ∈ P 1 (P ) is affine and An integration by parts leads to with q 1 ∈ P 1 (P ) and with the Cauchy inequality in the last step.Remark 2 and 3 on the Ritz projection, and the definition of (2.17) The function with (2.17) in the last step.Let φ c ∈ S 1 0 ( T (P )) := {w ∈ C 0 (P ) : w| T (E) ∈ P 1 (T (E)) for all E ∈ E(P )} denote the piecewise linear nodal basis function of the interior node c with respect to the triangulation

b for an illustration of T (P )). An inverse estimate
) for all f 1 ∈ P 1 ( T (P )) on the triangle T (E) := conv(E ∪{c}) holds with the universal constant C 1 .A constructive proof computes the mass matrices for T with and without the weight φ c to verify that the universal constant C 1 does not depend on the shape of the triangle T (E).This implies with an integration by parts for φ c q 1 ∈ H 1 0 (P ) and I h v in the last step.The mesh-size independent constant C 2 in the standard inverse estimate depends merely on the angles in the triangle T (E), E ∈ E(P ), and so exclusively on ρ.With . This and (2.19) lead to (2.20) The combination with (2.16)-(2.18)proves Proof of (b).The identity (2.17) reads Π 0 ∇(1 − I h )v = 0 and the triangle inequality results in This and the boundedness of the interpolation operator I h lead to with Remark 2 in the last step.The combination of (2.21) and (2.22) proves the first part of (b).
The identity follows from Lemma 2.1.a.This concludes the proof of (b).
Proof of (c).This is an immediate consequence of the part (b) with (2.12) and the Poincaré-Friedrichs inequality for v − I h v (from above) in Lemma 2.1.a.

Preliminary estimates
This subsection formulates the discrete problem along with the properties of the discrete bilinear form such as boundedness and a Gårding-type inequality.

The discrete problem
Choose the stability term S P (u h , v h ) as a symmetric positive definite bilinear form on V h (P ) × V h (P ) for a positive constant C s independent of P and h P satisfying For some positive constant approximation A P of A over P and the number N P := |E(P )| of the degrees of freedom (2.10) of V h (P ), a standard example of a stabilization term from [4], [36,Sec. 4.3] with a scaling coefficient A P reads Note that an approximation A P is a positive real number (not a matrix) and can be chosen as √ a 0 a 1 with the positive constants a 0 and a 1 from (A2).For f ∈ L 2 (Ω) and The sum over all the polygonal domains P ∈ T reads Remark 5 (polygonal mesh with small edges).The conditions (M1)-(M2) are well established and apply throughout the paper.The sub-triangulation T may not be shape-regular without the edge condition |E| ≥ ρh P for an edge E ∈ T (P ) and P ∈ T , but satisfies the maximal angle condition and the arguments employed in the proof of [8,Lemma 6.3] can be applied to show (2.20) in Theorem 2.8.a.For more general star-shaped polygon domains with short edges, the recent anisotropic analysis [8,15,18] indicates that the stabilization term has to be modified as well to avoid a logarithmic factor in the optimal error estimates.

Properties of the discrete bilinear form
The following proposition provides two main properties of the discrete bilinear form B h .
Proposition 3.1.There exist positive universal constants M, α and a universal nonnegative constant β depending on the coefficients A, b, γ such that Proof of (a).The upper bound of the coefficients from the assumption (A1), the Cauchy-Schwarz inequality, the stability (2.11) of Π 1 , and the definition (3.5) of the stabilization term imply the boundedness of B h with The details of the proof follow as in [6, Lemma 5.2] with the constant C F from Lemma 2.6.

Proof of (b). The first step shows that a
).This orthogonality, the assumption (A2), and the definition of the stability term (3.5) with the constant (3.9) The Cauchy-Schwarz inequality, (2.11), and the Young inequality lead to The combination of (3.9)-(3.10)proves This concludes the proof of (b) with α = α0 2 and β = ).The discrete space V h of the nonconforming VEM is endowed with the natural norm • h := a h (•, •) 1/2 induced by the scalar product a h .The boundedness of a h is proven in (a), while (3.9) shows the converse estimate in the equivalence

Consistency error
This subsection discusses the consistency error between the continuous bilinear form B and the corresponding discrete bilinear form B h .Recall the definition

(consistency). (a)
There exists a positive constant C cst (depending only on ρ) such that any v ∈ H 1 (Ω) and w h ∈ V h satisfy The definition of B P and B P h show The term T 1 in (3.13) is defined as the difference of the contributions from a P and a P h .Their definitions prove the equality (at the end of the first line below) and the definition of Π 1 prove the next equality in The last inequality follows from the Cauchy-Schwarz inequality, the Lipschitz continuity of A, and the stabilities ∇Π 1 v h L 2 (P ) ≤ ∇v h L 2 (P ) and ∇(1 − Π 1 )w h L 2 (P ) ≤ ∇w h L 2 (P ) from Remark 4. Similar arguments apply to T 2 from the differences of b P and b P h , and T 3 from those of c P and c P h in (3.13).This leads to The inequality for the last step in T 2 follows from the Cauchy-Schwarz inequality, the Lipschitz continuity of b, the estimate The combination of the above estimates shows (3.11).The proof of (3.12) adapts the arguments in the above analysis of T 3 and the definition of osc 1 (f, T ) in Subsection 2.1 for the proof of This concludes the proof.

Nonconformity error
Enrichment operators play a vital role in the analysis of nonconforming finite element methods [12].For any v h ∈ V h , the objective is to find a corresponding function Jv h ∈ H 1 0 (Ω).The idea is to map the VEM nonconforming space into the Crouzeix-Raviart finite element space with respect to the shape-regular triangulation T from Remark 1.Let ψ E be the edge-oriented basis functions of CR 1 0 ( T ) with ψ E (mid E) = 1 and ψ E (mid F ) = 0 for all other edges F ∈ E \ {E}.Define the interpolation operator The definition of  The following lemma describes the construction of a modified companion operator J : V h → H 1 0 (Ω), which is a right-inverse of the interpolation operator I h from Definition 2.7.

Lemma 3.3 (conforming companion operator).
There exists a linear map J : There exists an operator Recall that T (P ) is a shape-regular triangulation of P into a finite number of triangles.For each Let b T be extended by zero outside T and, for P ∈ T , define b P := 20 9 with − P b P dx = 1 and ∇b P L 2 (P ) h −1 P |P | 1/2 ≈ 1.Let v P ∈ P 1 (T ) be the Riesz representation of the linear functional P 1 (T ) → R defined by w 1 → (v h − v, w 1 ) L 2 (Ω) for w 1 ∈ P 1 (T ) in the Hilbert space P 1 (T ) endowed with the weighted L 2 scalar product (b P •, •) L 2 (P ) .Hence v P exists uniquely and satisfies Π 1 (v h − v) = Π 1 (b P v P ).Given the bubble-functions (b P : P ∈ T ) from (3.17) and the above functions (v P : P ∈ T ) for v h ∈ V h , define Proof of (a).Since b P vanishes at any x ∈ E ∈ E, it follows for any E ∈ E that where the definition of v = J v CR , (a'), and v CR = I CR v h lead to the second, third, and fourth equality.This proves (a).
Proof of (b).An integration by parts and (a) show, for all v h ∈ V h with Jv h from (3.18), that Since this holds for all P ∈ T , it proves (b).
Proof of (c).This is Π 1 v h = Π 1 Jv h and guaranteed by the design of J in (3.18).

Proof of (d).
This relies on the definition of J in (3.18) and J with (c').Since (a) allows for ∂P (v h − Jv h ) ds = 0, the Poincaré-Friedrichs inequality from Lemma 2.1.aimplies

Hence it remains to prove |v
show the first and second inequality in with (b') for in the last step.The equivalence of norms in the finite-dimensional space P 1 (P ) assures the existence of a positive constant C b , independent of h P , such that any χ ∈ P 1 (P ) satisfies the inverse inequalities The equality results from Π 1 (v h − v) = Π 1 (v P b P ) and v P ∈ P 1 (T ), while the last step is the Cauchy-Schwarz inequality.Consequently, Proof of I h J = id in V h .Definition 2.7 and Lemma 3.3.ashow, for all v h ∈ V h , that This concludes the proof of Lemma 3.3.
Since V h is not a subset of H 1 0 (Ω) in general, the substitution of discrete function v h in the weak formulation leads to a nonconformity error.
The first term on the right-hand side of (3.25) involves the factor The last inequality follows from the Lipschitz continuity of the coefficients A and b, and the estimate (2.12).Lemma 3.3.dleads to the estimates The substitution of the previous estimates in (3.25) with h max ≤ 1 (from δ ≤ 1 by assumption) and the regularity (1.5) show This concludes the proof of Lemma 3.4.a.

Proof of (b). The solution Φ
The arguments in the proof of (a) lead to the bound (3.24) with The remaining analogous details are omitted in the proof of Lemma 3.4.bfor brevity.

A priori error analysis
This section focuses on the stability, existence, and uniqueness of the discrete solution u h .The a priori error analysis uses the discrete inf-sup condition.

Existence and uniqueness of the discrete solution
Theorem 4.1 (stability).There exist positive constants δ ≤ 1 and C stab (depending on α, β, σ, ρ, and C F ) such that, for all T ∈ T(δ) and for all f ∈ L 2 (Ω), the discrete problem (3.8) has a unique solution u h ∈ V h and In the first part of the proof, suppose there exists some solution u h ∈ V h to the discrete problem (3.8) for some f ∈ L 2 (Ω).(This is certainly true for all f ≡ 0 ≡ u h , but will be discussed for all those pairs at the end of the proof and shall lead to the uniqueness of discrete solutions.)Since u h satisfies a Gårding-type inequality in Proposition 3.1.b, This, (2.14), and the definition of the dual norm in (3.12) lead to Given g := u h ∈ L 2 (Ω), let Φ ∈ V ∩ H 1+σ (Ω) solve the dual problem L * Φ = g and let I h Φ ∈ V h be the interpolation of Φ from Subsection 2.4.Elementary algebra shows Rewrite a part of the third term corresponding to diffusion on the right-hand side of (4.2) as The Cauchy-Schwarz inequality in the semi-scalar product S P (•, •), and (3.5) with the upper bound A ∞ for the coefficient A in a P (•, •) lead to the estimate with Theorem 2.8.b followed by (2.12) in the final step.This and Theorem 2.8 imply that The terms b P − b P h and c P − c P h are controlled by The combination of the previous four displayed estimates with Lemma 2.6 leads to an estimate for P .The sum over all polygonal domains P ∈ T reads with a universal constant C d .The bound for (4.2) results from Lemma 3.4.bfor the first term, the boundedness of B pw (with a universal constant ) and (2.15) for the second term, (4.4) for the third term, and Theorem 2.8.a for the last term on the right-hand side of (4.2).This shows This and the regularity estimate (1.5) lead to

The substitution of this in (4.1) proves
For all 0 < h max ≤ δ := ( is well-defined.This leads in (4.5) to In the last part of the proof, suppose f h ≡ 0 and let u h be any solution to the resulting homogeneous linear discrete system.The stability result (4.6) proves u h ≡ 0. Hence, the linear system of equations (3.8) has a unique solution and the coefficient matrix is regular.This proves that there exists a unique solution u h to (3.8) for any right-hand side f h ∈ V * h .The combination of this with (4.6) concludes the proof.
An immediate consequence of Theorem 4.1 is the following discrete inf-sup estimate.

Theorem 4.2 (discrete inf-sup).
There exist 0 < δ ≤ 1 and β 0 > 0 such that, for all T ∈ T(δ), Proof.Define the operator L h : The stability Theorem 4.1 can be interpreted as follows: For any f h ∈ V * h there exists u h ∈ V h such that L h u h = f h and The discrete problem B h (u h , •) = (f h , •) has a unique solution in V h .Therefore, f h and u h are in one to one correspondence and the last displayed estimate holds for any u h ∈ V h .The infimum over u h ∈ V h therein proves (4.7) with β 0 = C −1 stab .

A priori error estimates
This subsection establishes the error estimate in the energy norm | • | 1,pw and in the L 2 norm.The discrete inf-sup condition allows for an error estimate in the H 1 norm and an Aubin-Nitsche duality argument leads to an error estimate in the L 2 norm.Recall u ∈ H 1 0 (Ω) is a unique solution of (1.8) and u h ∈ V h is a unique solution of (3.8).Recall the definition of the bilinear form s h (•, •) from Section 3.1 and define the induced seminorm Theorem 4.3 (error estimate).Set σ := A∇u + bu ∈ H(div, Ω).There exist positive constants C 4 , C 5 , and δ such that, for all T ∈ T(δ), the discrete problem (3.8) has a unique solution Proof.
Step 1 (initialization).Let I h u ∈ V h be the interpolation of u from Definition 2.7.The discrete inf-sup condition (4.7) for Step 2 (error estimate for |u − u h | 1,pw ).Rewrite the last equation with the continuous and the discrete problem (1.8) and (3.8) as This equality is rewritten with the definition of B(u, v) in (1.7), the definition of B h (I h u, v h ) in Section 3.1, and with f h = Π 1 f .Recall v := Jv h ∈ V from Lemma 3.3 and recall ∇ pw Π 1 I h u = Π 0 ∇u from (2.17).This results in Abbreviate w := v − Π 1 v h and observe the orthogonalities ∇ pw w ⊥ P 0 (T ; R 2 ) in L 2 (Ω; R 2 ) and w ⊥ P 1 (T ) in L 2 (Ω) from Lemma 3.3.b-cand the definition of Π 1 with Π

and the Poincaré-Friedrichs inequality for
Elementary algebra and the above orthogonalities prove that
Step 3 (duality argument).To prove the bound for u − u h in the L 2 norm with a duality technique, let g := I h u − u h ∈ L 2 (Ω).The solution Φ ∈ H 1 0 (Ω) ∩ H 1+σ (Ω) to the dual problem (1.4) satisfies the elliptic regularity (1.5), Step 4 (error estimate for u − u h L 2 (Ω) ).Let I h Φ ∈ V h be the interpolation of Φ from Definition 2.7.Elementary algebra reveals the identity The bound (4.4) with g as the first argument shows This controls the third term in (4.14), Lemma 3.4.bcontrols the first term, the boundedness of B pw and the interpolation error estimate (2.15) control the second term on the right-hand side of (4.14).This results in It remains to bound B h (g, I h Φ).The continuous and the discrete problem (1.8) and (3.8) imply The definition of B h and Π 0 lead to The bound for the stability term as in (4.12) is Step 5 (oscillation).The last term in (4.16) is of optimal order O(h 1+σ max ), but the following arguments allow to write it as an oscillation.Recall the bubble-function b T | P := b P ∈ H 1 0 (P ) from (3.17) extended by zero outside P .Given Ψ := Φ − Π 1 I h Φ, let Ψ 1 ∈ P 1 (T ) be the Riesz representation of the linear functional P 1 (T ) → R defined by w 1 → (Ψ, w 1 ) L 2 (Ω) in the Hilbert space P 1 (T ) endowed with the weighted scalar product (b allow the rewriting of the latter identity as It remains to control the terms h −1 T (Ψ − b T Ψ 1 ) L 2 (Ω) and |b T Ψ| 1,pw .Since the definition of I h and the definition of Π ∇ 1 with Π 1 = Π ∇ 1 in V h imply ∂P Ψ ds = ∂P (Φ − Π 1 I h Φ) ds = 0, this allows the Poincaré-Friedrichs inequality for Ψ from Lemma 2.1.aon each P ∈ T .This shows with Theorem 2.8.b and (2.12) in the last inequality.Since b P Ψ 1 ∈ H 1 0 (P ) for P ∈ T , the Poincaré-Friedrichs inequality from Lemma 2.1.aleads to The first estimate in (3.20), the identity Π 1 (b T Ψ 1 ) = Π 1 Ψ, and the Cauchy-Schwarz inequality imply . The second estimate in (3.21) followed by the first estimate in (3.20) leads to the first inequality and the arguments as above lead to the second inequality in L 2 (P ) from above in the last step.The combination of the previous displayed estimate and (4.18)-(4.20)results with A lower bound of the stability term (3.5) and the assumption (A2) imply The analogous arguments for u−Π 1 u h L 2 (Ω) , (4.25), and the estimate for |u h | s prove the bound (4.8) for the term h −σ max u − Π 1 u h L 2 (Ω) .This concludes the proof of Theorem 4.3.

A posteriori error analysis
This section presents the reliability and efficiency of a residual-type a posteriori error estimator.

Residual-based explicit a posteriori error control
Recall u h ∈ V h is the solution to the problem (3.8), and the definition of jump [•] E along an edge E ∈ E from Section 2. For any polygonal domain P ∈ T , set , and Ξ T := ( P ∈T Ξ 2 P ) 1/2 .The following theorem provides an upper bound to the error u − u h in the H 1 and the L 2 norm.Recall the elliptic regularity (1.5) with the index 0 < σ ≤ 1, and recall the assumption h max ≤ 1 from Subsection 2.1.

Theorem 5.1 (reliability).
There exist positive constants C rel1 and C rel2 (both depending on ρ) such that The proof of this theorem in Subsection 5.3 relies on a conforming companion operator elaborated in the next subsection.The upper bound in Theorem 5.1 is efficient in the following local sense, where ω E := int(∪T (E)) denotes the patch of an edge E and consists of the one or the two neighbouring polygons in the set T (E) := {P ∈ T : E ⊂ ∂P } that share E. Recall σ = A∇u + bu from Subsection 4.2 and the data-oscillation osc 1 (f, P ) := h P (1 − Π 1 )f L 2 (P ) from Subsection 2.1.
Theorem 5.2 (local efficiency up to oscillation).The quantities η P , ζ P , Λ P , and Ξ P from Theorem 5.1 satisfy (5.3) The proof of Theorem 5.2 follows in Subsection 5.4.The reliability and efficiency estimates in Theorem 5.1 and 5.2 lead to an equivalence up to the approximation term Recall the definition of |u h | s from Subsection 4.2.In this paper, the norm | • | 1,pw in the nonconforming space V h has been utilised for simplicity and one alternative is the norm • h from Remark 6 induced by a h .Then it appears natural to have the total error with the stabilisation term as The point is that Theorem 4.3 assures that total error + apx converges with the expected optimal convergence rate.
This proves the first inequality in the assertion.Theorem 5.1, the estimates in Subsection 5.3.3.1, and the definition of |u h | s show total error estimator.The first of the terms in apx is The definition of σ and σ h plus the triangle and the Cauchy-Schwarz inequality show The upper bound is estimator as mentioned above.Since the term (1 Section 5 establishes the a posteriori error analysis of the nonconforming VEM.Related results are known for the conforming VEM and the nonconforming FEM. Remark 7 (comparison with nonconforming FEM).Theorem 5.1 generalizes a result for the nonconforming FEM in [19,Thm. 3.4] from triangulations into triangles to those in polygons (recall Example 2.2).The only difference is the extra stabilization term that can be dropped in the nonconforming FEM.
Remark 8 (comparison with conforming VEM).The volume residual, the inconsistency term, and the stabilization also arise in the a posteriori error estimator for the conforming VEM in [16,Thm. 13].But it also includes an additional term with normal jumps compared to the estimator (5.1).The extra nonconformity term in this paper is caused by the nonconformity V h ⊂ V in general.

Enrichment and conforming companion operator
The link from the nonconforming approximation u h ∈ V h to a global Sobolev function in H 1 0 (Ω) can be designed with the help of the underlying refinement T of the triangulation T (from Section 2).The interpolation I CR : V + V h → CR 1 0 ( T ) in the Crouzeix-Raviart finite element space CR 1 0 ( T ) from Subsection 3.4 allows for a right-inverse J .A companion operator J • I CR : Define an enrichment operator E pw : P 1 ( T ) → S 1 0 ( T ) by averaging nodal values: For any vertex z in the refined triangulation T , let T (z) = {T ∈ T : z ∈ T } denote the set of | T (z)| ≥ 1 many triangles that share the vertex z, and define for an interior vertex z (and zero for a boundary vertex z according to the homogeneous boundary conditions).This defines E pw v 1 at any vertex of a triangle T in T , and linear interpolation then defines E pw v 1 in T ∈ T , so that E pw v 1 ∈ S 1 0 ( T ).Huang et al. [31] design an enrichment operator by an extension of [32] to polygonal domains, while we deduce it from a sub-triangulation.The following lemma provides an approximation property of the operator E pw .
Lemma 5.4.There exists a positive constant C En that depends only on the shape regularity of T such that any v 1 ∈ P 1 (T ) satisfies

.7)
Proof.There exists a positive constant C En independent of h and v 1 [32, p. 2378] such that .
Note that any edge E ∈ E is unrefined in the sub-triangulation T .Since v 1|P ∈ H 1 (P ) is continuous in each polygonal domain P ∈ T and h T ≤ h P for all T ∈ T (P ), the above inequality reduces to (5.7).This concludes the proof.
Recall the L 2 projection Π 1 onto the piecewise affine functions P 1 (T ) from Section 2. An enrichment operator E pw • Π 1 : V h → H 1 0 (Ω) acts as displayed Abbreviate w := v − Π 1 I h v and σ h := A∇ pw Π 1 u h + bΠ 1 u h .This and (5.9) simplify with P ∇w dx = 0 for any P ∈ T from (2.17) in the last step.Recall the notation η P , Λ P , and ζ P from Subsection 5.1.The Cauchy-Schwarz inequality and Theorem 2.8.b followed by The upper bound A ∞ of the coefficient A, (3.5), and the Cauchy-Schwarz inequality for the stabilization term lead to the first inequality in The second inequality in (5.13) follows as in (4.3) and with (1 − Π 0 )∇v L 2 (P ) ≤ 1. Recall the boundedness constant M b of B pw from Subsection 4.1 and deduce from (5.7) and the definition of Ξ T from Subsection 5.1 that The substitution of (5.10)-(5.14) in (5.8) reveals that with The combination of (4.24), (5.15) and (5.7) leads in the triangle inequality

Reliable L 2 error control
Recall I CR from (3.14) and J from the proof of Lemma 3.3, and define The substitution of v := u − E 2 u h ∈ V in the dual problem shows The algebra in (5.8)-(5.10)above leads with v = Ψ to the identity (5.17) The definition of I CR and J proves the first and second equality in This and an integration by parts imply P ∇(u h − E 2 u h ) dx = 0 for all P ∈ T .Hence Definition 2.2 of Ritz projection Π ) and the definition of B pw in the last term of (5.17) result with elementary algebra in The triangle inequality and (c') from the proof of Lemma 3.3 imply the first inequality in (5. 19) The second estimate in (5.19) follows from E 1 u h ∈ V , the third is a triangle inequality, and eventually The second inequality in (5.20) follows from the Poincaré-Friedrichs inequality in Lemma 2.1.afor Π 1 u h − E 2 u h with ∂P (Π 1 u h − E 2 u h ) ds = 0 (from above); the constant (5.19) and h P ≤ h σ P (recall h max ≤ 1).Lemma 5.4 with v 1 = Π 1 u h and (4.24) in (5.20) show  The arguments in the proof of (5.20)-(5.21)also lead to P ∈T h P (ζ P + Ξ P ). (5.23) The combination of (4.25), (5.22)-(5.23)and the triangle inequality This concludes the proof of the L 2 error estimate in Theorem 5.1.

Estimator for u − Π 1 u h
The triangle inequality with (5.1) and (4.24) provide an upper bound for H 1 error The same arguments for an upper bound of the L 2 error in Theorem 5.1 show that The numerical experiments do not display C rel1 and C rel2 , and directly compare the error H1e := |u − Π 1 u h | 1,pw in the piecewise H 1 norm and the error L2e := u − Π 1 u h L 2 (Ω) in the L 2 norm with the upper bound H1µ and L2µ (see, e.g., Figure 6.3).

Motivation and discussion of apx
We first argue that those extra terms have to be expected and utilize the abbreviations σ := A∇u + bu and g := f − γu for the exact solution u ∈ H 1 0 (Ω) to (1.8), which reads (σ, ∇v) L 2 (Ω) = (g, v) L 2 (Ω) for all v ∈ H 1 0 (Ω).(5.24) Recall the definition of s h (•, •) from Subsection 3.1.The discrete problem (3.8) with the discrete solution u h ∈ V h assumes the form for σ h := A∇Π 1 u h + bΠ 1 u h , and g h := f − γΠ 1 u h .Notice that σ h and g h may be replaced in (5.25) by Π 0 σ h and Π 1 g h because the test functions ∇Π 1 v h and Π 1 v h belong to P 0 (T ; R 2 ) and P 1 (T ) respectively.In other words, the discrete problems (3.8) and (5.25) do not see a difference of σ h and g h compared to Π 0 σ h and Π 1 g h and so the errors σ h −Π 0 σ h and g h −Π 1 g h may arise in a posteriori error control.This motivates the a posteriori error term σ h − Π 0 σ h L 2 (Ω) = Λ T as well as the approximation terms σ − Π 0 σ and g − Π 1 g on the continuous level.The natural norm for the dual variable σ is L 2 and that of g is H −1 and hence their norms form the approximation term apx as defined in Subsection 5.1.
Example 5.1 (b = 0).The term (1 − Π 0 )σ may not be visible in case of no advection b = 0 at least if A is piecewise constant.Suppose A ∈ P 0 (T ; R 2×2 ) and estimate If A is not constant, there are oscillation terms that can be treated properly in adaptive meshrefining algorithms, e.g., in [27].
Example 5.2 (γ piecewise constant).While the data approximation term osc 1 (f, T ) [10] is widely accepted as a part of the total error in the approximation of nonlinear problems, the term osc is of higher order and may even be absorbed in the overall error analysis for a piecewise constant coefficient γ ∈ P 0 (T ).In the general case γ ∈ L ∞ (Ω)\P 0 (T ), however, osc 1 (u, T ) leads in particular to terms with γ − Π 0 γ L ∞ (Ω) .

Higher-order nonconforming VEM
The analysis applied in Theorem 5.1 can be extended to the nonconforming VEM space of higher order k ∈ N (see [17,Sec. 4] for the definition of discrete space).Since the projection operators ∇Π ∇ k and Π k−1 ∇ are not the same for general k, and the first operator does not lead to optimal order of convergence for k ≥ 3, the discrete formulation uses Π k−1 ∇ (cf.[6,Rem. 4.3] for more details).The definition and approximation properties of the averaging operator E pw extend to the operator E k : P k ( T ) → H 1 0 (Ω) (see [32, p. 2378] for a proof).The identity (5.9) does not hold in general, but algebraic calculations lead to The analysis developed for the upper bound of L 2 norm also extends to the general case.The model problem is chosen in 2D for the simplicity of the presentation.The results of this work can be extended to the three-dimensional case with appropriate modifications.The present analysis holds for any higher regularity index σ > 0 and avoids any trace inequality for higher derivatives.This is possible by a medius analysis in the form of companion operators [26].

Inhomogeneous boundary data
The error estimator for general Dirichlet condition u| ∂Ω = g ∈ H 1/2 (∂Ω) can be obtained with some modifications of [33] in Theorem 5.1.The only difference is in the modified jump contributions of the boundary edges in the nonconformity term

Proof of Theorem 5.2
Recall the notation σ = A∇u + bu and σ h = A∇Π 1 u h + bΠ 1 u h from Subsection 5.3.
Proof of (5.3).The upper bound (3.5) for the stabilisation term and the triangle inequality show This concludes the proof of (5.3).
Proof of (5.5).The definition of Λ P , Π 0 , and the triangle inequality lead to The upper bound A ∞ and b ∞ for the coefficients and the triangle inequality lead to ).This followed by (5.3) concludes the proof of (5.5).
Recall the bubble-function b T | P = b P supported on a polygonal domain P ∈ T from (3.17) as the sum of interior bubble-functions supported on each triangle T ∈ T (P ).
Proof of (5.4).Rewrite the term and denote R P := R| P and θ P := θ| P .The definition of B pw (u − Π 1 u h , v) and the weak formulation (5.29) Since b P R P belongs to H 1 0 (Ω) (extended by zero outside P ), v := b P R P ∈ V is admissible in (5.29).An integration by parts proves that (Π 0 σ h , ∇(b P R P )) L 2 (P ) = 0. Therefore, (5.29) shows The substitution of χ = R P = Π 1 (f − γΠ 1 u h )| P ∈ P 1 (P ) in (3.20) and the previous identity with the boundedness of B and the Cauchy-Schwarz inequality lead to the first two estimates in The last inequality follows from the definition of Λ P , and (3.21) with χ = R P .This proves that C . Recall η P from Subsection 5.1 and η P = h P f − γΠ 1 u h L 2 (P ) ≤ h P R P L 2 (P ) + h P θ P L 2 (P ) from the split in (5.28) and the triangle inequality.This and the previous estimate of h P R P L 2 (P ) show the first estimate in The second step results from the definition of θ P = (1 − Π 1 )(f − γΠ 1 u h )| P in (5.28) followed by the L 2 orthogonality of Π 1 , and the last step results from an elementary algebra with the triangle inequality and osc 1 (f − γu, P ) = h P (1 − Π 1 )(f − γu) L 2 (P ) from Subsection 5.1.The triangle inequality for the term u − Π 1 u h and the estimate of u h − Π 1 u h 1,P as in (5.27) lead to ).The combination of (5.3) and (5.5) in the last displayed estimate concludes the proof of (5.4).
Proof of (5.6).Recall for u ∈ H 1 0 (Ω) and u h ∈ V h that − E u ds and − E u h ds are well defined for all edges E ∈ E, and so the constant α E := − E (u − u h ) ds is uniquely defined as well.Since the jump of u − α E across any edge This and the triangle inequality show the first estimate in This with the mesh assumption h P ≤ ρ −1 |E| and (5.30) result in Since this holds for any edge E ∈ E(P ), the sum over all these edges and the bound (5.3) in the above estimate conclude the proof of (5.6).
Remark 9 (convergence rates of L 2 error control for 0 < σ ≤ 1).The efficiency estimates (5.4)-(5.6)with a multiplication of h 2σ P show that the local quantity h 2σ P (η 2 P + Λ 2 P + Ξ 2 P ) converges to zero with the expected convergence rate.
Remark 10 (efficiency up to stabilisation and oscillation for L 2 error control when σ = 1).For convex domains and σ = 1, there is even a local efficiency result that is briefly described in the sequel: The arguments in the above proof of (5.4)-(5.5)lead to The observation [Π 1 u h ] E = [Π 1 u h − u] E for the term Ξ P , the trace inequality, and the triangle inequality show, for any E ∈ E, that The bound (4.25) for the first term and the inverse estimate ∇χ L 2 (P ) ≤ C inv h −1 P χ L 2 (P ) for χ ∈ P k (P ) for the third term result in The mesh assumption (M2) implies that h . This and the above displayed inequality prove the efficiency estimate for h 2 P Ξ 2 P .

Numerical experiments
This section manifests the performance of the a posteriori error estimator and an associated adaptive mesh-refining algorithm with Dörfler marking [37].The numerical results investigate three computational benchmarks for the indefinite problem (1.1).

Adaptive algorithm
Input: initial partition T 0 of Ω.For = 0, 1, 2, . . .do 1.SOLVE.Compute the discrete solution u h to (3.8) with respect to T for = 0, 1, 2 . . .(cf. [5] for more details on the implementation).3. MARK.Mark the polygons P in a subset M ⊂ T with minimal cardinality and 4. REFINE -Refine the marked polygon domains by connecting the mid-point of the edges to the centroid of respective polygon domains and update T .(cf. Figure 6.1 for an illustration of the refinement strategy.)The adaptive algorithm is displayed for mesh adaption in the energy error H 1 .Replace estimator H1µ in the algorithm by L2µ (the upper bound in (5.2)) for local mesh-refinement in the L 2 error.Both uniform and adaptive mesh-refinement run to compare the empirical convergence rates and provide numerical evidence for the superiority of adaptive mesh-refinement.Note that uniform refinement means all the polygonal domains are refined.In all examples below, A P = 1 in (3.6).The numerical realizations are based on a MATLAB implementation explained in [35] with a Gauss-like cubature formula over polygons.The cubature formula is exact for all bivariate polynomials of degree at most 2n − 1, so the choice n ≥ (k + 1)/2 leads to integrate a polynomial of degree k exactly.The quadrature errors in the computation of examples presented below appear negligible for the input parameter n = 5.

Square domain (smooth solution)
This subsection discusses the problem (1.1) with the coefficients A = I, b = (x, y) and γ = x 2 +y 3 on a square domain Ω = (0, 1) 2 , and the exact solution   with f := Lu.Since the exact solution is not zero along the boundary ∂Ω, the error estimators are modified according to Subsection 5.3.3.4.Since γ − 1 2 div(b) = −5 < 0, the problem is noncoercive.Observe that with increase in number of iterations, refinement is more at the singularity as highlighted in Figure 6.4.Since the exact solution u is in H (5/3)− (Ω) for all > 0, from a priori error estimates the expected order of convergence in H 1 norm is 1/3 and in L 2 norm is at least 2/3 with respect to number of degrees of freedom for uniform refinement.Figure 6.5 shows that uniform refinement gives the sub-optimal convergence rate, whereas adaptive refinement lead to optimal convergence rates (1/2 for H 1 norm and 5/6 in L 2 norm).There is an internal layer around the circle centered at (0, 0) and of radius 0.25 where the second derivatives of u are large because of steep increase in the solution resulting in the large error at the beginning, and this gets resolved with refinement as displayed in Figure 6.6.

Conclusion
The three computational benchmarks provide empirical evidence for the sharpness of the mathematical a priori and a posteriori error analysis in this paper and illustrate the superiority of adaptive over uniform mesh-refining.The empirical convergence rates in all examples for the H 1 and L 2 errors coincide with the predicted convergence rates in Theorem 4.3, in particular, for the non-convex domain and reduced elliptic regularity.The a posteriori error bounds from Theorem 5.1 confirm these convergence rates as well.The ratio of the error estimator µ by the H 1 error e , sometimes called efficiency index, remains bounded up to a typical value 6; we regard this as a typical overestimation factor for the residual-based a posteriori error estimate.
Recall that the constant C reg has not been displayed so the error estimator µ does not provide a guaranteed error bound.Figure 6.8 and 6.9 display the four different contributions volume residual ( P η 2 P ) 1/2 , stabilization term ( P ζ 2 P ) 1/2 , inconsistency term ( P Λ 2 P ) 1/2 and the nonconformity term ( P Ξ 2 P ) 1/2 that add up to the error estimator µ .We clearly see that all four terms converge with the overall rates that proves that none of them is a higher-order term and makes it doubtful that some of those terms can be neglected.The volume residual clearly dominates the a posteriori error estimates, while the stabilisation term remains significantly smaller for the natural stabilisation (with undisplayed parameter one).The proposed adaptive mesh-refining algorithm leads to superior convergence properties and recovers the optimal convergence rates.This holds for the first example with optimal convergence rates in the large pre-asymptotic computational range as well as in the second with suboptimal convergence rates under uniform mesh-refining according to the typical corner singularity and optimal convergence rates for the adaptive mesh-refining.The third example with the Helmholtz equation and a moderate wave number shows certain moderate local mesh-refining in Figure 6.6 but no large improvement over the optimal convergence rates for uniform mesh-refining.The adaptive refinement generates hanging nodes because of the way refinement strategy is defined, but this is not troublesome in VEM setting as hanging node can be treated as a just another vertex in the decompostion of domain.However, an increasing number of hanging nodes with further mesh refinements may violate the mesh assumption (M2), but numerically the method seems robust without putting any restriction on the number of hanging nodes.The future work on the theoretical investigation of the performance of adaptive mesh-refining algorithm is clearly motivated by the successful numerical experiments.The aforementioned empirical observation that the stabilisation terms do not dominate the a posteriori error estimates raises the hope for a possible convergence analysis of the adaptive mesh-refining strategy with the axioms of adaptivity [20] towards a proof of optimal convergence rates: The numerical results in this section support this conjecture at least for the lowest-order VEM in 2D for indefinite non-symmetric second-order elliptic PDEs.

Figure 2 . 2 .
The constant C PF depends exclusively on the number m := |E(P )| of the edges in the polygonal domain P and the quotient of the maximal area divided by the minimal area of a triangle in the triangulation T (P ).

. 8 )(
The proof of (2.8) consists in the verification of (2.7):The equation (2.7.a) follows from Remark 2 with an integration by parts.The equation (2.7.b) follows from the definition of mid(∂P ) as the barycenter of ∂P .)

. 21 )
These estimates are completely standard on shape-regular triangles[2, p. 27] or[37]; so they hold on each T ∈ T and, by definition of b P , their sum is (3.20)-(3.21).The analysis of the term |v − Jv h | 1,pw starts with one P ∈ T and (3.18) for |v − Jv h | 1,P = |v P b P | 1,P ≤ C b h −1 P v P L 2 (P ) (3.22) with (3.21) in the last step.The estimate (3.20) leads to the first inequality in pw with ∂P (v − v h ) ds = 0 from (a) and hence the Poincaré-Friedrichs inequality for v − v h from Lemma 2.1.ain the last step.Recall |v−v h | 1,pw |v h | 1,pw from (3.19) to conclude |v−Jv h | 1,pw |v h | 1,pw from the previous displayed inequality.This concludes the proof of (d).

8 )
and the piecewise version B pw of B in the last step.The definition of B h from Subsection 3.1 and the discrete problem (3.8) with v h = I h v imply This and the regularity (5.16) result in u − E 2 u h L 2 (Ω) ≤ C 9 C * reg P ∈T h σ P (η P + ζ P + Λ P + Ξ P ).(5.22)

u
= 16x(1 − x)y(1 − y) arctan(25x − 100y + 50) with f = Lu.Since γ − 1 2 div(b) = x 2 + y 3 − 1 is not always positive on Ω, this is an indefinite problem.Initially, the error and the estimators are large because of an internal layer around the line 25x − 100y + 50 = 0 with large first derivative of u resolved after few refinements as displayed in Figure 6.2.

Figure 6 . 3 :
Figure 6.3: Convergence history plot of estimator µ and error e := u − Π 1 u h in the (a) piecewise H 1 norm (b) L 2 norm vs number ndof of degrees of freedom for both uniform and adaptive refinement

Figure 6 . 5 :
Figure 6.5: Convergence history plot of estimator µ and error e := u − Π 1 u h in the (a) piecewise H 1 norm (b) L 2 norm vs number ndof of degrees of freedom for both uniform and adaptive refinement

Figure 6 . 8 :Figure 6 . 9 :
Figure 6.8: Estimator components corresponding to the error H1e = |u − Π 1 u h | 1,pw of the adaptive refinement presented in Subsection 6.2-6.4 E denote the jump across an edge E ∈ E: For two neighboring polygonal domains P + and P − sharing a common edge E ∈ E(P + ) ∩ E(P − ), [v h ] E := v h|P + − v h|P − , where P + denote the adjoint polygonal domain with n P + |E = n E and P − denote the polygonal domain with n and any v ∈ H 1+σ (P ) with the local interpolationI h v| P ∈ V h (P ) satisfy v − I h v L 2 (P ) + h P |v − I h v| 1,P ≤ C I h 1+σ P |v| 1+σ,P .