Functional a posteriori error estimates for boundary element methods

Functional error estimates are well-estabilished tools for a posteriori error estimation and related adaptive mesh-refinement for the finite element method (FEM). The present work proposes a first functional error estimate for the boundary element method (BEM). One key feature is that the derived error estimates are independent of the BEM discretization and provide guaranteed lower and upper bounds for the unknown error. In particular, our analysis covers Galerkin BEM as well as collocation, which makes our approach of particular interest for people working in engineering. Numerical experiments for the Laplace problem confirm the theoretical results.


Introduction
Let Ω ⊂ R d , d ≥ 2, be a bounded Lipschitz domain with polygonal boundary Γ := ∂Ω.We consider the Poisson problem with inhomogeneous Dirichlet boundary data g, i.e., ∆u = 0 in Ω, u = g on Γ. (1) Throughout the paper, we assume that d ∈ {2, 3}.However, all results can easily be extended to higher dimensions.In case of d = 2, we suppose that diam(Ω) < 1, which can always be achieved by scaling.For the numerical solution of (1), we employ the boundary element method (BEM); see, e.g., [Ste08a,SS11,GS18].For the ease of presentation, let us consider an indirect ansatz based on the single-layer potential with unknown integral density φ, where G(z) = − 1 2π log |z| for d = 2 resp.G(z) = 1 4π |z| −1 for d = 3 denotes the fundamental solution of the Laplacian.Taking the trace on Γ, the potential ansatz leads to the weakly-singular integral equation where the integral representation of g = V φ coincides with that of u = V φ (at least for bounded densities) but is now evaluated on Γ (instead of inside Ω).Given a triangulation F Γ h of the boundary Γ, the latter equation is solved by the lowest-order BEM and provides some piecewise constant approximation φ h , i.e., where the precise discretization (e.g., Galerkin BEM, collocation, etc.) will not be exploited by our analysis.However, as a BEM inherent characteristics, we obtain an approximation of the potential u ≈ u h := V φ h , which satisfies the Laplace problem ∆u h = 0 in Ω.The latter is the key argument for the error identity where • , • Γ denotes the extended L 2 (Γ) scalar product (see Theorem 3 below).The identities (5) generate a posteriori error estimates of the functional type that are independent of the discretization and provide fully guaranteed lower and upper bounds for the unknown error without any constants at all.In general, these functional type a posteriori estimates involve only constants in basic functional inequalities associated with the concrete problem (e.g., Poincaré-Friedrichs type or trace inequalities) and are applicable for any approximation from the admissible energy class (see [Rep00,AP19] or the monograph [Rep08] and the references cited therein).In particular, the equations ( 5) have also been used in [Rep08] for the analysis of errors arising in the Trefftz method.From (5), constant-free (i.e., with known constant 1) lower and upper bounds for the unknown potential error ∇(u − u h ) L 2 (Ω) can be obtained by choosing arbitrary instances of τ and w.In the present work, we compute these bounds by solving problems in a suitable boundary layer S ⊂ Ω along Γ by use of the finite element method (FEM).Moreover, these bounds are then employed to drive an adaptive mesh-refinement for the triangulation F Γ h of Γ and, as novelty, also quantify the accuracy ||∇(u − u h )|| L 2 (Ω) of the BEM induced potential u h = V φ h in each step of the adaptive algorithm.In particular, the latter quantification is essentially constant-free (up to data oscillations terms arising for the FEM majorant) and can thus also be used as reasonable stopping criterion for adaptive BEM computations.Especially for practical applications, this is an important step forward, since there exist neither a posteriori error estimates with constant 1 nor estimates for physically relevant errors.While available results focus on the density φ h (see, e.g., [CS95, Car97, MSW98, CMS01, EH06, FLP08] for some prominent results or the surveys [CF01, FFH + 15] and the references therein), estimating rather the energy error of u h circumvents, in particular, BEM-natural challenges like the localization of non-integer Sobolev norms.
Outline.The remainder of this work is organized as follows: In Section 2, we collect the necessary notations as well as the fundamental properties of (Galerkin) BEM.In Section 3, we formulate our approach for functional a posteriori error estimation.Theorem 3 states the error identity (5).Theorem 4 provides a computable upper bound on (5) by means of an H 1 -conforming FEM approach as well as a computable lower bound on (5) by means of an H(div)-conforming mixed FEM approach.Section 4 shows how these findings can be used to steer an adaptive mesh-refinement.Algorithm 9 formulates such a strategy with reliable error control on ||∇(u − u h )|| L 2 (Ω) .In Section 5, we employ the proposed adaptive algorithm to underpin our theoretical findings by some numerical experiments in 2D.Section 6 concludes the work with possible extensions of our analysis like higher-order BEM, alternative BEM discretizations like collocation, direct BEM formulations, and error control for exterior domain problems (where Ω is unbounded).The final Section 7 summarizes the contributions of the present work and addresses possible topics for future research.

Preliminaries and notation
2.1.Domains and function spaces.Throughout this paper, let Ω ⊂ R d , d ∈ {2, 3}, be a bounded Lipschitz domain (i.e., locally below the graph of some Lipschitz function) with boundary Γ = ∂Ω and normal vector field n.For all numerical results involving discretisations, we assume that Γ is a polygon.We denote by • , • L 2 (Ξ) and || • || L 2 (Ξ) the standard inner product and norm in L 2 (Ξ), respectively, where, e.g., Ξ ∈ {Ω, Γ}.Based on L 2 (Ω), we define the Hilbert spaces The corresponding inner products and (induced) norms are together with its range H 1/2 (Γ) and equipped with the natural quotient norm, our analysis also employs which is a closed subspace of H 1 (Ω).With the help of an extension operator (•) : , which comes together with the trace operator as a right inverse, we can construct a (unique) harmonic extension Finally, we need the dual space , where the is a Gelfand triple and refer to [BBF13] for the fact that H −1/2 (Γ) can also be characterised as the range of normal traces n Definition 1.A subset S ⊂ Ω is called a boundary layer, if it is a Lipschitz domain with Γ ⊂ ∂S, which admits a conforming triangulation T S h into simplices.We then define In particular, we define the corresponding induced triangulation of Γ by For q ∈ N 0 and P q being the space of polynomials of degree q, we define Moreover, for p ∈ N, we employ the standard H 1 -conforming FEM spaces Let F S h denote the set of all interior faces, i.e., all F ∈ F S h admit unique T + , T − ∈ T S h with F = T + ∩ T − .For q ∈ N 0 , we define the H(div)-conforming Raviart-Thomas space where n F is a normal vector for the face Remark 2. In the proofs of Section 3 below, we exploit that for arbitrary v h ∈ S p Γ c (T S h ) and σ h ∈ RT q Γ c (T S h ) the definitions provide conforming extensions vh ∈ H 1 (Ω) and σh ∈ H(div, Ω).In particular, we will implicitly identify v h (resp.σ h ) with its zero-extension vh (resp.σh ).

Functional a posteriori BEM error estimation
In this section, we prove the error identity (5) and provide efficiently computable upper and lower bounds for the potential error ||∇(u − u h )|| L 2 (Ω) , where u ∈ H 1 (Ω) solves (9) and u h := V φ h is defined in (2).
3.1.Functional error identity.The fact that the error u − u h satisfies (9a) exactly is a powerful tool.However, the consideration of the potential u h from a BEM comes with a drawback: it is not a discrete function and lacks further a priori knowledge like the Galerkin orthogonality, which is obviously never available for any approximation u h := V φ h ≈ u ∈ H 1 (Ω).Functional a posteriori error estimates are eminently suitable for the BEM, since they do not require any such a priori assumption.On top of that, for homogeneous problems, they provide constant-free error identities by means of so-called primal and dual problems.For the Laplacian, the key argument is the Dirichlet principle: Harmonic functions are minimisers of the Dirichlet energy ||∇w|| 2 L 2 (Ω) .Note that the boundary residual g − u h | Γ ∈ H 1/2 (Γ) is essential for both the majorant M and the minorant M, and comprises all relevant information about the error.
Proof.The proof is split into two parts.
• Upper bound: Let w ∈ H 1 (Ω) with w| Γ = u| Γ = g.Since we have ∆(u − v) = 0 and u − w ∈ H 1 0 (Ω), integration by parts shows that This yields . The substitution w := w − v proves that The unique infimum is attained at w = u − v.
• Lower bound: In any Hilbert space H, it holds that where the maximum is unique and attained for b = a.Since we have In particular, the maximum is attained for τ = ∇(u − v).This concludes the proof.
To ease the readability, the remainder of this chapter focusses on our numerical setup.For the functional analytic framework in a Sobolev space setting, which might be of independent interest, we refer to Appendix A below.

Computable error bounds.
We aim at error bounds obtained by solving FEM problems on a boundary layer S ⊂ Ω.For the maximization problem in (12), the constraint div τ = 0 can be realized by a mixed formulation (see also Lemma 15 in Appendix A).However, the boundary condition w| Γ = g − v| Γ cannot be satisfied exactly by any piecewise polynomial solution w h corresponding to (12).Therefore, the upper bound involves an additional oscillation term given by a discretisation operator J h .

Theorem 4 (Computable bounds via boundary layer
For q ∈ N 0 , let the pair Then, it holds that Proof.It is well-known that (13) admits a unique solution w h ∈ S p (T S h ), being the natural FEM discretization of an homogeneous Dirichlet-Laplace problem with inhomogeneous Dirichlet conditions; see, e.g., [BCD04, SV06, AFK + 13].To prove the upper bound (16), let f h ∈ H 1 (Ω) be the (unique) harmonic extension of , and Theorem 3 leads to Since the zero-extension of w h belongs to H 1 (Ω) according to Remark 2, this proves that and hence verifies the computable upper bound (16).
In order to circumvent the implementation of the constraint div τ = 0, it is also an option to reformulate the maximization problem in (12) by means of potentials.While the 3D case involves vector potentials, for 2D such an approach is particularly attractive due to the possible use of scalar potentials.In the following, we thus concentrate on d = 2 (and refer to Appendix A for d = 3).
To this end, we recall the definitions of the 2D curl operators Note that div curl ϕ = 0.For ϕ ∈ H 1 (Ω), we thus have curl ϕ ∈ H(div, Ω) so that the Neumann trace

Corollary 5 (Computable lower bound via boundary layer -H
Then, it holds that Proof.It is well-known that (17) admits a unique solution w h ∈ S p Γ c (T S h ) being the natural FEM discretization of a mixed Dirichlet-Neumann-Laplace problem; see, e.g., [BCD04].According to Remark 2, the zero-extension of w h belongs to H 1 (Ω) and hence . The claim thus follows from Theorem 3.

Adaptive algorithm
4.1.Triangulations and mesh-refinement.In our numerical experiments, we start from a conforming simplicial triangulation T h such that Γ ⊂ T ∈T h T ⊆ Ω.We obtain the boundary layer S ⊂ Ω as the second-order patch of Γ with respect to T h , i.e., Moreover, recall the BEM mesh ).These definitions are illustrated in Figure 1.
For (local) mesh-refinement, we employ newest vertex bisection; see, e.g., [Ste08b,KPP13].The adaptive strategy will only mark elements of T S h , but refinement will be done with respect to the full triangulation T h .In particular, we stress that the secondorder patch S will generically change, if the triangulation T h is refined; see, e.g., Figure 2. In this way, we guarantee that the number of degrees of freedom with respect to T S h will increase proportionally to those with respect to F Γ h ; see also Tables 2-5 below.

Data oscillations. The upper bound (16) in Theorem 4 involves the data approximation term (1
, where u h := V φ h .Besides the fact that we still have to specify the operator , we note that the nonlocal nature of the H 1/2 (Γ)-norm makes this term hardly computable.
In the following, we choose see [CT87].We note that these conditions are automatically satisfied for known for low-order FEM (on the 2D manifold Γ); see [KPP13] for p = 1 and [GHS16] for p ∈ {1, . . ., 12}.We recall the following result from [AFF + 15]: where h ∈ L ∞ (Ω) is the local mesh-width function defined by h| F := diam(F ) for all F ∈ F Γ h .The constant C osc > 0 depends only on C stab and the shape regularity of T h .Provided that the given Dirichlet boundary data satisfy g ∈ H 1 (Γ), the foregoing lemma allows to dominate the data approximation term by where osc h is, in fact, computable, while the constant C osc is generic and hardly accessible.With v = u h , the upper bound (16) becomes where w h ∈ S p (T S h ) solves (13).For the use in the adaptive algorithm, we note that Remark 7. In our numerical experiments, we will consider p = 1 as well as p = 2 to compute the uppermost bound in (21).Since the lower bound (15) is independent of the data approximation, we did only implement the lowest-order case q = 0.
Remark Algorithm 9. Let p ∈ N and let 0 < θ ≤ 1 be a fixed marking parameter.Let T h be a conforming initial triangulation of Ω.Let ε > 0 be the tolerance for the energy error ∇(u − u h ) L 2 (Ω) with u h = V φ h .Then, perform the following steps (i) -(ix): h ) of (13) for the majorant (16).(vi) Compute the error indicators Remark 10.We note that a reliable adaptive algorithm mandatorily relies on a computable upper bound.Optionally, one can furthermore include the computation of a lower bound to provide a guaranteed confidence interval for the error: Computing the FEM solution (τ h , φ h ) ∈ RT 0 Γ c (T S h ) × P 0 (T S h ) of (14), we can also assemble the discrete minorant from Theorem 4. This is also done in the numerical experiments in Section 5.While
Even though u as well as its Dirichlet data g = u| Γ are smooth, we note that the sought integral density φ ∈ H −1/2 (Γ) of the indirect formulation (3) has no physical meaning and usually lacks smoothness (by inheriting the generic singularities from the interior as well as the exterior domain problem).Consequently, one may expect that uniform mesh-refinement (on the boundary) will not reveal the optimal convergence behavior  16), and the width of the confidence interval given by the quotient of majorant and minorant.16), and the width of the confidence interval given by the quotient of majorant and minorant.
), where N = #F Γ h is the number of elements of a uniform mesh F Γ h of Γ and 3/2 is the best possible convergence rate for a piecewise constant approximation φ ≈ φ h ∈ P 0 (F Γ h ).The initial meshes and some adaptively generated meshes are visualized in Figure 2.
Figure 3 shows the resulting potential error and the computed minorant (15) and majorant (16), as well as the corresponding data oscillations (20) for p = 1 resp.p = 2. Here, the potential error is computed by numerical quadrature.More precisely, we employ the P 2 -nodal interpolant I h : C(Ω) → S 2 (T unif h ) on a (three times) uniform refinement T unif h of the finest adaptive mesh T h .We stress that the plot neglects the non-accessible constant C osc from (20).The results for p = 1 and p = 2 are similar.For uniform mesh-refinement, we obtain the expected reduced order of convergence.For adaptive mesh-refinement, we regain the optimal order of convergence.Moreover, for adaptive mesh-refinement, we see that the majorant is, in fact, a sharp estimate for the (in general unknown) potential error.The computed minorant is less accurate.However, we stress that the minorant is always computed on the same boundary layer as the majorant (which is obtained by adaptivity driven by the majorant).The empirical values for uniform (resp.adaptive) mesh-refinement are also provided in Table 1 (resp.Table 2).The accuracy of the lower bound could be improved by appropriately adapting a second boundary layer (not displayed).Figure 4 compares the numerical results for different choices of the adaptivity parameter θ ∈ {0.2, 0.4, 0.6, 0.8}.We observe that any choice of θ regains, in fact, the optimal convergence rate.
Example 5.2 (Smooth potential in L-shaped domain).We consider problem (1) with prescribed exact solution on the L-shaped domain Ω with diameter diam(Ω) = 1/2.We start Algorithm 9 with an initial triangulation T h of Ω into #T 0 = 384 right triangles.As in Section 5.1, the potential u is smooth, but the sought density φ of the indirect BEM formulation lacks regularity.The initial meshes as well as some adaptively generated meshes are visualized in Figure 5. Figure 6 visualizes some numerical results for uniform and adaptive mesh-refinement, where we proceed as in Section 5.1.Since p = 1 and p = 2 lead to similar results (not displayed), we only report the results for p = 1.
As expected from theory, the shape of Ω does not impact the functional error estimates: Overall, the results obtained correspond to those from Section 5.1, where uniform meshrefinement leads to a suboptimal convergence behavior, which is cured by means of the proposed adaptive strategy.16), the data oscillations osc h from (20), and the minorant M(τ h ) 1/2 from (25) for uniform (dashed) and adaptive mesh-refinement (solid) with θ = 0.4.
Right: We compare the potential error (solid) and the majorant (dashed) for adaptive mesh-refinement for various choices of θ.
= 0 Example 5.3 (Non-smooth potential in L-shaped domain).We consider problem (1) with prescribed exact solution given in standard polar coordinates x = x(r, ϕ) on the L-shaped domain Ω with diameter diam(Ω) = 1/2.We start Algorithm 9 with an initial triangulation T 0 of Ω into #T 0 = 384 right triangles.
Right: We compare the potential error (solid) and the majorant (dashed) for adaptive mesh-refinement for various choices of θ.
convergence results are visualized in Figure 8.Moreover, Table 3 provides some empirical values for adaptive mesh-refinement.Our observations are the same as in Section 5.1 and Section 5.2 and underline that the functional error bounds do not rely on any a priori smoothness of the unknown potential u: While uniform mesh-refinement leads to a suboptimal convergence behavior, the proposed adaptive strategy regains the optimal convergence rate.

Extension of the analysis
So far, we have considered functional a posteriori error estimation for an indirect BEM formulation (10) discretized by Galerkin BEM (11).The following sections address some obvious extensions of our analysis.While the subsequent numerical experiments (as well as those from Section 5) focus on d = 2, we again stress that the theoretical results also apply to arbitrary dimensions, in particular to d = 3.However, 3D experiments are beyond the scope of this work and left to future research.

Collocation BEM.
It is worth noting that all results of Section 3 hold, in particular, for any v = V φ h with arbitrary φ h ∈ H −1/2 (Γ).Consequently, the computable bounds of Theorem 4 (resp.Corollary 5) hold for any approximation φ h ≈ φ.In particular, Algorithm 9 can also be applied to (e.g., lowest-order) collocation BEM, where  16), and the width of the confidence interval given by the quotient of majorant and minorant.
where x F ∈ F is an appropriate collocation node (e.g., the center of mass).We stress that well-posedness of collocation BEM is non-obvious (see, e.g., [CPS92,CPS93,MP96]).However, this does not affect our developed functional a posteriori error bounds.
6.2.Other BEM ansatz spaces.With the same argument as for collocation BEM, one can replace the discrete BEM ansatz space P 0 (F Γ h ) φ h by an arbitrary discrete space P h ⊆ H −1/2 (Γ) (e.g., higher-order piecewise polynomials, splines, isogeometric NURBS, etc.).For r ∈ N 0 and P h = P r (F Γ h ), we expect that the choices p = r + 1 and q = r will lead to accurate computable upper and lower bounds in Theorem 4. The numerical validation of this expectation is, however, beyond the scope of the present work.6.3.Direct BEM approach.The indirect BEM approach makes ansatz (10) for the unknown solution of (9).Unlike this, the direct BEM approach is based on the Green's third identity: Any solution of (9) can be written as the sum of a single-layer and a double-layer potential, i.e., where g = u| Γ ∈ H 1/2 (Γ) is the trace of u (i.e., the Dirichlet data) and φ = n • ∇u| Γ ∈ H −1/2 (Γ) is the normal derivative (i.e., the Neumann data).Taking the trace of this identity and respecting the jump properties of the double-layer potential (see,  16) based on P 1 -FEM, the data oscillations osc h from (34), and the minorant M(τ h ) 1/2 from (25).Adaptive mesh-refinement with θ = 0.4 in Example 6.4.We focus on the degrees of freedom, the potential error ||∇(u − u h )|| L 2 (Ω) , the accuracy of the P 1 -FEM majorant ||∇w h || L 2 (Ω) from ( 16), and the width of the confidence interval given by the quotient of majorant and minorant.e.g., [McL00, Ste08a, HW08, SS11, GS18]), one sees that Table 5. Adaptive mesh-refinement with θ = 0.4 in Example 6.5.We focus on the degrees of freedom, the potential error ||∇(u − u h )|| L 2 (Ω) , the accuracy of the P 1 -FEM majorant ||∇w h || L 2 (Ω) from ( 16), and the width of the confidence interval given by the quotient of majorant and minorant.
where K formally coincides with K, but is evaluated for x ∈ Γ instead.Elementary calculations then lead to the variational formulation We stress that the factor 1/2 is only valid almost everywhere on Γ and hence correct for the variational formulation and Galerkin BEM, while collocation BEM would require a modification at corners (and additionally along edges in 3D); see [McL00,Ste08a,HW08]. Usual implementations approximate g ≈ g h ∈ S p (F Γ h ) so that the integral operators in (31) are only evaluated for discrete functions.Overall, the lowest-order Galerkin BEM formulation then reads As above, the Lax-Milgram lemma proves that (31) (resp.(32)) admit unique solutions φ ∈ H −1/2 (Γ) (resp.φ h ∈ P 0 (F Γ h )).Moreover, the computed density φ h is now indeed an approximation of the Neumann data n one obtains an approximation u h of the solution u = V φ − Kg of (9) (resp.(30)).We stress that u h | Γ = V φ h +(1/2−K)g h so that the data oscillation term in the upper bound of Theorem 4 reads where see Section 4.2.In our implementation, we also employed g h = J h g ∈ S 1 (F Γ h ).Example 6.4 (Direct BEM for smooth potential in square domain).We consider the setting (26) from Section 5.1.Applying the direct BEM approach (31), we know that φ h ≈ φ = n • ∇u| Γ , where φ (as well as the potential u) is smooth.In this particular situation, we know that uniform mesh-refinement would already lead to the optimal convergence behavior (not displayed).The same is empirically observed for the proposed adaptive strategy, where we even observe that the majorant ||∇w h || L 2 (Ω) from (16) as well as the minorant M(τ h ) 1/2 from (25) provide sharp error bounds for the potential error ||∇(u − u h )|| L 2 (Ω) ; see Figure 9 (left) as well as Table 4.
Example 6.5 (Direct BEM for non-smooth potential in L-shaped domain).We consider the setting (28) from Section 5.3.Applying the direct BEM approach (31), we know that φ h ≈ φ = n • ∇u| Γ , where φ (as well as the potential u) is only nonsmooth with a singularity at (0, 0).Also for this case, the proposed adaptive strategy regains the optimal convergence rate; see Figure 9 (right) as well as Table 5.Even though the quotient ||∇w h || L 2 (Ω) /M(τ h ) 1/2 of the computable upper and lower bound is larger than for the smooth problem of Section 6.4, we observe that the lower bound is, in fact, much more accurate for the direct BEM than for the indirect BEM computations from Section 5.
where Ω is the L-shaped domain from Section 5.2.We consider (35) with constant Dirichlet data where K is the double-layer integral operator.Consequently, the corresponding indirect BEM formulation (10) turns out to be a direct BEM formulation for the exterior domain problem [McL00,SS11], where all data oscillation terms vanish.Thus, one can expect that the sought density φ ∈ H −1/2 (Γ) has singularities at the convex corners of Ω (but not at the reentrant corner).We employ Algorithm 9 (with Galerkin BEM).The initial mesh T h with #T h = 416 right triangles is a triangulation of (−1/4, 3/4) 2 \Ω ⊂ Ω c ; see Figure 10.Some numerical results are shown in Figure 11.Since the exact potential u is unknown, we cannot compute the potential error ||∇(u − u h )|| L 2 (Ω) .However, adaptive mesh-refinement leads to the optimal convergence behavior of majorant and minorant (and hence also of ||∇(u − u h )|| L 2 (Ω) ).

Conclusion
We have presented, for the first time, functional error estimates for BEM.Not only that the presented estimates are independent of the specific discretization method (i.e., Galerkin or collocation), they also provide guaranteed upper and lower bounds for the unknown energy error.This is in contrast to existing techniques, which usually contain generic constants.The error bounds are obtained by solving auxiliary variational problems by a finite element method on a strip.In the paper, we consider the Dirichlet problem of the Laplace equation, but the approach is expected to generalize to other boundary value problems.In the considered case, the upper error bound is based on the Dirichlet principle, while the lower error bound is based either on a variational problem in terms of a potential (scalar stream function in 2D and vector potential in 3D) or a mixed problem (in 2D and 3D).The upper bound is localized and drives an adaptive refinement of the boundary mesh.Since the finite element strip contains always two layers of elements, it geometrically shrinks towards the boundary during refinement.This way, the ratio between the degrees of freedom (DoF) for obtaining the error estimates and the BEM DoF remains bounded.We have examined various 2D test problems on square and L-shaped domains, with and without singular potential, including exterior problems.The proposed adaptive algorithm exhibited excellent performance.In all cases, the optimal convergence rates could be achieved.
Suggested further research work concerns an implementation of the algorithm in 3D and the extension to electromagnetic scattering problems.

Appendix A. Some remarks on the analysis
We recall the notations introduced in Section 2.1, in particular, the harmonic extension operator (•) : H 1/2 (Γ) → H 1 (Ω) from (6).Moreover, we add the definitions of see, e.g., [BPS16], for the density result, and .
The latter space generalizes the (partial) homogeneous boundary condition n • σ| Γ c = 0 to functions σ ∈ H(div, S).Note that the natural trace n • σ| ∂S can be 'restricted' to, e.g., Γ in the following sense.
Remark 11.Functions in H1 Γ c (S) vanish at Γ c .In particular, any function ϕ ∈ H 1 Γ c (S) can be extended by zero to a function ϕ ∈ H 1 (Ω).Analogously, vector fields in H Γ c (div, S) have vanishing normal component at Γ c in a weak sense.In particular, any vector field σ ∈ H Γ c (div, S) can be extended by zero to a vector field σ ∈ H(div, Ω).Hence, the normal trace n • σ| ∂S ∈ H −1/2 (∂S) of σ may be identified with a well defined element n • σ| Γ ∈ H −1/2 (Γ) vanishing on Γ c in a weak sense.
Let us discuss the minimiser w = u − v of the upper bound and the maximiser τ = ∇w of the lower bound from Theorem 3 in more detail.Note that Moreover, by replacing w with w + εϕ and by replacing τ with τ + εσ in (12), where ϕ ∈ H 1 0 (Ω) and σ ∈ H(div, Ω) with div σ = 0 as well as ε ∈ R, we obtain the variational formulations Let ψ ∈ H 1/2 (Γ) and let ψ ∈ H 1 (Ω) be its harmonic extension.Testing the second variational formulation (38b) with σ = ∇ ψ shows Thus, additionally to the scalar and tangential boundary conditions for w and τ in (37), respectively, we have also found a normal boundary condition for τ = ∇w, namely This shows that there are different options for computing w and τ .
Remark 12.Note that is the well known Dirichlet-to-Neumann operator for the homogeneous Laplacian.Moreover, the normal trace of τ in (39) does not depend on the harmonic extension as Recalling Theorem 3 and Theorem 4, we note the following.
Remark 13 (Minimiser of the upper bound).The unique minimiser w of the upper bound is the unique harmonic extension of g − v| Γ to Ω, i.e., w ∈ H 1 (Ω) is the unique solution of the Dirichlet-Laplace problem Note that, at least analytically, both formulations can also be used to find the unique maximiser τ = ∇w of the upper bound.For numerical purposes the Dirichlet-Laplace problem is the better and easier choice to compute w.
Next we want to find equations and variational formulations for τ not involving w.As τ ∈ H(div, Ω) with div τ = 0, by (38b) the vector field τ solves for all ω ∈ L 2 (Ω) the mixed problem for all (σ, ψ) ∈ H(div, Ω) × L 2 (Ω) with div σ = 0. On the other hand, especially for numerical reasons, we want to skip the solenoidal conditions, which leads to the following mixed variational saddle point formulation (cf.Theorem 4).
Remark 16.The mixed formulation in Lemma 15 is the mixed formulation of the Dirichlet-Laplace problem from Remark 13.For numerical purposes the latter mixed formulation is only a good choice to compute τ , since the numerical approximations will only satisfy (τ , ω) ∈ H(div, Ω) × L 2 (Ω) but in general not ω ∈ H 1 (Ω).

Figure 2 .Figure 3 .
Figure 2. Adaptively generated meshes in Example 5.1 for p = 1 and θ = 0.6.We indicate the boundary layer S (blue), the boundary Γ (red), and the interior boundary Γ c = ∂S \ Γ (green).The triangles T ∈ T S h ⊂ T h are indicated in blue.The triangles T ∈ T h \T S h are indicated in gray.

Figure 5 .
Figure 5. Adaptively generated meshes in Example 5.2 for p = 1 and θ = 0.6; see Figure 2 for the color code.

Figure 7 .
Figure 7. Adaptively generated meshes in Example 5.3 for p = 1 and θ = 0.6; see Figure 2 for the color code.
necessarily satisfied, this is empirically observed throughout all experiments below.Moreover, if one aims for sharper error bounds for a fixed approximate solution u h , one can solve (13) and (14) (resp.(17)) by adaptive FEM on separate (generically) different boundary layers [Seb19].