{'fraction_non_alphanumeric': 0.047738362125807725, 'fraction_numerical': 0.01157853092443624, 'mean_word_length': 3.9286364227219552, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Residual based Error Estimate and Quasi-Interpolation on Polygonal Meshes for High Order BEM-based FEM

Only a few numerical methods can treat boundary value problems on polygonal and polyhedral meshes. The BEM-based Finite Element Method is one of the new discretization strategies, which make use of and benefits from the flexibility of these general meshes that incorporate hanging nodes naturally. The article in hand addresses quasi-interpolation operators for the approximation space over polygonal meshes. To prove interpolation estimates the Poincar\'e constant is bounded uniformly for patches of star-shaped elements. These results give rise to the residual based error estimate for high order BEM-based FEM and its reliability as well as its efficiency are proven. Such a posteriori error estimates can be used to gauge the approximation quality and to implement adaptive FEM strategies. Numerical experiments show optimal rates of convergence for meshes with non-convex elements on uniformly as well as on adaptively refined meshes.


Introduction
Adaptive finite element strategies play a crucial role in nowadays efficient implementations for the numerical solution of boundary value problems. Due to a posteriori error control, the meshes are only refined locally to reduce the computational cost and to increase the accuracy. For classical discretization techniques, as the Finite Element Method (FEM), triangular or quadrilateral (2D) and tetrahedral or hexahedral (3D) elements are applied and one has to take care to preserve the admissibility of the mesh after local refinements [29]. Polygonal (2D) and polyhedral (3D) meshes are, therefore, very attractive for such local refinements. They allow a variety of element shapes and they naturally treat hanging nodes in the discretization, which appear as common nodes on a more general element. But only a few methods can handle such meshes which also provide a more flexible and direct way to meshing for complex geometries and interfaces.
The BEM-based Finite Element Method is one of the new promising approaches applicable on general polygonal and polyhedral meshes. It was first introduced in [10] and analysed in [16,15]. The method makes use of local Trefftz-like trial functions, which are defined implicitly as local solutions of boundary value problems. These problems are treated by means of Boundary Element Methods (BEM) that gave the name. The BEM-based FEM has been generalized to high order approximations [25,32,34], mixed formulations with H(div) conforming approximations [14] as well as to convection-adapted trial functions [18]. Furthermore, the strategy has been applied to general polyhedral meshes [26] and time dependent problems [33]. For efficient computations fast FETI-type solvers have been developed for solving the resulting large linear systems of equations [17]. Additionally, the BEM-based FEM has shown its flexibility and applicability on adaptively refined polygonal meshes [31,35].
In the last few years, polygonal and polyhedral meshes attracted a lot of interest. New methods have been developed and conventional approaches were mathematically revised to handle them. The most prominent representatives for the new approaches beside of the BEM-based FEM are the Virtual Element Method [4] and the Weak Galerkin Method [30]. Strategies like discontinuous Galerkin [12] and the mimetic discretization techniques [5] are also considered on polygonal and polyhedral meshes. But there are only a few references to adaptive strategies and a posteriori error control. A posteriori error estimates for the discontinuous Galerkin method are given in [20]. To the best of our knowledge there is only one publication for the Virtual Element Method [11] and one for the Weak Galerkin Method [7]. The first mentioned publication deals with a residual a posteriori error estimate for a C 1 -conforming approximation space, and the second one is limited to simplicial meshes. For the mimetic discretization technique there are also only few references which are limited to low order methods, see the recent work [3].
The article in hand contributes to this short list of literature. In Section 2, some notation and the model problem are given. The regularity of polygonal meshes is discussed and useful properties are proven. Additionally, the BEM-based FEM with high order approximation spaces is reviewed. Section 3 deals with quasiinterpolation operators of Clément type for which interpolation estimates are shown on regular polygonal meshes with non-convex elements. These operators are used in Section 4 to derive a residual based error estimate for the BEM-based FEM on general meshes which is reliable and efficient. The a posteriori error estimate can be applied to adaptive mesh refinement that yields optimal orders of convergence in the numerical examples in Section 5.

Model problem and discretization
Let Ω ⊂ R 2 be a connected, bounded, polygonal domain with boundary Γ = Γ D ∪Γ N , where Γ D ∩ Γ N = ∅ and |Γ D | > 0. We consider the diffusion equation with mixed Figure 1: Two neighbouring elements in a polygonal mesh, nodes are marked with dots Dirichlet-Neumann boundary conditions where n Ω denotes the outer unite normal vector to the boundary of Ω. For simplicity, we restrict ourselves to piecewise constant and scalar valued diffusion a ∈ L ∞ (Ω) with 0 < a min ≤ a ≤ a max almost everywhere in Ω. Furthermore, we assume that the discontinuities of a are aligned with the initial mesh in the discretization later on. The treatment of continuously varying coefficients is discussed in [25], where the coefficient is approximated by a piecewise constant function or by a proper interpolation in a BEM-based approximation space of lower order. The usual notation for Sobolev spaces is utilized, see [1,21]. Thus, for an open subset ω ⊂ Ω, which is either a two dimensional domain or a one dimensional Lipschitz manifold, the Sobolev space H s (ω), s ∈ R is equipped with the norm · s,ω = · H s (ω) and semi-norm | · | s,ω = | · | H s (ω) . This notation includes the Lebesgue space of square integrable functions for s = 0, since H 0 (ω) = L 2 (ω). The L 2 (ω)-inner product is abbreviated to (·, ·) ω . For g D ∈ H 1/2 (Γ D ), g N ∈ L 2 (Γ N ) and f ∈ L 2 (Ω), the well known weak formulation of problem (1) reads and admits a unique solution, since the bilinear form b(u, v) = (a∇u, ∇v) Ω is bounded and coercive on Here, g D + V ⊂ H 1 (Ω) denotes the affine space which incorporates the Dirichlet datum. Note that the same symbol g D is used for the Dirichlet datum and its extension into H 1 (Ω). The domain Ω is decomposed into a family of polygonal discretizations K h consisting of non-overlapping, star-shaped elements K ∈ K h such that Figure 2: Auxiliary triangulation T h (K) of star-shaped element K, altitude h a of triangle T E ∈ T h (K) perpendicular to E and angle β Their boundaries ∂K consist of nodes and edges, cf. Figure 1. An edge E = z b z e is always located between two nodes, the one at the beginning z b and the one at the end z e . We set N (E) = {z b , z e }, these points are fixed once per edge and they are the only nodes on E. The sets of all nodes and edges in the mesh are denoted by N h and E h , respectively. The sets of nodes and edges corresponding to an element K ∈ K h are abbreviated to N (K) and E(K). Later on, the sets of nodes and edges which are in the interior of Ω, on the Dirichlet boundary Γ D and on the Neumann boundary Γ N are needed. Therefore, we decompose N h and E h into The length of an edge E and the diameter of an element K are denoted by h E and h K = sup{|x − y| : x, y ∈ ∂K}, respectively.

Polygonal mesh
In order to proof convergence and approximation estimates the meshes have to satisfy some regularity assumptions. We recall the definition gathered from [32] which is needed in the remainder of the presentation. Definition 1. The family of meshes K h is called regular if it fulfills: 1. Each element K ∈ K h is a star-shaped polygon with respect to a circle of radius ρ K and midpoint z K .
2. The aspect ratio is uniformly bounded from above by σ K , i.e.
3. There is a constant c K > 0 such that for all elements K ∈ K h and all its edges The circle in the definition is chosen in such a way that its radius is maximal, cf. Figure 1. If the position of the circle is not unique, its midpoint z K is fixed once per element. Without loss of generality, we assume h K < 1, K ∈ K h that can always be satisfied by scaling of the domain.
In the following, we give some useful properties of regular meshes. An important analytical tool is an auxiliary triangulation T h (K) of the elements K ∈ K h . This triangulation is constructed by connecting the nodes on the boundary of K with the point z K of Definition 1, see Figure 2. Consequently, T h (K) consists of the triangles T E for E = z b z e ∈ E(K), which are defined by the points z b , z e and z K . Lemma 1. Let K be a polygonal element of a regular mesh K h . The auxiliary triangulation T h (K) is shape-regular in the sense of Ciarlet [8], i.e. neighbouring triangles share either a common node or edge and the aspect ratio of each triangle is uniformly bounded by some constant σ T , which only depends on σ K and c K .
Proof. Let T E ∈ T h (K) be a triangle with diameter h T E and let ρ T E be the radius of the incircle. It is known that the area of T E is given by where h a is the altitude of the triangle perpendicular to E, see Figure 2. Since the element K is star-shaped with respect to a circle of radius ρ K , the line through the side E ∈ E h of the triangle does not intersect this circle. Thus, h a ≥ ρ K and we have the estimate |T E | ≥ 1 2 h E ρ K . Together with Definitions 1, we obtain In the previous proof, we discovered and applied the estimate for the area of the auxiliary triangle. This inequality will be of importance once more. We may also consider the auxiliary triangulation T h (Ω) of the whole domain Ω which is constructed by gluing the local triangulations T h (K). Obviously, T h (Ω) is also shape-regular in the sense of Ciarlet. Furthermore, the angles in the auxiliary triangulation T h (K) next to ∂K can be bounded from below. This gives rise to the following result.
Lemma 2. Let K h be a regular mesh. Then, there is an angle α K with 0 < α K ≤ π/3 such that for all elements K ∈ K h and all its edges E ∈ E(K) the isosceles triangle T iso E with longest side E and two interior angles α K lies inside T E ∈ T h (K) and thus inside the element K, see Figure 2. The angle α K only depends on σ K and c K .
Proof. Let T E ∈ T h (K). We bound the angle β in T E next to E = z b z e from below, see Figure 2. Without loss of generality, we assume that β < π/2. Using the projection y of z K onto the straight line through the edge E, we recognize Consequently, it is β ≥ arcsin σ −1 K . Since this estimate is valid for all angles next to ∂K of the auxiliary triangulation, the isosceles triangles T iso E , E ∈ E(K) with common angle α K = min{π/3, arcsin σ −1 K } lie inside the auxiliary triangles T E and therefore inside K.
for z ∈ N h , E ∈ E h and K ∈ K h , see Figure 3. An important role play the neighbourhoods ω z and ω z of a node z. Their diameters are denoted by h ωz and h ωz . Furthermore, ω z is of comparable size to K ⊂ ω z as shown in Lemma 3. Let K h be regular. Then, the mesh fulfils: 1. The number of nodes and edges per element is uniformly bounded, i.e. |N (K)| = |E(K)| ≤ c, ∀K ∈ K h .
3. Each element is covered by a uniformly bounded number of neighbourhoods of elements, i.e. |{K ∈ K h : The generic constant c > 0 only depends on σ K and c K from Definition 1.
Proof. The properties 2, 3 and 4 are proven as in [31] and rely on Lemma 2, which guaranties that all interior angles of the polygonal elements are bounded from below by some uniform angle α K .
To see the first property, we exploit the regularity of the mesh. Let K ∈ K h . In 2D it is obviously |N (K)| = |E(K)|. With the help of (3), we obtain

Trefftz-like trial space
In order to define a conforming, discrete space V k h ⊂ V over the polygonal decomposition K h , we make use of local Trefftz-like trial functions. We give a brief review of the approach in [32], which yields an approximation space of order k ∈ N. The discrete space is constructed by prescribing its basis functions that are subdivided into nodal, edge and element basis functions. Each of them is locally (element-wise) the unique solution of a boundary value problem.
The nodal functions ψ z , z ∈ N h are defined by ψ z is linear on each edge of the mesh.
Let p E,0 = ψ z b | E and p E,1 = ψ ze | E for E = z b z e and let {p E,i : i = 0, . . . , k} be a basis of the polynomial space P k (E) of order k over E. Then the edge bubble functions Finally, the element bubble functions ψ K,i,j for i = 0, . . . , k − 2 and j = 0, . . . , i, where {p K,i,j : i = 0, . . . , k − 2 and j = 0, . . . , i} is a basis of the polynomial space P k−2 (K) of order k − 2 over K.
Due to the regularity of the local problems each basis function ψ fulfills ψ ∈ C 2 (K) ∩ C 0 (K) for convex K ∈ K h . In the case of non-convex elements, the definitions are understood in the weak sense, but we still have ψ ∈ H 1 (K) ∩ C 0 (K). Consequently, it is ψ ∈ H 1 (Ω) and we set contains the locally (weakly) harmonic trial functions and contains the trial functions vanishing on all edges. It can be seen that the restriction of V k h to an element K ∈ K h , which does not touch the Dirichlet boundary, fulfills

Discrete variational formulation
With the help of the conforming approximation space V k h ⊂ V over polygonal meshes, we can give the discrete version of the variational formulation (2). Thus, we obtain For simplicity, we assume that g D ∈ P k pw (∂Ω), such that its extension into H 1 (Ω) can be chosen as interpolation using nodal and edge basis functions.
A closer look at the trial functions enables a simplification of (6). Since the nodal and edge basis functions ψ are (weakly) harmonic on each element, they fulfill and especially (∇ψ, Since the diffusion coefficient is assumed to be piecewise constant, the variational formulation decouples. Thus, the discrete problem (6) and Moreover, (9) reduces to a problem on each element since the support of the element bubble functions are restricted to one element. Consequently, u h,2 can be computed in a preprocessing step on an element level. A further observation is, that in the case of vanishing source term f , equation (9) yields u h,2 = 0. Therefore, it is sufficient to consider the discrete variational formulation (8) for the homogeneous diffusion equation.
In [32], it has been shown that the discrete variational formulation with trial space V k h yields optimal rates of convergence for uniform mesh refinement under classical regularity assumptions on the boundary value problem. More precisely, we have and additionally, for H 2 -regular problems, where h = max{h K : K ∈ K h } and the constant c only depends on the mesh parameters given in Definition 1 as well as on k.

BEM-based approximation
In this subsection we briefly discuss how to deal with the implicitly defined trial functions. More details can be found in the publications [25,32], for example. We focus on a local problem of the following form where p K ∈ P k−2 (K), p ∂K ∈ P k pw (∂K) and K ∈ K h is an arbitrary element. Let γ K 0 : H 1 (K) → H 1/2 (∂K) be the trace operator and denote by γ K 1 the conormal derivative, which maps the solution of (10) to its Neumann trace on the boundary ∂K, see [21]. For v ∈ H 1 (K) with ∆v in the dual of H 1 (K), γ K 1 v is defined as unique function in H −1/2 (∂K) which satisfies Greens first identity such that where n K (x) denotes the outer normal vector of the domain K at x. Without loss of generality, we only consider the Laplace equation, i.e. p K = 0 in (10). In the case of a general ψ ∈ V k h with −∆ψ = p K , we can always construct a polynomial q ∈ P k (K), see [19], such that −∆(ψ −q) = 0 in K and ψ −q ∈ P k pw (∂K) on ∂K. Due to the linearity of γ K 1 , it is with γ K 1 q ∈ P k−1 pw (∂K). Thus, γ K 1 ψ can be expressed as conormal derivative of a (weakly) harmonic function and a piecewise polynomial. Consequently, we reduced the general case to the Laplace problem.
It is well known, that the solution of (10) with p K = 0 has the representation where U * is the fundamental solution of minus the Laplacian with see, e.g., [21]. The Dirichlet trace γ K 0 ψ = p ∂K is known by the problem, whereas the Neumann trace γ K 1 ψ is an unknown quantity. Taking the trace of (11) yields the boundary integral equation with the single-layer potential operator and the double-layer potential operator To obtain an approximation of the Neumann trace a Galerkin scheme is applied to (12), which only involves integration over the boundary of the element, see, e.g., [27]. Thus, we approximate γ K 1 ψ by γ K 1 ψ ∈ P k−1 pw,d (∂K) such that where The discrete formulation (13) has a unique solution. The approximation described here is very rough, since no further discretization of the edges is performed. For k = 1, the Neumann trace is approximated by a constant on each edge E ∈ E(K), for example. The resulting boundary element matrices are dense, but, they are rather small, since the number of edges per element is bounded by Lemma 3. Consequently, the solution of the system of linear equations is approximated by an efficient LAPACK routine. The boundary element matrices of different elements are independent of each other. Therefore, they are computed in parallel during a preprocessing step for the overall simulation. Furthermore, once computed, they are used throughout the computations.

Approximated discrete variational formulation
It remains to discuss the approximations of the terms b(ψ, ϕ), (f, ϕ) Ω and (g N , ϕ) Γ N for ψ, ϕ ∈ V k h in the variational formulation (8)- (9). For the approximation of the first term b(ψ, ϕ), we apply Greens identity locally.
The approximation of this bilinear form is defined by where γ K 1 ψ ∈ P k−1 pw,d (∂K) is the approximation coming from the BEM and v is an appropriate polynomial approximation of v over K as, i.e., the averaged Taylor polynomial in Lemma (4.3.8) of [6]. If b h (ψ, ϕ) is evaluated for ψ, ϕ ∈ V k h,1 , we obviously end up with only boundary integrals, where polynomials are integrated. For ψ, ϕ ∈ V k h,2 , the volume integrals remain. These integrals can also be evaluated analytically, since the integrand is a polynomial of order 2k − 2.
The volume integral (f, ϕ) Ω is treated in a similar way as (∆ψ, v) K . Let ϕ be a piecewise polynomial approximation of ϕ, then we use as approximation (f, ϕ) Ω . The resulting integral is treated over each polygonal element by numerical quadrature. It is possible to utilize an auxiliary triangulation for the numerical integration over the elements or to apply a quadrature scheme for polygonal domains directly, see [22].
The last integral (g N , ϕ) Γ N is treated by means of Gaussian quadrature over each edge within the Neumann boundary. The data g N is given and the trace of the shape functions is known explicitly by their definition.
Finally, the approximated discrete variational formulation for u h = u h,1 + u h,2 ∈ V k h is given as and Remark 1. The approximation of (∆ψ, v) K by (∆ψ, v) K and (f, v) K by (f, v) K with a polynomial v can be interpreted as numerical quadrature. In consequence, we can directly approximate (∆ψ, v) K by an appropriate quadrature rule in the computational realization without the explicit construction of v. In this case the representation formula (11) is used to evaluate the shape functions inside the elements.
Remark 2. The representation (15) is advantageous for the a posteriori error analysis in the following, where we indeed consider b h (ψ, v) for v ∈ V . In the computational realization, however, we are only interested in v = ϕ ∈ V k h . Consequently, one might use for ψ, ϕ ∈ V k h,2 with supp (ψ) = supp (ϕ) = K homogenization. Therefore, let q ∈ P k (K) such that ∆ϕ = ∆q. This time, Greens identity yields since (∇(ϕ − q), ∇ψ) K = 0 due to (7). In this formulation, the volume integral can be evaluated analytically and there is no need to approximate it. But, we have used the properties of the basis functions. For the approximation of the boundary integral (γ K 1 ψ, γ K 0 q) ∂K one might even use an improved strategy compared to the one described above, see [25].

Quasi-interpolation
In the case of smooth functions like in H 2 (Ω), it is possible to use nodal interpolation. Such interpolation operators are constructed and studied in [32], and they yield optimal orders of approximation. The goal of this section, however, is to define interpolation for general functions in H 1 (Ω). Consequently, quasi-interpolation operators are applied, which utilizes the neighbourhoods ω z and ω z defined in (4) and (5). Let ω ⊂ Ω be a simply connected domain. With the help of the L 2 -projection Q ω : L 2 (ω) → R into the space of constants, we introduce the quasi-interpolation for v ∈ V . The definition of I h is a direct generalization of the quasi-interpolation operator in [31] and I h is a small modification. The definition is very similar to the one of Clément [9]. The major difference is the use of non-polynomial trial functions on polygonal meshes with non-convex elements. The quasi-interpolation operator I h has been studied in [31] on polygonal meshes with convex elements. In the following, the theory is generalized to meshes with star-shaped elements, which are in general non-convex. Before we give approximation results for the quasi-interpolation operators, some auxiliary results are reviewed and extended. If no confusion arises, we write v for both the function and the trace of the function on an edge. Since Lemma 2 guaranties the existence of the isosceles triangles with common angles for non-convex elements in a regular mesh, we can directly use the following lemma proven in [31].
with the isosceles triangle T iso E ⊂ K from Lemma 2, where c only depends on α K , and thus, on the regularity parameters σ K and c K .
Next, we prove an approximation property for the L 2 -projection. Lemma 5. Let K h be a regular mesh. There exist uniform constants c, which only depend on the regularity parameters σ K and c K , such that for every z ∈ N h , it holds v − Q ωz v 0,ωz ≤ ch ωz |v| 1,ωz for v ∈ H 1 (ω z ), where h ωz and h ωz denote the diameter of ω z and ω z , respectively.
This result is of interest on its own. It is known that these inequalities hold with the Poincaré constant which depends on the shape of ω, see [28]. For convex ω, the authors of [23] showed C P (ω) < 1/π. In our situation, however, ω z is a patch of non-convex elements which is itself non-convex in general. The lemma says, that we can bound C(ω z ) and C( ω z ) by a uniform constant c depending only on the regularity parameters of the mesh. The main tool in the forthcoming proof is Proposition 2.10 (Decomposition) of [28]. As preliminary of this proposition, an admissible decomposition {ω i } n i=1 of ω with pairwise disjoint domains ω i and ω = n i=1 ω i is needed. Admissible means in this context, that there exist triangles {T i } n i=1 such that T i ⊂ ω i and for every pair i, j of different indices, there is a sequence i = k 0 , . . . , k = j of indices such that for every m the triangles T k m−1 and T km share a complete side. Under these assumptions, the Poincaré constant of ω is bounded by Proof (Lemma 5). The second estimate in the lemma can be seen easily. The neighbourhood ω z is a patch of triangles, see (5). Thus, we choose : z ∈ T } for the admissible decomposition of ω z . According to Lemma 3, n is uniformly bounded. Furthermore, it is C P (ω i ) < 1/π, | ω z | < h 2 ωz and h 2 ω i /|T i | = h 2 T i /|T i | ≤ c, because of the shape-regularity of the auxiliary triangulation proven in Lemma 1. Consequently, applying Proposition 2.10 (Decomposition) of [28] yields C P ( ω z ) < c that proves the second estimate in the lemma.
The first estimate is proven in three steps. First, consider an element K, which fulfills the regularity assumptions of Definition 1. Remembering the construction of T h (K), K can be interpreted as patch of triangles ω z K corresponding to the point z K . Arguing as above gives C P (K) < c.
Next, we consider a neighbourhood ω z and apply Proposition 2.10 of [28]. In the second step, to simplify the explanations, we assume that the patch consists of only one element, i.e. ω z = K ∈ K h , and let E 1 , E 2 ∈ E(K) with z = E 1 ∩ E 2 . We decompose ω z , or equivalently K, into ω 1 and ω 2 such that n = 2. The decomposition is done by splitting K along the polygonal chain through the points z, z K and z , where z ∈ N (K) is chosen such that the angle β = ∠zz K z is maximized, see Figure 4 left. It is β ∈ (π/2, π], since K is star-shaped with respect to a circle centered at z K . The triangles {T i } n i=1 are chosen from the auxiliary triangulation in Lemma 1 as Figure 4: Construction of admissible decomposition for K and ω z from Figure 3 The element K is star-shaped with respect to a circle of radius ρ K and we have split this circle into two circular sectors during the construction of ω i , i = 1, 2. A small calculation shows that ω i is also star-shaped with respect to a circle of radius which lies inside the mentioned circular sector and fulfills ρ K /(1+ √ 2) < ρ ω i ≤ ρ K /2, see Figure 4 (left). Thus, the aspect ratio of ω i is uniformly bounded, since Furthermore, we observe that h ω i ≤ h K ≤ σ K ρ K ≤ σ K |zz K | and accordingly h ω i ≤ σ K |z z K |. Consequently, ω i , i = 1, 2 is a regular element in the sense of Definition 1 and, thus, we have already proved that C P (ω i ) ≤ c. Additionally, we obtain by (3) and by the regularity of the mesh that This yields together with |ω z | ≤ h 2 ωz and Proposition 2.10 (Decomposition) of [28] that and thus, a uniform bound in the case of ω z = K and n = 2.
In the third step, the general case, the patch ω z is a union of several elements, see (4) and Figure 3. In this situation, we proceed the construction for all neighbouring elements of the node z as in the second step, see Figure 4 (right). Consequently, n is two times the number of neighbouring elements. This number is uniformly bounded according to Lemma 3. The resulting decomposition {ω i } n i=1 is admissible with n i=1 T i = ω z and the estimate of [28] yields C P (ω z ) ≤ c, where c only depends on σ K and c K .
Finally, we formulate approximation properties for the quasi-interpolation operators I h and I h . The proof is skipped since it is analogous to [31] with the exception, that the generalized Lemmata 2, 3, 4 and 5 are applied instead of their counterparts.
where the constants c > 0 only depend on the regularity parameters σ K and c K , see Definition 1.

Residual based error estimate
In this section, we formulate the main results for the residual based error estimate and prove its reliability and efficiency. This a posteriori error estimate bounds the difference of the exact solution and the Galerkin approximation in the energy norm · b associated to the bilinear form, i.e. · 2 b = b(·, ·). Among others, the estimate contains the jumps of the conormal derivatives over the element edges. Such a jump over an internal edge E ∈ E h,Ω is defined by where K, K ∈ K h are the neighbouring elements of E with E ∈ E(K) ∩ E(K ). The element residual is given by

and the edge residual by
Theorem 1 (Reliability). Let K h be a regular mesh. Furthermore, let u ∈ g D + V and u h ∈ g D + V k h be the solutions of (2) and (16)- (17), respectively. Then the residual based error estimate is reliable, i.e.
where the error indicators are defined by The constant c > 0 only depends on the regularity parameters σ K , c K , see Definition 1, the approximation order k and on the diffusion coefficient a.
The term δ K measures the approximation error in the Neumann traces of the basis functions of V k h coming from the boundary element method. To state the efficiency, we introduce the notation · b,ω for ω ⊂ Ω, which means that the energy norm is only computed over the subset ω. More precisely, it is v 2 b,ω = (a∇v, ∇v) ω for our model problem.
Theorem 2 (Efficiency). Under the assumptions of Theorem 1, the residual based error indicator is efficient, i.e.
, where f and g N are piecewise polynomial approximations of the data f and g N , respectively. The constant c > 0 only depends on the regularity parameters σ K , c K , see Definition 1, the approximation order k and on the diffusion coefficient a.
The terms involving the data approximation f − f 0, ω K and g N − g N 0,E are often called data oscillations. They are usually of higher order. Furthermore, the approximation of the Neumann traces by the boundary element method appear in the right hand side. This term is related to δ K .

Remark 3.
Under certain conditions on the diffusion coefficient it is possible to get the estimates in Theorems 1 and 2 robust with respect to a, see, i.e., [24].

Reliability
We follow the classical lines in the proof of the reliability, see, i.e., [29]. However, we have to take care on the polygonal elements and the quasi-interpolation operators.
Proof (Theorem 1). The bilinear form b(·, ·) is a scalar product on V due to its boundedness and coercivity, and thus, V is a Hilbert space together with b(·, ·) and · b . The Riesz representation theorem yields To see the reliability of the residual based error estimate, we use this representation of the error in the energy norm and reformulate and estimate the term |R(v)| in the following. Let v h ∈ V 1 h , after a few manipulations using (2) and (16)- (17), we obtain The formulas (14) and (15) for b(·, ·) and b h (·, ·) lead to After some rearrangements of the sums, we obtain The Cauchy-Schwarz and the triangular inequality yield We have introduced the notation v for the polynomial approximation of v. Consequently, we obtain for the approximation v h of v h according to Lemma (4.3.8) in [6] v with a constant c only depending on σ K and k. Let I h,C be the usual Clément interpolation operator over the auxiliary triangulation T h (Ω), which maps into the space of piecewise linear and globally continuous functions, see [9]. Due to the construction of I h , it is harmonic on K, it minimizes the energy such that where an interpolation error estimate for the Clément operator has been applied, see [9]. Proposition 1 and |N (K)| < c, according to Lemma 3 Estimating v − v h 0,E with Proposition 1 gives where we have used again Cauchy-Schwarz and Lemma 3. We combine the previous estimates and apply the trace inequality in such a way that the residual in (19) is bounded by The last estimate is valid since each element is covered by a finite number of patches, see Lemma 3. The norm · 1,Ω and the semi-norm | · | 1,Ω are equivalent on V and a/a min > 1. Therefore, we conclude which finishes the proof.

Efficiency
We adapt the bubble function technique to polygonal meshes. Therefore, let φ T and φ E be the usual polynomial bubble functions over the auxiliary triangulation T h (Ω), see [2,29]. Here, φ T is a quadratic polynomial over the triangle T ∈ T h (Ω), which vanishes on Ω \ T and in particular on ∂T . The edge bubble φ E is a piecewise quadratic polynomial over the adjacent triangles in T h (Ω), sharing the common edge E, and it vanishes elsewhere. We define the new bubble functions over the polygonal mesh as Lemma 6. Let K ∈ K h and E ∈ E(K). The bubble functions satisfy and fulfill for p ∈ P k (K) the estimates The constants c > 0 only depend on the regularity parameters σ K , c K and on the approximation order k.
Proof. Similar estimates are valid for φ T and φ E on triangular meshes, see [2,29]. By the use of Cauchy-Schwarz inequality and the properties of the auxiliary triangulation T h (Ω) the estimates translate to the new bubble functions. The details of the proof are omitted.
With these ingredients the proof of Theorem 2 can be addressed. The arguments follow the line of [2].
Proof (Theorem 2). Let R K ∈ P k (K) be a polynomial approximation of the element residual

Lemma 6 gives
We thus get , and by the lower triangular inequality Next, we consider the edge residual. Let R E ∈ P k (E) be an approximation of R E , with E ∈ E h,Ω . The case E ∈ E h,N is treated analogously.
Applying Lemma 6 and the previous formula leads to Therefore, it is By the lower triangular inequality, h −1 K ≤ h −1 E and the previous estimate for R K 0,K we obtain Let f and g N be piecewise polynomial approximations of f and g N , respectively. We Consequently, we have R K ∈ P k (K) and R E ∈ P k (E). Finally, the estimates for R K 0,K and R E 0,E yield after some applications of the Cauchy-Schwarz inequality and Lemma 3

Application on uniform meshes
The residual based error estimate can be used as stopping criteria to check if the desired accuracy is reached in a simulation on a sequence of meshes. However, it is well known that residual based estimators overestimate the true error a lot. But, because of the equivalence of the norms · 1,Ω and · b on V , we can still use η R to verify numerically the convergence rates for uniform mesh refinement when h → 0.

Application in adaptive FEM
The classical adaptive finite element strategy proceeds in the following steps Sometimes an additional COARSENING step is introduced. We have all ingredient to formulate an adaptive BEM-based FEM on polygonal meshes. In the SOLVE step, we approximate the solution of the boundary value problem on a given polygonal mesh in the high order trial space V k h with the help of the BEM-based FEM.
The ESTIMATE part is devoted to the computation of local error indicators which are used to gauge the approximation accuracy over each element. Here, we use the term η K from the residual based error estimate.
In MARK , we choose some elements according to their indicator η K for refinement. Several marking strategies are possible, we implemented the Dörflers marking, see [13].
Finally in REFINE , the marked elements are refined, and thus, a problem adapted mesh is generated. In the implementation each marked element is split into two new ones, as proposed in [31], and the regularity conditions of the mesh are checked. At this point, we like to stress that the polygonal meshes are beneficial in this context, since hanging nodes are naturally included. Consequently, there is no need for an additional effort avoiding hanging nodes or treating them as conditional degrees of freedom. Also a COARSENING of the mesh is very easy. This can be achieved by simply gluing elements together.

Numerical experiments
In the following we present numerical examples on uniform and adaptive refined meshes, see Figure 5. For the convergence analysis, we consider the error with respect to the mesh size h = max{h K : K ∈ K h } for uniform refinement. This makes no sense for adaptive strategies. Since the relation holds for the number of degrees of freedom (DoF) on uniform meshes, we study the convergence of the adaptive BEM-based FEM with respect to them. where f ∈ L 2 (Ω) is chosen in such a way that u(x) = sin(πx 1 ) sin(πx 2 ) for x ∈ Ω is the exact solution. The solution in smooth, and thus, we expect optimal rates of convergence for uniform mesh refinement. The problem is treated with the BEMbased FEM for different approximation orders k = 1, 2, 3 on a sequence of meshes with L-shaped elements of decreasing diameter, see Figure 5 left. In Figure 6, we give the convergence graphs in logarithmic scale for the value η R /|u| 1,Ω , which behaves like the relative H 1 -error, and the relative L 2 -error with respect to the mesh size h. The example confirms the theoretical convergence rates stated in Subsection 2.3.

Adaptive refinement strategy
Let Ω = (−1, 1) × (−1, 1) ⊂ R 2 be split into two domains, Ω 1 = Ω \ Ω 2 and Ω 2 = (0, 1) × (0, 1). Consider the boundary value problem −div (a∇u) = 0 in Ω, u = g on Γ, Using polar coordinates (r, ϕ), we choose the boundary data as restriction of the global function g(x) = r λ cos(λ(ϕ − π/4)) for x ∈ R 2 + , β cos(λ(π − |ϕ − π/4|)) else, This problem is constructed in such a way that u = g is the exact solution in Ω. Due to the ratio of the jumping coefficient it is u ∈ H 2 (Ω) with a singularity in the origin of the coordinate system. Consequently, uniform mesh refinement does not yield optimal rates of convergence. Since f = 0, it suffices to approximate the solution in V k h,1 with the variational formulation (16). We implemented an adaptive strategy according to Subsection 4.4, where the introduced residual based error indicator is utilized in the ESTIMATE step. Starting from an initial polygonal mesh, see Figure 5 middle, the adaptive BEM-based FEM produces a sequence of locally refined meshes. The approach detects the singularity in the origin of the coordinate system and polygonal elements appear naturally during the local refinement, see Figure 5 right. In Figure 7, the energy error u − u h b as well as the error estimator η R are plotted with respect to the number of degrees of freedom in logarithmic scale. As expected by the theory the residual based error estimate represents the behaviour of the energy error very well. Furthermore, the adaptive approach yields optimal rates of convergence in the presence of a singularity.

Conclusion
This publication contains one of the first results on a posteriori error control for conforming approximation methods on polygonal meshes. Such general meshes are very flexible and convenient especially in adaptive mesh refinement strategies. But, their potential application areas are not yet fully exploited. The presented results are a building block for future developments of efficient numerical methods involving polygonal discretizations.