Local convergence of the boundary element method on polyhedral domains

The local behavior of the lowest order boundary element method on quasi-uniform meshes for Symm’s integral equation and the stabilized hyper-singular integral equation on polygonal/polyhedral Lipschitz domains is analyzed. We prove local a priori estimates in \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}L2 for Symm’s integral equation and in \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}H1 for the hyper-singular equation. The local rate of convergence is limited by the local regularity of the sought solution and the sum of the rates given by the global regularity and additional regularity provided by the shift theorem for a dual problem.


Introduction
The boundary element method (BEM) for the discretization of boundary integral equations is an established numerical method for solving partial differential equations on (un)bounded domains. As an energy projection method, the Galerkin BEM is, like the finite element method (FEM), (quasi-)optimal in some global norm. However, often the quantity of interest is not the error on the whole domain, but rather a local error on part of the computational domain. For the FEM, the analysis of local errors goes back at least to [18]; advanced versions can be found in [10,30]. For (1.1) where u is the exact solution, u h the finite element approximation from a space X h of piecewise polynomials, and B 0 B 1 are open subsets of Ω with R:=dist(B 0 , ∂ B 1 ). Thus, the local error in the energy norm is bounded by the local best approximation on a larger domain and the error in the weaker L 2 -norm. The local best approximation allows for convergence rates limited only by the local regularity; the L 2 -error is typically controlled with a duality argument and limited by the regularity of the dual problem as well as the global regularity of the solution. Therefore, if the solution is smoother locally, we can expect better rates of convergence for the local error.
Significantly fewer works study the local behavior of the BEM. The case of smooth two dimensional curves is treated in [5,21,28], in [27] three dimensional screen problems are studied, and [14] discusses local error estimates on polygons. [19,20] provide estimates in the L ∞ -norm on smooth domains. Local error estimators for the BEM are presented in [23]. However, for the case of piecewise smooth geometries such as polygonal and polyhedral domains, sharp local error estimates that exploit the maximal (local) regularity of the solution are not available. Moreover, the analyses of [14,21,27,28] are tailored to the energy norm and do not provide optimal local estimates in stronger norms, whereas [5] imposes additional global regularity.
In this article, we obtain sharp local error estimates for lowest order discretizations on quasi-uniform meshes for Symm's integral equation in the L 2 -norm and for the (stabilized) hyper-singular integral equation in the H 1 -seminorm on polygonal/polyhedral domains. Structurally, the local estimates are similar to (1.1): The local error is bounded by a local best approximation error and a global error in a weaker norm. More precisely, our local convergence rates depend only on the local regularity and the sum of the rates given by the global regularity and the additional regularity of the dual problem on polygonal/polyhedral domains. Numerical examples show the sharpness of our analysis. As discussed in Remark 2.5 below, our results improve [21,27,28] as estimates in L 2 (for Symm's equation) and H 1 (for the hyper-singular equation) are obtained there from local energy norm estimates with the aid of inverse estimates, thereby leading to a loss of h −1/2 . In contrast, we avoid using an inverse inequality to go from the energy norm to a stronger norm.
The paper is structured as follows. We start with some notations and then present the main results for both Symm's integral equation and the hyper-singular integral equation in Sect. 2. In Sects. 3 and 4 we are concerned with the proofs of these results. First, some technical preliminaries that exploit the additional regularity on piecewise smooth geometries to prove some improved a priori estimates for solutions of Poisson's equation as well as for the boundary integral operators are presented. Then, we prove the main results, first for Symm's equation, then for the stabilized hyper-singular equation. In principle, the proofs take ideas from [30], but important modifications of the arguments are necessary due to the nonlocal character of the integral operators. As in [30] a key ingredient are interior regularity estimates, which were provided recently in [11,12], and to exploit some additional smoothing properties of commutators that arise in a localization step. Finally, Sect. 5 provides numerical examples that underline the sharpness of our theoretical local a priori estimates.

Notation on norms
For domains ω ⊂ R d , we define the integer order Sobolev spaces H k (ω), k ∈ N 0 , in the standard way [15, p. 73ff]. The fractional Sobolev spaces H k+s (ω), k ∈ N 0 , s ∈ (0, 1) are defined by the Slobodeckii norm as described in [15, p. 73ff]. For open sets ω = ∪ m i=1 ω i consisting of finitely many components of connectedness ω i , the Sobolev spaces H k+s (ω) are understood in a piecewise way with norm u 2 H k+s (ω) = i u 2 H k+s (ω i ) . The spaces H s (ω), s ≥ 0, consist of those function whose zero extension to R d is in H s (R d ). The spaces H −s (ω), s ≥ 0, are taken to be the dual space of H s (ω). We will make use of the fact that for bounded Lipschitz domains ω H s (ω) = H s (ω) ∀s ∈ [0, 1/2  [15,Thm. 3.37], [17,Lem. 2.6]) assert the existence of a continuous lifting L in the range 0 < s < 1 as well so that (1.3) is an equivalent norm for 0 < s < 1 as well.
(iii) For polygonal (in 2D) and polyhedral (in 3D) Lipschitz domains the spaces H s (Γ ) in the range s ∈ (1, 3/2) can be characterized alternatively as follows: Let Γ i , i = 1, . . . , N , be the affine pieces of Γ , which may be identified with an interval (for the 2D case) or a polygon (for the 3D case). Then The equivalence (1.5) gives rise to yet another norm equivalence for the space is a compatibility condition. More generally, for s > 3/2 similar, more complicated compatibility conditions can be formulated to describe the space H s (Γ ) in terms of piecewise Sobolev spaces.
We will also need local norms on the boundary. For an open subset Γ 0 ⊂ Γ and s ≥ 0, we define local negative norms by In the following, we write γ int 0 for the interior trace operator, i.e., the trace operator from the inside of the domain and γ ext 0 for the exterior trace operator. For the jump of the trace of a function u we use the notation [γ 0 u]:=γ ext 0 u − γ int 0 u. In order to shorten notation, we write γ 0 for the trace, if the interior and exterior trace are equal, i.e., [γ 0 u] = 0. We denote the interior and exterior conormal derivative by γ int 1 u:=γ int 0 ∇u · n, γ ext 1 u:=γ ext 0 ∇u · n, where n denotes the normal vector pointing into R d \Ω. The jump of the normal derivative across the boundary is defined by [∂ n u]:=γ ext 1 u − γ int 1 u, and we write ∂ n u for the normal derivative if [∂ n u] = 0.
We will call axis-parallel squares/cubes "boxes".
centered at the origin contains Ω. The parameter α D ∈ (0, 1/2) is such that for every ε ∈ (0, α D ] there is C ε > 0 such that the a priori bound Recall that the norms · H s (B R Ω (0)\Γ ) , s > 0 are understood as the sum of the norm on Ω and B R Ω (0)\Ω, i.e.,

Remark 2.2
The condition on the parameter α D in Assumption 2.1 can be described in terms of two Dirichlet problems, one posed on Ω and one posed on B R Ω (0)\Ω. For each of these two domains, a shift theorem is valid, and α D is determined by the more stringent of the two conditions. It is worth stressing that the type of boundary In the case d = 2 the parameter α D is determined by the extremal angles of the polygon Ω. Specifically, let 0 < ω j < 2π , j = 1, . . . , J , be the interior angles of the polygon Ω. Then, Assumption 2.1 is valid for any α D > 0 that satisfies (Note that ω j = π for all j so that the right inequality is indeed strict.) We consider Symm's integral equation in its weak form: Here, the single-layer operator V is given by where, with the surface measure |S d−1 | of the Euclidean sphere in R d , we set The single layer operator V is a bounded linear operator in L(H −1/2+s (Γ ), H 1/2+s (Γ )) for |s| ≤ 1 2 , [22,Thm. 3.1.16]. It is elliptic for s = 0 with the usual proviso for d = 2 that diam(Ω) < 1, which we may assume by scaling.
Let T h = {T 1 , . . . , T N } be a quasi-uniform, regular and γ -shape regular triangulation of the boundary Γ with mesh-width h := max T ∈T h diam(T ). By S 0,0 (T h ):={u ∈ L 2 (Γ ) : u| T j is constant ∀T j ∈ T h } we denote the space of piecewise constant functions on the mesh T h . The Galerkin formulation of (2.3) reads: The following theorem is one of the main results of this paper. It estimates the Galerkin error in the L 2 -norm on a subdomain by the local best approximation error in L 2 on a slightly larger subdomain and the global error in a weaker norm.

Theorem 2.3 Let Assumption 2.1 hold and let
The constant C > 0 depends only on Γ, Γ 0 , Γ , d, R, and the γ -shape regularity of T h .
If we additionally assume higher local regularity as well as some (low) global regularity of the solution φ, this local estimate implies that the local error converges faster than the global error, which is stated in the following corollary.
In the results of [18,30] singularities far from the domain of interest still have a weak influence on the local convergence of the FEM. Corollary 2.4 shows that this is similar in the BEM: The a priori estimate shows the effect of singularities of the solution (represented by α) and those induced by the geometry (represented by α D ) affect the local convergence.

Remark 2.5
In comparison to [27], Corollary 2.4 gives a better result for the rate of convergence of the local error in the case where the convergence is limited by the global error in the weaker norm. More precisely, for the case φ ∈ H 1/2 ( Γ ) ∩ L 2 (Γ ), [27] obtains the local rate of 1/2, which coincides with our local rate. However, if φ ∈ H 1 ( Γ ), we obtain rate 1 in the L 2 -norm, whereas the rate in [27] remains 1/2.

Remark 2.6
Even for smooth functions f , the solution φ of (2.3) is, in general, not better than H α (Γ ) with α = 1 2 + α D . Recall from Remark 2.2 that α D is determined by the mapping properties for both the interior and the exterior Dirichlet problem. A special situation therefore arises if Symm's integral equation is obtained from reformulating an interior (or exterior) Dirichlet problem. To be specific, consider again the case d = 2 of a polygon Ω with interior angles ω j , j = 1, . . . , J . We rewrite the boundary value problem −Δu = 0 in Ω with u| Γ = g as the integral equation for the unknown function φ = ∂ n u with the double layer operator K defined by Then, φ ∈ H α (Γ ) for any α with α < 1/2 + min j π ω j .

The hyper-singular integral equation
For the Neumann problem, we assume an extended shift theorem as well.
, is a bounded Lipschitz domain whose boundary consists of finitely many affine pieces (i.e., Ω is the intersection of finitely many halfspaces). R Ω > 0 is such that the open ball B R Ω (0) ⊂ R d contains Ω. The parameter α N ∈ (0, α D ], where α D is the parameter from Assumption 2.1, is such that for The condition on the parameter α N again can be described in terms of two problems, a pure Neumann problem posed in Ω, for which we need a compatibility condition, and a mixed Dirichlet-Neumann problem posed on B R Ω (0)\Ω, which is uniquely solvable without the need to impose a solvability condition for f, g.
The parameter α N again depends only on the geometry and the corners/edges that induce singularities. In fact, on polygonal domains, i.e., d = 2, α D = α N , see, e.g., [9].
Studying the inhomogeneous Neumann boundary value problem −Δu = 0, ∂ n u = g, leads to the boundary integral equation of finding ϕ ∈ H 1/2 (Γ ) such that W ϕ = f with f ∈ H −1/2 (Γ ) satisfying the compatibility condition f, 1 L 2 (Γ ) = 0, and the hyper-singular integral operator W ∈ L(H 1/2 (Γ ), H −1/2 (Γ )) defined by We additionally assume that Γ is connected, so that the hyper-singular integral operator has a kernel of dimension one consisting of the constant functions. Therefore, the boundary integral equation is not uniquely solvable. Employing the constraint ϕ, 1 L 2 (Γ ) = 0 leads to the stabilized variational formulation which has a unique solution ϕ ∈ H 1/2 (Γ ), see, e.g., [26]. For the Galerkin discretization we employ lowest order test and trial functions in The following theorem is the analog of Theorem 2.3 for the hyper-singular integral equation. The local error in the H 1 -seminorm is estimated by the local best approximation error and the global error in a weak norm. Theorem 2.8 Let Assumption 2.7 hold and let T h be a quasi-uniform, γ -shape regular triangulation. Let ϕ ∈ H 1/2 (Γ ) and ϕ h ∈ S 1,1 (T h ) satisfy the Galerkin orthogonality condition 12 with a fixed constant C α N depending only on α N . Assume that ϕ ∈ H 1 ( Γ ). Then, we have The constant C > 0 depends only on Γ, Γ 0 , Γ , d, R, and the γ -shape regularity of T h .
Again, assuming additional regularity, the local estimate of Theorem 2.8 leads to a higher rate of local convergence of the BEM for the stabilized hyper-singular integral equation.

Shift theorems
The following two sections are dedicated to the proofs of Theorem 2.3 and Corollary 2.4 for Symm's integral equation as well as Theorem 2.8 and Corollary 2.9 for the hyper-singular integral equation. We start with some technical results that are direct consequences of the assumed shift theorems from Assumption 2.1 for the Dirichlet problem and Assumption 2.7 for the Neumann problem. The shift theorem of Assumption 2.1 implies the following shift theorem for Dirichlet problems: (i) There is a constant C > 0 depending only on Ω and α D such that Here, the constant C > 0 additionally depends on dist(B, ∂ B ).
Then, in view of (1.2), we have Integration by parts on Ω and B R Ω (0)\Ω and the boundary condition γ 0 v = 0 lead to We split the polygonal/polyhedral boundary Γ = m =1 Γ into its (smooth) faces Γ and prolong each face Γ to the hyperplane Γ ∞ , which decomposes R d into two half spaces Ω ± . Let χ ∈ L 2 (Γ ) be the characteristic function for Γ . Since the normal vector on a face does not change, we may use the trace estimate (note: 0 < α D < 1/2) facewise, to estimate As the boundary ∂ B R Ω (0) is smooth, standard elliptic regularity yields . This leads to Proof of (ii): With the lifting operator L : With the shift theorem from Assumption 2.1 we get which proves the second statement.
The following lemma collects mapping properties of the single-layer operator V that exploits the present setting of piecewise smooth geometries: Proof: Proof of (i): The case s ∈ (−1/2, 1/2) is shown in [22, Thm. 3.1.16], and for s = − 1 2 we refer to [29]. For the case s ∈ [1/2, 1), we exploit that Γ is piecewise smooth. We split the polygonal/polyhedral boundary Γ = m =1 Γ into its (smooth) faces Γ . Let χ ∈ L 2 (Γ ) be the characteristic function for Γ . Then, for We prolong each face Γ to the hyperplane Γ ∞ , which decomposes R d into two half spaces Ω ± . Due to s < 1, we have χ ϕ ∈ H −1/2+s (Γ ∞ ). Since the half spaces Ω ± have smooth boundaries, we may use the mapping properties of V on smooth geometries, see, e.g., [15,Thm. 6.13] to estimate Proof of (ii): In addition to the single layer operator V , we will need to understand localized versions of these operators, i.e., the properties of commutators. For a smooth cut-off function η, we define the commutators Since the singularity of the Green's function at x = y is smoothed by η(x) − η(y), we expect that the commutators C η , C η η have better mapping properties than the singlelayer operator; this is stated in the following lemma.
, Ω, and s. Furthermore, the operator is skew-symmetric (with respect to the L 2 (Γ )-inner product). (ii) The commutator C η η is a symmetric and continuous mapping C Here, the continuity constant depends only on η W 1,∞ (R d ) , Ω, and the constants appearing in Assumption 2.1.
Proof: Proof of (i): 1. step: We show the boundedness for the case 0 < s < 1/2. Let φ ∈ H −1+s (Γ ), and set Since the volume potential V φ is harmonic and in view of the jump relations The definition of C η and the definition of the norm · H 1+s (Γ ) from (1.3) prove the mapping properties of C η for 0 < s < 1/2. The mapping properties of the Newton potential ( see, e.g., [22,Thm. 3.1.2]) also lead to (3.10) With the mapping property C η : Proof of (ii): Again, the function v and the Newton potential of the right-hand side of (3.12) have the same decay for |x| → ∞, and the mapping properties of the Newton potential as well as the previous estimate (3.10) for (3.14) We apply Lemma 3.
, and we can estimate this term by an arbitrary negative norm of φ on Γ to obtain The mapping properties of V of Lemma 3.2, (ii) and the symmetry of V imply Inserting this in (3.14) together with the definition of the H 1+α D (Γ )-norm in (1.3), proves the lemma.
The shift theorem for the Neumann problem from Assumption 2.7 implies the following shift theorem.

Lemma 3.4 Let Assumption 2.7 be valid, and let u be the solution of the inhomogeneous problems
(i) There is a constant C > 0 depending only on Ω and α N such that (3.16) Here, the constant C > 0 depends on Ω, α N , and dist(B, ∂ B ).
Proof: Proof of (i): Integration by parts on Ω and B R Ω (0)\Ω and the boundary conditions of v lead to The definition of the norm (1.3) implies and the same estimate holds for γ ext

This leads to
The shift theorem from Assumption 2.7 and the trace inequality which proves the second statement.
The following lemma collects mapping properties of the double-layer operator K and the hyper-singular operator W that exploit the present setting of piecewise smooth geometries: Proof: Proof of (i): With the mapping properties of the single layer potential V ∈ which finishes the proof for the case s ∈ (1/2, For a smooth function η, we define the commutators By the mapping properties of W , both operators map H 1/2 (Γ ) → H −1/2 (Γ ). However, C η is in fact an operator of order 0 and C η η is an operator of positive order: The constant C depends only on η W 1,∞ (R d ) , Ω, and s. Furthermore, the operator is skew-symmetric (with respect to the extended L 2 -inner product). (ii) The commutator C η η is a symmetric and continuous mapping C , Ω, and the constants appearing in Assumption 2.7.
Proof: Proof of (i): 1. step: We show (3.20) for the range 0 < s < 1/2. For ϕ ∈ H 1/2 (Γ ), consider the operator with the single layer potential V and the double layer potential K from (3.17). Using the jump conditions The decay of u -the dominant part is the single-layer potential -and the Newton potential N (2∇η · ∇ K ϕ + (Δη) K ϕ) for |x| → ∞ are the same, which allows us to The trace estimate applied facewise as in the proof of Lemma 3.1, and estimates (3.3), Similarly, we obtain with Lemma 3.2, (i) 3. step: The operator C η is skew-symmetric: The operator W maps H 1/2 (Γ ) → H −1/2 (Γ ) and is symmetric. The skew-symmetry of C η then follows from a direct calculation.

step:
The skew-symmetry of C η allows us to extend (in a unique way) the operator as an operator H −s (Γ ) → H −s (Γ ) for 0 < s < 1/2 by the following argument: For ϕ, ψ ∈ H 1/2 (Γ ) we compute Proof of (ii): Let ϕ ∈ H −α N (Γ ). The argument leading to the first inequality in (3.22), i.e., the mapping properties of N , also shows Again, the decay of v and the Newton potential applied to the right-hand side of the equation are the same. Hence, the mapping properties of the Newton potential provide We apply Lemma 3.4 to K ϕ − K ϕ, with K ϕ:= 1 |Ω| K ϕ, 1 L 2 (Ω) . Since we assumed dist(Γ, ∂ B R Ω (0)) > 0, we have that K ϕ is smooth on ∂ B R Ω (0), and we can estimate this term by an arbitrary negative norm of ϕ on Γ to obtain The mean value K ϕ can be estimated with r 2 = |x| 2 , the observation Δr 2 = 2d, and integration by parts by where the last step follows since K is a bounded operator mapping H α N (Γ ) → H α N (Γ ) by Lemma 3.2. Using the mapping properties of W of Lemma 3.5, (iii) and inserting (3.29) in (3.28) leads together with a facewise trace estimate to The computation the mapping properties of V and the commutator of K (as normal trace of the commutator C η from Lemma 3.3, cf. (3.9)) prove the lemma.

Proof of main results
With the consequences of the shift theorems from the previous section, we can prove our main results, the local error estimates for Symm's integral equation and the hypersingular integral equation.

Symm's integral equation (proof of Theorem 2.3)
The main tools in our proofs are the Galerkin orthogonality and a Caccioppoli-type estimate for discrete harmonic functions that satisfy the orthogonality  As a consequence of this interior regularity estimate and Lemma 3.1, we get an estimate for the jump of the normal derivative of a discrete harmonic single-layer potential. (4.5) The constant C > 0 depends only on Ω, d, d, the γ -shape regularity of T h , η W 1,∞ (R d ) , and the constants appearing in Assumption 2.1.

Proof:
We split the function u = u far + u near , where the near field u near and the far field u far solve the Dirichlet problems We first consider γ int 1 u near -the case γ ext 1 u near is treated analogously. Let η be another cut-off function satisfying η ≡ 1 on Γ and supp η ⊂ B. The multiplicative trace inequality, see, e.g., [16,Thm. A.2], implies for any ε ≤ 1/2 that Since u near ∈ H h (B ), we use the interior regularity estimate (4.2) for the first term on the right-hand side of (4.6), and the second term of (4.6) can be estimated using (3.2). In total, we get for ε ≤ α D < 1/2 that (4.7) Let (4.8) With the classical a priori estimate for the inhomogeneous Dirichlet problem in the H 1 -norm, the commutator C η , and Lemma 3.3, we estimate We apply (3.1), (since η ≡ 0 on ∂ B R Ω (0) only the boundary terms for Γ appear) together with Young's inequality ab ≤ a p / p + b q /q applied with p = (1 + 2ε)/2ε, q = 1 + 2ε to obtain (4.8) ηV ζ h 2ε/(1+2ε) Similarly, we get for the second term in (4.7) Inserting everything in (4.7) and using h 1 gives Applying the same argument for the exterior Dirichlet boundary value problem leads to an estimate for the jump of the normal derivative It remains to estimate the far field u far , which can be treated similarly to the near field using a trace estimate and Lemma 3.1. Applying Lemma 3.1 with a cut-off function η satisfying η ≡ 1 on B ∩ Γ and supp η ⊂ B , the boundary term in (3.2) disappears since η(1 − η) ≡ 0, which simplifies the arguments: which proves the lemma.
We use the Galerkin projection Π V : H −1/2 (Γ ) → S 0,0 (T h ), which is defined by We denote by I h the L 2 (Γ )-orthogonal projection given by The operator I h has the following super-approximation property, [18]: For any discrete function ψ h ∈ S 0,0 (T h ) and a cut-off function η, we have (with implied constants depending on η W 1,∞ ) The following lemma provides an estimate for the local Galerkin error and includes the key steps to the proof of Theorem 2.3.

Lemma 4.3 Let the assumptions of Theorem 2.3 hold. Let Γ 0 be an open subset of
The constant C > 0 depends only on Γ, Γ 0 , d, R, and the γ -shape regularity of T h . With an inverse inequality and the L 2 -orthogonal projection I h , which satisfies the super-approximation property (4.12) for η 5 φ h , we get where the last estimate follows from Céa's lemma and super-approximation. The same argument leads to In fact, this argument shows L 2 -stability of Π V : The bounds (4.15), (4.16) together imply

(4.18)
For the first term on the right-hand side of (4.14), we want to use Lemma 4.2. Since [∂ n V ζ h ] = −ζ h ∈ S 0,0 (T h ) for any discrete function ζ h ∈ S 0,0 (T h ), we need to construct a discrete function satisfying the orthogonality condition (4.2). Using the Galerkin orthogonality with test functions ψ h with support supp ψ h ⊂ Γ 4 and noting that η 5 ≡ 1 on supp ψ h , we obtain with the commutator C η 5 defined in (3.5) Thus, defining we get on the volume box B 4 ⊂ R d a discrete harmonic function The correction ξ h can be estimated using the L 2 -stability (4.17) of the Galerkin projection, the mapping properties of V −1 , C η 5 , C The definition of ζ h , the bound (4.21), and the H −1/2 -stability of the Galerkin projection lead to With the L 2 -stability (4.17) of the Galerkin projection and (4.21) we get We use the orthogonality (4.19) satisfied by ζ h on Γ 4 , the L 2 -orthogonal projection I h and the properties of the commutator C η 5 given by Lemma 3.3 to arrive at Inserting (4.25)-(4.27) in (4.24) and using h 1, we arrive at Combining (4.14), (4.22) with (4.18), (4.23), (4.28), and finally (4.13), we get Since we only used the Galerkin orthogonality as a property of the error φ − φ h , we , and we have proven the inequality claimed in Lemma 4.3.
In order to prove Theorem 2.3, we need a lemma:

Lemma 4.4 For every δ > 0 there is a bounded linear operator J
with the following properties: Proof: Operators with such properties are obtained by the usual mollification procedure (on a length scale O(δ) for domains in R d ). This technique can be generalized to the present setting of surfaces with the aid of localization and charts. We also mention [1,7] where similar operators mapping into S 1,1 (T h ) are constructed.
We are in position to prove our main result, a local estimate for the Galerkinboundary element error for Symm's integral equation in the L 2 -norm.

Proof (of Theorem 2.3): Starting with Lemma 4.3, it remains to estimate the two terms
We start with the latter. Let η ∈ C ∞ (R d ) be a cut-off-function with η ≡ 1 on (1) such that the operator J ch of Lemma 4.4 has the support property supp J ch ( η) ⊂ B Γ 0 R/2+h . We will employ the operator I h • J ch : H −1 (Γ ) → S 0,0 (Γ ) with the L 2 -orthogonal projection I h . It is easy to see that we may assume that Concerning the approximation properties, we have With the definition of the commutators C η , C η η , the Galerkin orthogonality satisfied by e, and the fact that V : (4.31) The first term on the right-hand side of (4.31) can be treated in the same way as the term h α D /(1+2α D ) e L 2 ( Γ 1 ) on the right-hand side of Lemma 4.3, which is treated by iterating the L 2 -estimate of the statement of Theorem 4.3. That is, we set m: The assumption C α D h R ≤ 1 12 allows us to define m nested Since the term h α D /(1+2α D ) e L 2 ( Γ 1 ) again contains a local L 2 -norm, we may use Lemma 4.3 and (4.31) again on the larger set Γ 2 Γ to estimate Inserting this in the initial estimate of Lemma 4.3 (using h 1) leads to Now, the L 2 -term on the right-hand side is multiplied by h 2α D /(1+2α D ) , i.e., the square of the initial factor. Iterating this argument m − 2-times, provides the fac-tor h mα D /(1+2α D ) , and by the choice of m, we have h 1+α D ≤ h mα D /(1+2α D ) . Together with an inverse estimate, we obtain which proves the theorem.
where the second estimate is the standard global error estimate for the BEM, see [22]. It remains to estimate e H −1−α D (Γ ) , which is treated with a duality argument: We note that Assumption 2.1 and the jump relations imply the following shift theorem for Γ ) . Hence, with the Galerkin projection Π V , we estimate Therefore, the term of slowest convergence is of order O(h min{1/2+α+α D ,β} ), which proves the corollary.

Remark 4.5
The term of slowest convergence in the case of high local regularity is the global error in the negative H −1−α D (Γ )-norm, which is treated with a duality argument that uses the maximum amount of additional regularity on the polygonal/polyhedral domain. Therefore, further improvements of the convergence rate cannot be achieved with our method of proof. In fact, the numerical examples in the next section confirm the sharpness of this observation, i.e., that the best possible convergence is The trivial estimate ηe H −1/2 (Γ ) ηe L 2 (Γ ) immediately implies that the local convergence in the energy norm is at least of order O(h 1/2+α+α D ) as well. Again, analyzing the proof of Lemma 4.3, we observe that an improvement is impossible, since the limiting term is once more the error in the negative H −1−α D (Γ )-norm. Remark 4.6 Remark 4.5 states that the local rate of convergence is limited by the shift theorem of Assumption 2.1. If the geometry Ω is smooth, then elliptic shift theorems for the Dirichlet problem hold in a wider range, e.g., if f ∈ H 1/2 (Ω), we may get u ∈ H 5/2 (Ω). It can be checked that in this setting, an estimate of the form is possible since the commutator C η 5 η 5 in (4.21) then maps H −2 (Γ ) → H 1 (Γ ). If an even better shift theorem holds, then the H −2 -norm can be further weakened by using commutators of higher order. The best possible achievable local rates are then

The hyper-singular integral equation (proof of Theorem 2.8)
We start with the Galerkin orthogonality and a Caccioppoli-type estimate on D ⊂ R d for functions characterized by the orthogonality for some μ ∈ R. Here, we define the space of discrete harmonic functions H N h (D, μ) for an open set D ⊂ R d and μ ∈ R as  We use the Galerkin projection Π W : The following lemma collects approximation properties of the Galerkin projection that will be applied in both Lemmas 4.10 and 4.11 below. (4.36), and let η, η ∈

Lemma 4.8 Let Π W be the Galerkin projection defined in
The constant C > 0 depends only on Ω, the γ -shape regularity of T h , and η W 2,∞ (R d ) .
Proof: Let J h be a quasi-interpolation operator with approximation properties in the H s -seminorm, e.g., the Scott-Zhang-projection, [24]. We use super-approximation similarly to (4.12). Since ϕ h ∈ S 1,1 (T h ), we have to use the piecewise H 2 -norm, and an inverse inequality leads to where, in the last step, the assumption on η was used. Similarly, the H 1 -norm estimate ηϕ h − J h (ηϕ h ) H 1 (Γ ) h s ηϕ h H s (Γ ) holds. Interpolation finally leads to a superapproximation result in H s With an inverse inequality, see, e.g., [13,Thm. 3.2], as well as Céa's lemma this implies A similar argument leads to and consequently to the H 1 -stability of the Galerkin-projection.
In the following, we need stability and approximation properties of the Scott-Zhang projection J h in the space H 1+α N (Γ ) provided by the following lemma. [24]. Then, for s ∈ [0, 3/2) we have (4.39) and therefore, for every 0 ≤ t ≤ s < 3/2

Lemma 4.9 Let J h be the Scott-Zhang projection defined in
The constants C s , C s,t > 0 depend only on Ω, the γ -shape regularity of T h , and s, t.
Proof: We start with the proof of (4.39). The stability for the case s = 1 is given in [24] and the stability for the case s = 0 (note that Γ is a closed surface without boundary) is discussed in [3,Lemma 7]. By interpolation, (4.39) follows for 0 < s < 1. The starting point for the proof of (4.39) for s ∈ (1, 3/2) is that, by Remark 1.1, (iii), we may focus on a single affine piece Γ i of Γ and can exploit that the notion of H s (Γ i ) coincides with the standard notion on intervals (in 1D) and polygons (in 2D).
In particular, H s (Γ i ) can be defined as the interpolation space between H 1 (Γ i ) and It therefore suffices to show J h u H s (Γ i ) ≤ C u H s (Γ i ) . Since H s (Γ i ) is an interpolation space between H 1 (Γ i ) and H 2 (Γ i ), we can find (cf. [4]), for every t > 0, a function u t ∈ H 2 (Γ i ) with Let I h be an approximation operator with the simultaneous approximation property see, e.g., [4], [6,Thm. 14.4.2]. With an inverse inequality, cf. [8,Appendix], the H 1 -stability of the Scott-Zhang projection, and (4.41), (4.42), we estimate Choosing t = O(h), we get the H s (Γ i )-stability of J h and thus also the H s (Γ )stability of J h .
We only prove the approximation property (4.40) for s ∈ (1, 3/2) as the case s ∈ [0, 1] is covered by standard properties of the Scott-Zhang operator.
Case 1 ≤ t ≤ s < 3/2: we observe with the stability properties of J h and the approximation properties of I h

(4.44)
Case 0 < t < 1: The remaining cases are obtained with the aid of an interpolation inequality: which concludes the proof.
The following lemma is similar to Lemma 4.2. Here, we obtain an estimate for the jump of the trace of a discrete harmonic double-layer potential.
The constant C > 0 depends only on Ω, d, the γ -shape regularity of T h , and the constants appearing in Assumption 2.7.

Proof: Step 1 (Splitting into near and far-field):
Let η ∈ C ∞ 0 (R d ) satisfy η ≡ 1 on B ∩ Γ and supp η ⊂ B . Define the near-field u near and the far field u far as potentials u near : Here, z is a function with z ≡ μ on Γ ∩ B such that the compatibility condition ηW ζ h − ηz, 1 = (η − 1)W ζ h − ηz, 1 = 0 holds. Since W ζ h , 1 = 0 such a function exists. More precisely, we choose z ∈ L 2 (Γ ) to be the piecewise constant function The definition of z and η ≡ 1 on B ∩ Γ lead to Consequently, we obtain The last inequality follows from the orthogonality of W ζ h to discrete functions in S 1,1 (T h ) on B and the arguments shown in (4.47) below (specifically: go through the arguments of (4.47) with z ≡ μ).
Step 2 (Approximation of the near field): Let J h denote the Scott-Zhang projection. The ellipticity of W on H 1/2 (Γ )/R and the orthogonality (4.33) (4.47) With the same arguments and Lemma 4.9 we may estimate Together with the mapping properties of W from Lemma 3.5, v h , 1 = 0, the definition of v h , and the stability and approximation properties of J h from Lemma 4.9, we obtain (4.48), (4.46) h 1+αN W ζ h L 2 (Γ ) + |μ| . (4.49) With the mapping properties of W from Lemma 3.5, an inverse estimate, and (4.47) we obtain for 0 We first consider γ int 0 u near ; the case γ ext 0 u near is treated analogously. By construction of u near , we have 0). Let η be another cut-off function satisfying η ≡ 1 on Γ and supp η ⊂ B. The multiplicative trace inequality, see, e.g., [16,Thm. A.2], implies for any ε ≤ 1/2 that (4.52) Since u near ∈ H N h (B , 0), we may use the interior regularity estimate (4.35) with μ = 0 for the first term on the right-hand side of (4.52). The second factor of (4.52) can be estimated using (3.16) of Lemma 3.4. In total, we get for ε ≤ α N < 1/2 that The mapping properties of K imply with (4.47) and (4.50) (4.47) h 2ε/(1+2ε)+1/2 ζ h H 1 (Γ ) + |μ| , We apply (3.15) (note: u near has mean zero) and since . Together with (4.50), (4.49), and Young's inequality this leads to (3.15), (4.50) (4.49) Similarly, with (4.54) we get for the second term in (4.53) Inserting everything in (4.53) and choosing ε = α N gives Applying the same argument for the exterior trace leads to an estimate for the jump of the trace Step 3 (Approximation of the far field): We define the function ν ∈ H 1/2 (Γ ) as the solution of Let u far := K ν − K ν where K ν:= 1 |Ω| K ν, 1 L 2 (Ω) and η be another cut-off function with η ≡ 1 on Γ and supp η ⊂ B. Then, with the Galerkin projection Π W , the triangle inequality and the jump conditions of K imply The smoothness of K ν on ∂ B R Ω (0) and the coercivity of W on H 1/2 (Γ )/R lead to We apply Lemma 3.4 with a cut-off function η satisfying η ≡ 1 on B ∩ Γ and supp η ⊂ B . Then η ≡ 1 and z ≡ μ on B ∩ Γ imply η(1 − η) ≡ 0 and ηηz = ημ.
The H 1 -stability of the Galerkin projection from Lemma 4.8, a facewise trace estimate, and similar estimates as for the near field imply (4.46) ζ h H 1/2 (Γ ) + |μ| . (4.56) It remains to estimate the first term on the right-hand side of (4.55). With an inverse estimate and Lemma 4.8 we get We use the abbreviation e ν :=ν − ν h . The ellipticity of W on H 1/2 (Γ )/R and the definition of the Galerkin projection Π W imply With the commutator C η we get W ( ηe ν ), Π W ( ηe ν ) = ηW (e ν ) + C η e ν , Π W ( ηe ν ) .
For the term involving C η in (4.58), we get with Lemma 3.6 A duality argument implies e ν H −α N (Γ ) h 1/2+α N ν H 1/2 (Γ ) , for details we refer to the proof of Corollary 2.9. Inserting everything in (4.57) leads to Finally, this implies with (4.55) and (4.56) that which proves the lemma.
We want to use Lemma 4.10. Since [γ 0 K ζ h ] = ζ h ∈ S 1,1 (T h ) for any discrete function ζ h ∈ S 1,1 (T h ), we need to construct a discrete function satisfying the orthogonality (4.33). Using the Galerkin orthogonality with test functions with supp ψ h ⊂ Γ 2 and noting that η 3 ≡ 1 on supp ψ h , we obtain with the commutator C η 3 defined in (3.18), the abbreviation η 3 C η 3 e = 1 |Γ | η 3 C η 3 e, 1 , and the Galerkin projection Π W from (4.36) (4.59) Here and below, we understand the inverse W −1 as the inverse of the bijective operator W : H no additional terms in the orthogonality (4.59) appear. Thus, defining we get on a volume box B 2 ⊂ R d a discrete harmonic function Using the H 1 -stability of the Galerkin projection Π W , the mapping properties of W −1 and C η 3 as well as Lemma 3.6, the correction ξ h can be estimated by For the second term on the right-hand side of (4.60) we have η 1 ζ h H 1 (Γ ) η 1 ∇ζ h L 2 (Γ ) + ζ h L 2 (Γ ) . We apply Lemma 4.10 to u = K ζ h ∈ H N h (B 2 , μ) and obtain Since we only used the Galerkin orthogonality as a property of the error e, we may write ϕ − ϕ h = (ϕ − χ h ) + (χ h − ϕ h ) for arbitrary χ h ∈ S 1,1 (T h ) with supp χ h ⊂ Γ and we have proven the claimed inequality. where the second estimate is the standard global error estimate for the Galerkin BEM applied to the hyper-singular integral equation, see [22]. For the remaining term in Theorem 2.8, we use a duality argument. Let ψ solve W ψ = w − w ∈ H α N (Γ ), ψ, 1 = 0, where w = 1 |Γ | w, 1 . Then ψ ∈ H 1+α N (Γ ), and since e, 1 = 0, we get with the Scott-Zhang projection J h and Lemma 4.9 Therefore, the term of slowest convergence is of order O(h min{1/2+α+α N ,β} ), which proves the corollary.

Numerical examples
In this section we provide numerical examples to illustrate the theoretical results of Sect. 2 and indicate their sharpness. We only consider Symm's integral equation on quasi-uniform meshes. Provided the right-hand side and the geometry are sufficiently smooth, it is well-known that the lowest order boundary element method in two dimensions converges in the energy norm as O(N −3/2 ), where N denotes the degrees of freedom. In our examples we will consider problems, where the rate of convergence with uniform refinement is reduced due to singularities. In order to compute the error between the exact solution and the Galerkin approximation, we prescribe the solution u(r, θ) = r α cos(αθ ) of Laplace's equation in polar coordinates. Then, the normal derivative φ = ∂ n u of u is the solution of V φ = (K + 1/2)γ 0 u.
The regularity of φ is determined by the choice of α. In fact, we have φ ∈ H −1/2+α−ε (Γ ), ε > 0, and locally φ ∈ H 1 ( Γ ) for all subsets Γ ⊂ Γ that are a positive distance away from the singularity at the origin.
The lowest order Galerkin approximation to φ is computed using the MATLABlibrary HILBERT [2], where the errors in the L 2 -norm are computed using two point Gauss-quadrature. The error in the local H −1/2 -norm is computed as χ e 2 H −1/2 (Γ ) ∼ V (χ e), χ e , where χ is the characteristic function for a union of elements Γ 0 ⊂ Γ .

Example 1: L-shaped domain
On the L-shaped domain depicted in Fig. 1 (left), the dual problem permits solutions of regularity H 1/6−ε (Γ ) for arbitrary ε > 0; that is, we have α D = 1 6 − ε.    Figure 2 shows the global convergence in the energy norm (blue squares) as well as the local convergence on the fat, red part of the boundary (Γ 0 , union of elements) in the L 2 -norm (red stars) as well as the H −1/2 -norm (brown triangles). The black dotted lines mark the reference curves of order N −β for various β > 0.
In the left plot of Fig. 2 we chose α = 1 3 , which leads to α +α D = 1 2 −ε and, indeed, we observe convergence in the local L 2 -norm of almost order 1, which coincides with the theoretical rate obtained in Corollary 2.4. The error in the local H −1/2 -norm is smaller than the error in the L 2 -norm, but does converge with the same rate, i.e., an improvement of Theorem 2.3 in the energy norm is not possible. The right plot in Fig. 2 shows the same quantities for the choice α = 1 8 . Obviously, in this case the rates of convergence are lower, and the local L 2 -error does not converge with the best possible rate of one, but rather with the expected rate of N −19/24 = N −1/2−α−α D , as predicted by Corollary 2.4.

Example 2: Z-shaped domain
We consider the Z-shaped geometry depicted in Fig. 1 (right). Here, the dual problem permits solutions of regularity H α D (Γ ) with α D = 1 14 − ε. Again, we observe the