Adaptive BEM for elliptic PDE systems, part II: Isogeometric analysis with hierarchical B-splines for weakly-singular integral equations

We formulate and analyze an adaptive algorithm for isogeometric analysis with hierarchical B-splines for weakly-singular boundary integral equations. We prove that the employed weighted-residual error estimator is reliable and converges at optimal algebraic rate. Numerical experiments with isogeometric boundary elements for the 3D Poisson problem confirm the theoretical results, which also cover general elliptic systems like linear elasticity.

1. Introduction 1.1. State of the art. The central idea of isogeometric analysis (IGA) is to use the same ansatz functions for the discretization of the partial differential equation (PDE) at hand, as are used for the representation of the problem geometry. Usually, the problem geometry Ω is represented in computer aided design (CAD) by means of non-uniform rational B-splines (NURBS), T-splines, or hierarchical splines. This concept, originally invented in [HCB05] for finite element methods (IGAFEM) has proved very fruitful in applications; see also the monograph [CHB09].
Usually, CAD programs only provide a parametrization of the boundary ∂Ω instead of the domain Ω itself. In particular, for isogeometric FEM, the parametrization needs to be extended to the whole domain Ω, which is non-trivial and still an active research topic. The boundary element method (BEM) circumvents this difficulty by working only on the CAD provided boundary ∂Ω. However, compared to the IGAFEM literature, only little is found for isogeometric BEM (IGABEM). Based on a collocation approach, the latter was first considered in [PGK + 09] for 2D and in [SSE + 13] for 3D. We refer to the monograph [BMD20] for an introduction on IGABEM and to [DHK + 18, DHK + 20] for an efficient IGABEM implementation based on the fast multipole method. However, to the best of our knowledge, a posteriori error estimation and adaptive mesh-refinement for IGABEM, have only been considered for simple 2D model problems in the own works [FGP15, FGHP16, FGHP17, FGPS19, GPS20], which employ B-splines on the one-dimensional boundary, and the recent work [FGK + 19], which employs hierarchical B-splines. In particular, [FGP15] also appears to be the first work that considers Galerkin IGABEM.
For standard BEM with (dis)continuous piecewise polynomials, a posteriori error estimation and adaptive mesh-refinement are well understood. We refer to [FFH + 15] for an overview on available a posteriori error estimation strategies. Moreover, optimal convergence of mesh-refining adaptive algorithms has been proved for polyhedral boundaries [FFK + 14, FFK + 15, FKMP13] as well as smooth boundaries [Gan13]. The work [AFF + 17] allows to transfer these results to piecewise smooth boundaries; see also the discussion in the review article [CFPP14].
In the frame of adaptive IGAFEM, a rigorous error and convergence analysis is first found in [BG16b], which proves linear convergence for some adaptive IGAFEM with hierarchical splines for the Poisson equation. Optimal algebraic rates have been proved independently from each other in [BG17,GHP17] for IGAFEM with (truncated) hierarchical B-splines. In particular, [GHP17] especially identified certain abstract properties on the mesh-refinement and the used ansatz spaces that automatically guarantee optimal convergence. In [GP20b], we have recently verified these properties for IGAFEM with T-splines as well.
Unlike adaptive IGAFEM, the rigorous error and convergence analysis of adaptive IGA-BEM is completely open beyond the 2D results mentioned above and thus the focus of the present work. In Part 1 of this work [GP20a], we followed the idea of [GHP17] and identified an abstract framework, which guarantees optimal convergence of adaptive BEM for weaklysingular integral equations. In the present work (Part 2), we now show that IGABEM with hierarchical splines is indeed covered by this abstract framework.
1.2. Outline. The remainder of this manuscript is roughly organized as follows: Section 2 recalls conforming BEM for weakly-singular integral equations; see (2.17) for the precise model problem. Section 3 defines hierarchical meshes and hierarchical splines on the boundary Γ and introduces some local mesh-refinement rule (Algorithm 3.7) which preserves admissibility in the sense of a certain mesh-grading property. The main result of Section 3 is Theorem 3.10, which states reliability (3.48) of the weighted-residual error estimator (3.43) as well as linear convergence (3.49) with optimal algebraic rates (3.50) of a standard adaptive algorithm (Algorithm 3.9) applied to the model problem at hand. Remark 3.13 extends the result to rational hierarchical splines. The proof of Theorem 3.10 (and Remark 3.13) is given in Section 5, which recalls the abstract framework from [GP20a] and subsequently verifies the corresponding properties in the considered isogeometric setting. Two numerical experiments in Section 4 underpin the theoretical results, but also demonstrate the limitations of hierarchical splines in the frame of adaptive BEM when the solution φ exhibits edge singularities.

Preliminaries
In this section, we fix some general notation, recall Sobolev spaces on the boundary, and precisely state the considered problem. Throughout the work, let Ω ⊂ R d be a bounded Lipschitz domain as in [McL00,Definition 3.28] with boundary Γ := ∂Ω.
We write A B to abbreviate A ≤ CB with some generic constant C > 0, which is clear from the context. Moreover A ≃ B abbreviates A B A. Throughout, mesh-related quantities have the same index, e.g., X • is the ansatz space corresponding to the mesh T • . The analogous notation is used for meshes T • , T ⋆ , T ℓ , etc.
. If σ > 0, and ω ⊆ Γ is an arbitrary measurable set, we define v H σ (ω) and |v| H σ (ω) similarly. With the definition (2.8) As before we abbreviate (2.10) The spaces H σ (Γ) can be also defined as trace spaces or via interpolation, where the resulting norms are always equivalent with constants which depend only on the dimension d and the boundary Γ. More details and proofs are found, e.g., in the monographs [McL00, SS11, Ste08].
2.3. Model problem. We consider a general second-order linear system of PDEs on the d-dimensional bounded Lipschitz domain Ω with partial differential operator Moreover, we assume that P is coercive on H 1 0 (Ω) D , i.e., the sesquilinear form is elliptic up to some compact perturbation. This is equivalent to strong ellipticity of the matrices A ii ′ in the sense of [McL00, page 119].
Let G : R d \ {0} → C D×D be a corresponding (matrix-valued) fundamental solution in the sense of [McL00,page 198], i.e., a distributional solution of PG = δ, where δ denotes the Dirac delta function. For ψ ∈ L ∞ (Γ) D , we define the single-layer operator as (2.13) According to [McL00,page 209 and 219-220] and [HMT09, Corollary 3.38], this operator can be extended for arbitrary σ ∈ (−1/2, 1/2] to a bounded linear operator (2.14) For σ = 0, [McL00, Theorem 7.6] states that V is always elliptic up to some compact perturbation. We assume that it is elliptic even without perturbation, i.e., For instance, this is satisfied for the Laplace problem or for the Lamé problem, where the case d = 2 requires an additional scaling of the geometry Ω; see, e.g., [Ste08, Chapter 6]. We stress that V is not necessarily self-adjoint; see [McL00]. The sesquilinear form (V · ; ·) is continuous due to (2.14), i.e., it holds with C cont : Given a right-hand side f ∈ H 1 (Γ) D , we consider the boundary integral equation Such equations arise from the solution of Dirichlet problems of the form Pu = 0 in Ω with u = g on Γ for some g ∈ H 1/2 (Γ) D ; see, e.g., [McL00,] for more details. The Lax-Milgram lemma provides existence and uniqueness of the solution φ ∈ H −1/2 (Γ) D of the equivalent variational formulation of (2.17) (Vφ ; ψ) = (f ; ψ) for all ψ ∈ H −1/2 (Γ) D . (2.18) In particular, we see that V : H −1/2 (Γ) D → H 1/2 (Γ) D is an isomorphism. In the Galerkin boundary element method, the test space H −1/2 (Γ) D is replaced by some discrete subspace Again, the Lax-Milgram lemma guarantees existence and uniqueness of the solution Φ • ∈ X • of the discrete variational formulation and Φ • can in fact be computed by solving a linear system of equations. We note the Galerkin orthogonality along with the resulting Céa-type quasi-optimality Note that (2.14) implies that VΨ • ∈ H 1 (Γ) D for arbitrary Ψ • ∈ X • . The additional regularity f ∈ H 1 (Γ) D instead of f ∈ H 1/2 (Γ) D is only needed to define the residual error estimator (3.43) below. For a more detailed introduction to boundary integral equations, the reader is referred to the monographs [McL00,SS11,Ste08].

Boundary element method with hierarchical splines
In this section, we consider Ω ⊂ R d with d ≥ 2. We introduce hierarchical splines on the boundary Γ and propose a local mesh-refinement strategy. To this end, we assume the existence of a mesh Γ m : m = 1, . . . , M of Γ, i.e., Γ = M m=1 Γ m with interfaces Γ m ∩ Γ m ′ having (d − 1)-dimensional measure zero for m = m ′ , such that each surface Γ m can be parametrized over Γ m := [0, 1] d−1 . We formulate an adaptive algorithm (Algorithm 3.9) for conforming BEM discretizations of our model problem (2.17), where adaptivity is driven by the residual a posteriori error estimator (see (3.43) below). The main result of this section is Theorem 3.10, which states that the employed estimator is reliable and converges linearly at optimal rate. The proof of Theorem 3.10 is given in Section 5.
3.1. Parametrization of the boundary. We assume that for all m ∈ {1, . . . , M}, the surface Γ m ⊆ Γ can be parametrized via a bi-Lipschitz mapping In particular, γ m is almost everywhere differentiable, and there exists a constant C γ ≥ 1 such that and the Gram determinant satisfies that see, e.g., [Gan17, Lemma 5.2.1]. We define the set of nodes We suppose that there are no (initial) hanging nodes, i.e., the intersection Γ m ∩Γ m ′ with m = m ′ is either empty or a common (transformed) lower-dimensional hyperrectangle and V, V ′ the linear spans of at most d − 2 of the unit vectors e 1 , . . . e d−1 ∈ R d−1 . Moreover, with the node patch π γ (z) := Γ m : z ∈ Γ m ∧ m = 1, . . . , M for z ∈ N γ , we suppose the following compatibility assumption for the different parametrizations: For all nodes z ∈ N γ , there exists a polytope π γ (z) ⊂ R d−1 , i.e., an interval for d = 2, a polygon for d = 3, and a polyhedron for d = 4, and a bi-Lipschitz mapping such that γ −1 z • γ m is an affine bijection for all m ∈ {1, . . . , M} with Γ m ⊆ π γ (z). Put into words, this means that each node patch π γ (z) can be flattened, where the corresponding bi-Lipschitz mapping restricted to any contained surface Γ m essentially coincides with the parametrization γ m . In particular, this prohibits the case π γ (z) = Γ. The same assumption is also made in [SS11, Assumption 4.3.25] for curvilinear triangulations. It particularly implies that the parametrizations essentially coincide at the boundary of the surfaces, i.e., for all m = m ′ with non-empty intersection E := Γ m ∩ Γ m ′ = ∅, it holds that frag replacements x 1 x 2 x 3 x 1 x 2 x 1 PSfrag replacements x 1 x 2 x 3 x 1 x 2 x 1 PSfrag replacements x 1 x 2 x 3 x 1 x 2 x 1 x 2 Figure 3.1. A mesh Γ m : m = 1, . . . , 24 of some boundary Γ is depicted, where the patch π γ (z) for z = (1/20, 1/10) is highlighted in gray (left). The corresponding patch π γ (z) is depicted (middle). The mapping γ z maps the light, medium, and dark gray affine elements to the respective elements on Γ.
The parameter domain Γ m of the light gray element is depicted (right). The mapping γ −1 z • γ m maps Γ m affinely onto the light gray element in π γ (z). where By possibly enlarging C γ from (3.2), we can assume that for (in j) monotonously increasing knots t i(0,m),j ∈ [0, 1] such that 0 = t i(0,m),0 = · · · = t i(0,m),p i,m < t i(0,m),p i,m +1 and t i(0,m),N i,m −1 < t i(0,m),N i,m = · · · = t i(0,m),N i,m +p i,m = 1. Additionally, we suppose for the multiplicity (within the knot vector K i(0,m) ) of all interior knots t i(0,m),j ∈ (0, 1) that (3.12) The so-called B-splines form a basis of the latter space and the tensor-product splines A corresponding basis is given by the set of tensor-product B-splines For i ∈ {1, . . . , d − 1}, we set K i(uni(0),m) := K i(0,m) and recursively define K i(uni(k+1),m) for k ∈ N 0 as the uniform h-refinement of (3.17) i.e., it is obtained by inserting the midpoint of any non-empty knot span [t i(uni(k),m),j , t i(uni(k),m),j+1 ] with multiplicity one to the knots K i(uni(k),m) . In particular, the old knots K i(uni(k),m) are not changed in K i(uni(k+1),m) . Replacing the index 0 by uni(k), one can define the corresonding meshes, splines, and B-splines ,m as before. This yields a nested sequence of tensor-product spline spaces where the last relation follows from the assumption for the multiplicity of the interior knots.
For an illustrative example of a hierarchical mesh, see Figure 3.3. In particular, any uniformly refined tensor mesh T uni(k),m with k ∈ N 0 is a hierarchical mesh. For a hierarchical mesh T •,m , we define a corresponding nested sequence ( Ω k It holds that •,m are highlighted in black, red, blue and green. Indeed, in the literature, one usually assumes that one is given the sequence ( Ω k •,m ) k∈N 0 and defines the corresponding hierarchical mesh via (3.22). Note that, for The space of (D-dimensional) hierarchical splines on the parameter domain Γ m = [0, 1] d−1 is just defined as the product of the resulting span It is easy to check that if T •,m is a tensor mesh and hence coincides with some T uni(k),m , then the hierarchical basis and the standard tensor-product B-spline basis are the same. Thus, the notation is consistent with the notation from (3.16). One can prove that B •,m is linearly independent; see, e.g., [VGJS11, •,m . The hierarchical basis B •,m and the mesh T •,m are compatible in the sense that (3.26) see, e.g., [GHP17, Section 3.2] for the elementary argumentation. Finally, we say that a hierarchical mesh T •,m is finer than T •,m if T •,m is obtained from T •,m via iterative dyadic bisection. Formally, this can be stated as Ω k •,m ⊆ Ω k •,m for all k ∈ N 0 . In this case, the corresponding hierarchical spline spaces are nested according to [VGJS11, Proposition 4], i.e., (3.27) For a more detailed introduction to hierarchical meshes and splines, we refer to, e.g., [VGJS11,GJS12,BG16a,SM16].
3.4. Admissible hierarchical meshes on the parameter domains. We recall the definition and some properties of admissible meshes introduced in [GHP17, Section 3.5]. Let m ∈ {1, . . . , M} and T •,m be an arbitrary hierarchical mesh on Γ m . We define the set of all neighbors of an element T ∈ T •,m as (3.28) The following lemma from [GHP17, Lemma 5.2] provides a relation between the set of neighbors and the patch of an element T ∈ T •,m .
Lemma 3.2. There holds that Let T m be the set of all admissible hierarchical meshes on Γ m . Clearly, all tensor meshes T uni(k),m , k ∈ N 0 , belong to T m . Moreover, admissible meshes satisfy the following interesting properties, which are also important for an efficient implementation of finite or boundary element methods with hierarchical splines; see [GHP17, Proposition 3.2] for a proof. [BG16b] studied related admissible hierarchical meshes. Whereas our admissibility is chosen such that Proposition 3.3 holds for hierarchical B-splines, their weaker definition of admissibility only guarantees Proposition 3.3 for truncated hierarchical B-splines; see [GHP17, Remark 3.2]. We also note that our admissibility guarantees Proposition 3.3 for both standard and truncated hierarchical B-splines. This follows immediately from (5.21) below. Moreover, the computation of the truncated hierarchical B-splines simplifies greatly on our admissible meshes; see [GHP17, Proposition 5.2].
Remark 3.5. For T •,m ∈ T m , Lemma 3.2 immediately yields that (3.31) We note that in the proof of Lemma 3.2, one exploits that all interior knot multiplicities are smaller or equal than the corresponding polynomial degree p i,m . Actually, this is the only place, where we need this assumption, since we do not require continuous ansatz functions to discretize H −1/2 (Γ) D . However, this lemma is of course essential as it implies for example local quasi-uniformity (3.31) for admissible meshes. If one drops the additional assumption on the knot multiplicities and allows multiplicities up to p i,m + 1 as well as lowest-order polynomial degrees p i,m = 0, one could define the neighbors of an element T in an arbitrary hierarchical mesh T •,m differently as and adapt the definition of admissibility. Then, Lemma 3.2 and Proposition 3.3 remain valid. Moreover, newly inserted knots can also have multiplicity p i,m + 1, i.e., the choice q i,m = p i,m + 1 in Remark 3.1 is possible.
Finally, it remains to define hierarchical meshes and splines on Γ itself. This can be done by gluing the previous definitions for the single surfaces Γ m together. For all m ∈ {1, . . . , M}, let T •,m be a hierarchical mesh on Γ m . We define the corresponding hierarchical mesh on Γ by . We call T • admissible if the mesh satisfies the following two properties: • All partial meshes are admissible, i.e., T •,m ∈ T m for all m ∈ {1, . . . , M}.
• There are no hanging nodes on the boundary of the surfaces Γ m , i.e., the intersection T ∩T ′ for T ∈ T •,m , T ′ ∈ T •,m ′ with m = m ′ is either empty or a common (transformed) lower-dimensional hyperrectangle. We define the set of all admissible hierarchical meshes on Γ as T, and suppose that the initial mesh on Γ is admissible, i.e., (3.34) For an arbitrary hierarchical mesh T • on Γ, the corresponding hierarchical splines read Remark 3.6. (a) The property that there are no hanging nodes implies local quasi-uniformity at the boundaries ∂Γ m , i.e., if Since, the ansatz space only needs to be a subset of H −1/2 (Γ) D , the property that there are no hanging nodes can actually be replaced by local quasi-uniformity at the boundaries ∂Γ m as in (a), which is sufficient for the following analysis. However, for the hyper-singular integral equation which appears for the Neumann problem Pu = 0 in Ω with conormal derivative D ν u = φ on Γ for some given φ ∈ H −1/2 (Γ) D (see, e.g., [McL00, pages 229-231]), the ansatz functions must be in H 1/2 (Γ) D . In this case, the natural choice for the ansatz space is X • ∩ C 0 (Γ) D , which is even a subset of H 1 (Γ) D . If one supposes that there are no hanging nodes, one can define a local basis with the help of [GHP17, Proposition 3.1], which states that hierarchical B-splines on The knowledge of such a basis is not only essential for an efficient implementation, but is also needed, e.g., for the definition of a quasi-interpolation operator.
3.6. Refinement of hierarchical meshes. In this section, we present a concrete refinement algorithm to specify the setting of Section 5.2. We start in the parameter domain. Recall that a hierarchical mesh T •,m is finer than another hierarchical mesh T •,m if Ω k •,m ⊆ Ω k ; see (3.42) below. In this case, the corresponding spaces are nested according to (3.27). To transfer this definition onto the surface Γ m for m ∈ {1, . . . , M}, we essentially just drop the symbol ·. We say that a hierarchical mesh T •,m on Γ m is finer than another hierarchical mesh T •,m on Γ m , if the corresponding meshes in the parameter domain satisfy this relation, i.e., if T •,m is finer than T •,m . In this case, it holds that (3.36) Finally, we call a hierarchical mesh T • on Γ finer than another hierarchical mesh T • on Γ, if the corresponding partial meshes satisfy this relation, i.e., if T •,m is finer than T •,m for all m ∈ {1, . . . , M}. In this case, it holds that Let T • be a hierarchical mesh and T ∈ T • , i.e., T ∈ T •,m for some m ∈ {1, . . . , M}. Let T = γ −1 m (T ) be the corresponding element in the parameter domain. We define the sets of neighbors is the set from (3.28), and Further, we define the sets of bad neighbors Clearly, refine(T • , M • ) is finer than T • . For any hierarchical mesh T • on Γ, we define refine(T • ) as the set of all hierarchical meshes T • on Γ such that there exist hierarchical meshes T (0) , . . . , T (J) and marked elements M (0) , . . . , The following proposition characterizes the set refine(T • ). In particular, it shows that refine(T 0 ) = T. The proof works as for IGAFEM in [GHP17, Proposition 5.1], where the assertion is stated for Γ = Γ 1 ⊂ R d−1 ; see [Gan17, Proposition 5.4.3] for further details.
Proposition 3.8. If T • ∈ T, then refine(T • ) coincides with the set of all admissible hierarchical meshes T • that are finer than T • . 3.7. Error estimator. Let T • ∈ T. Due to the regularity assumption f ∈ H 1 (Γ) D , the mapping property (2.14), and X • ⊂ L 2 (Γ) D , the residual satisfies that f − VΨ • ∈ H 1 (Γ) D for all Ψ • ∈ X • . This allows to employ the weighted-residual a posteriori error estimator where, for all T ∈ T • , the local refinement indicators read This estimator goes back to the works [CS96,Car97], where reliability (3.48) is proved for standard 2D BEM with piecewise polynomials on polygonal geometries, while the corresponding result for 3D BEM is first found in [CMS01].
3.9. Optimal convergence. We set and, for all s > 0, We say that the solution φ ∈ H −1/2 (Γ) D lies in the approximation class s with respect to the estimator if φ A est s < ∞. (3.47) By definition, φ A est s < ∞ implies that the error estimator η • on the optimal meshes T • decays at least with rate O (#T • ) −s . The following main theorem of our work states that each possible rate s > 0 will indeed be realized by Algorithm 3.9. The proof is given in Section 5. There, we will verify certain abstract properties to employ our recent abstract counterpart [GP20a, Theorem 3.4] of Theorem 3.10, which itself builds on the so-called axioms of adaptivity from [CFPP14]. Such an optimality result was first proved in [FKMP13] for the Laplace operator P = −∆ on a polyhedral domain Ω. As ansatz space, [FKMP13] considered piecewise constants on shape-regular triangulations. [FFK + 14] in combination with [AFF + 17] extends the assertion to piecewise polynomials on shape-regular curvilinear triangulations of some piecewise smooth boundary Γ. Independently, [Gan13] proved the same result for globally smooth Γ and arbitrary symmetric and elliptic boundary integral operators.
Theorem 3.10. Let (T ℓ ) ℓ∈N 0 be the sequence of meshes generated by Algorithm 3.9. Then, there hold the following statements (i)-(iii): (i) The residual error estimator (3.43)satisfies reliability, i.e., there exists a constant C rel > 0 such that (3.48) (ii) For arbitrary 0 < θ ≤ 1 and C min ∈ [1, ∞], the estimator converges linearly, i.e., there exist constants 0 < ρ lin < 1 and C lin ≥ 1 such that There exists a constant 0 < θ opt ≤ 1 such that for all 0 < θ < θ opt and C min ∈ [1, ∞), the estimator converges at optimal rate, i.e., for all s > 0, there exist constants c opt , C opt > 0 such that All involved constants C rel , C lin , q lin , θ opt , and C opt depend only on the dimensions d, D, the coefficients of the differential operator P, the boundary Γ, the (fixed) number M of boundary parts Γ m , the parametrizations γ m and γ z , the initial meshes T 0,m , and the polynomial orders Remark 3.11. If the sesquilinear form (V · ; ·) is Hermitian, then C lin , ρ lin , and C opt are independent of (Φ ℓ ) ℓ∈N 0 ; see [GP20a, Remark 4.14].
Remark 3.12. Under the assumption that h ℓ L ∞ (Ω) → 0 as ℓ → ∞, one can show that for the elementary proof. This observation allows to follow the ideas of [BHP17] and to show that the adaptive algorithm yields optimal convergence even if the sesquilinear form (V · ; ·) is only elliptic up to some compact perturbation, provided that the continuous problem is well-posed. This includes, e.g., adaptive BEM for the Helmholtz equation; see [Ste08, Section 6.9]. For details, the reader is referred to [BHP17,BBHP19].We note that h ℓ L ∞ (Ω) → 0 can easily be algorithmically enforced (without influencing the optimality result) by additionally marking the largest elements in Algorithm 3.9 (iii); see [BHP17, Proposition 16].
Remark 3.13. (a) Theorem 3.10 is still valid if one replaces the ansatz space X • by rational hierarchical splines, i.e., by the set We will prove this generalization in Section 5.4. In this case, the constants depend additionally on W 0 . (b) Moreover, Theorem 3.10 still holds true if newly inserted knots have a multiplicity higher than one, i.e., if one uses the uniformly refined knots of Remark 3.1 with 1 ≤ q i,m ≤ p i,m to define (rational) hierarchical splines. The corresponding proof works verbatim.
(c) Finally, if one defines for an element T of a hierarchical mesh T •,m its neighbours N •,m ( T ) as in Remark 3.5, and adapts the definition of admissibility and refine(·, ·) accordingly, one can also allow for lowest-order polynomial degrees p i,m ∈ N 0 as well as full knot multiplicities q i,m = p i,m + 1.

Numerical experiments with hierarchical splines
In this section, we empirically investigate the performance of Algorithm 3.9 in two typical situations: In Section 4.1, the solution is generically singular at the edges of Γ = ∂Ω. In Section 4.2, the solution is nearly singular at one point. We consider the 3D Laplace operator P := −∆. The corresponding fundamental solution reads The integral representation (4.5) is satisfied for both examples. Indeed, the surfaces Γ m of the boundary Γ = M m=1 Γ m are parametrized via rational splines, i.e., for each m ∈ {1, . . . , M}, there exist polynomial orders p 1(γ,m) , p 2(γ,m) ∈ N, a two-dimensional vector K γ,m = ( K 1(γ,m) , K 2(γ,m) ) of p i(γ,m) -open knot vectors with multiplicity smaller or equal to p i(γ,m) for the interior knots, and a positive spline weight function W γ,m ∈ S (p 1(γ,m) ,p 2(γ,m) ) ( K γ,m ) such that the parametrization γ m : Γ m → Γ m satisfies that Based on the knots K γ,m for the geometry, we choose the initial knots K 0,m for the discretization. As basis for the considered ansatz spaces of (non-rational) hierarchical splines, we use the basis given in (5.24). To (approximately) calculate the Galerkin matrix and the right-hand side vector, we proceed as in [SS11,Chapter 5] where all singular integrals are transformed via Duffy transformations and then computed with tensor Gauss quadrature. To calculate the weighted-residual error estimator 3 (3.43), we employ formula (2.2) for the surface gradient and use again tensor Gauss quadrature. To this end, we approximate ∇((f − VΦ ℓ ) • γ m ) on an element T ∈ T ℓ,m by the gradient of the polynomial interpolation of the residual f − VΦ ℓ as in [Kar12, Section 7.1.5]. In particular, we have to evaluate the residual at some quadrature points which can be done (approximately) using appropriate Duffy transformations and tensor Gauss quadrature as in [Gan14, Sections 5.1-5.2].
To (approximately) calculate the energy error, we proceed as follows: Let Φ ℓ ∈ X ℓ be the Galerkin approximation of the ℓ-th step with the corresponding coefficient vector c ℓ . Further, let V ℓ be the Galerkin matrix. With Galerkin orthogonality (2.20) and the energy norm φ 2 V = Vφ ; φ obtained by Aitken's ∆ 2 -extrapolation, we can compute the energy error as  This means that the considered integral equation stems from an exterior Laplace-Dirichlet problem (4.2). In particular, we expect singularities at the non-convex edges of R 3 \ Ω, i.e., at all edges of the cube Ω. We choose the parameters of Algorithm 3.9 as θ = 0.5 and C min = 1, where we use the refinement strategy of Remark 3.13 (c) in the lowest-order case p = 0. For comparison, we also consider uniform refinement, where we mark all elements in each step, i.e., M ℓ = T ℓ for all ℓ ∈ N 0 . This leads to uniform bisection of all elements. In Figure 4.1, one can see some adaptively generated hierarchical meshes. In Figure 4.2 and Figure 4.3, we plot the energy error φ−Φ ℓ V and the error estimator η ℓ against the number of elements #T ℓ . All values are plotted in a double-logarithmic scale. The experimental convergence rates are thus visible as the slope of the corresponding curves. Although we only proved reliability (3.48) of the employed estimator, the curves for the error and the estimator are parallel in each case, which numerically indicates reliability and efficiency. The uniform approach always leads to the suboptimal convergence rate O((#T ℓ ) −1/3 ) due to the edge singularities. Independently on the chosen polynomial degree p, the adaptive approach leads approximately to the rate O((#T ℓ ) −1/2 ). For smooth solutions φ, one would expect the rate O((#T ℓ ) −3/4−p/2 ); see [SS11, Corollary 4.1.34]. However, according to Theorem 3.10, the achieved rate is optimal if one uses the proposed refinement strategy and the resulting hierarchical splines. The reduced optimal convergence rate is probably due to the edge singularites. A similar convergence behavior is also witnessed in [FL07, Section 5.2] for the lowest-order case p = 0. [FL07] additionally considers anisotropic refinement which recovers the optimal convergence rate O((#T ℓ ) −3/4 ). We mention that for graded meshes the theoretical convergence rates are explicitly known on polyhedral domains; see [VPS90] or [GS18, Chapter 7].
We prescribe the exact solution of the interior Laplace-Dirichlet problem (4.6) as the shifted fundamental solution with y 0 := 10 −1 (0.95·2 −3/2 , 0.95·2 −3/2 , 1/2) ∈ R 3 \Ω. Although u is smooth on Ω, it is nearly singular at the midpoint y 0 := 10 −1 (2 −3/2 , 2 −3/2 , 1/2) of Γ 1 . We consider the corresponding integral equation (4.7). The normal derivative φ = ∂ ν u of u reads (4.14)   . . 0, 1, . . . , 1) for all m ∈ {1, . . . , 6}, where the multiplicity of 0 and 1 is p + 1. We choose the parameters of Algorithm 3.9 as θ = 0.5 and C min = 1, where we use the refinement strategy of Remark 3.13 (c) in the lowest-order case p = 0. For comparison, we also consider uniform refinement, where we mark all elements in each step, i.e., M ℓ = T ℓ for all ℓ ∈ N 0 . In Figure 4.4, one can see some adaptively generated hierarchical meshes. In Figure 4.5 and Figure 4.6, we plot the energy error φ − Φ ℓ V and the error estimator η ℓ against the number of elements #T ℓ . All values are plotted in a double-logarithmic scale and the experimental convergence rates are visible as the slope of the corresponding curves. In all cases, the lines of the error and the error estimator are parallel, which numerically indicates reliability and efficiency. Since the solution φ is smooth, the uniform and the adaptive approach both lead to the optimal asymptotic convergence rate O((#T ℓ ) −3/4−p/2 ). However, φ is nearly singular at y 0 , which is why adaptivity yields a much better multiplicative constant.  Uniform (for p = 2) and adaptive (θ = 0.5 for p ∈ {0, 1, 2}) refinement is considered.

Proof of Theorem 3.10
In [GP20a, Section 3], we have identified the crucial properties of the underlying meshes (see (M1)-(M5) below), the mesh-refinement (see (R1)-(R3) below), as well as the boundary element spaces (see (S1)-(S6) below) which ensure that the residual error estimator fits into the general framework of [CFPP14] and which hence guarantee reliability 3.48 and linear convergence at optimal rate (3.49)-(3.50); see [GP20a,Theorem 3.4]. The goal of this section is to recapitulate these properties and subsequently verify them for the present isogeometric setting; see also [Gan17, Section 5.2-5.3] for more details. In Section 5.4, we further briefly comment on how to verify the mentioned properties for rational hierarchical splines; see Remark 3.13. In order to ease notation, we introduce for any mesh T • the corresponding mesh-width function For ω ⊆ Γ, we define the patches of order q ∈ N 0 inductively by The corresponding set of elements is To abbreviate notation, we set π • (ω) := π 1 • (ω) and Π • (ω) := Π 1 • (ω). If ω = {z} for some z ∈ Γ, we write π q • (z) := π q • ({z}) and Π q • (z) := Π q • ({z}), where we skip the index for q = 1 as before. For S ⊆ T • , we define π q • (S) := π q • ( S) and Π q • (S) := Π q • ( S), and the superscript is omitted for q = 1.
In the following Sections (5.1.1)-(5.1.5), we show the existence of constants C patch , C locuni , C shape , C cent , C semi > 0 such that the following properties are satisfied for all admissible hierarchical meshes T • ∈ T: (M1) Bounded element patch: For all T ∈ T • , there holds that i.e., the number of elements in a patch is uniformly bounded.
The following proposition shows that (M5) is actually always satisfied. However, in general the multiplicative constant depends on the shape of the point patches. Proposition 5.1. Let ω ⊂ R d−1 be a bounded and connected Lipschitz domain and γ ω : ω → ω ⊆ Γ bi-Lipschitz. In particular, there exists a constant C lipref > 0 such that Then, there exists a constant C semi ( ω) > 0 such that The constant C semi ( ω) > 0 depends only on the dimension d, the set ω, and C lipref .
Step 1: According to Remark 3.5 and Remark 3.6, admissibility T • ∈ T shows that |level(T ) − level(T ′ )| ≤ 1 for all T ′ ∈ Π • (T ). Since we only use dyadic bisection, there exists an upper bound of the number of possible configurations of T and π • (T ) depending only on the initial meshes T 0,m ′ and (an upper bound of) level(T ). In particular, this implies that diam(T ) dist(T, Γ \ π • (T )), but the hidden constant still depends on (an upper bound for) level(T ). We see that it only remains to consider small elements T with high level.
5.2. Mesh-refinement. In the following Section 5.2.1, we show that there exist C child ≥ 2 and 0 < ρ child < 1 such that for all admissible hierarchical meshes T • ∈ T and refinements T • := refine(T • , M • ) with arbitrary marked elements M • ⊆ T • , the following elementary properties (R1)-(R3) are satisfied: (R1) Child estimate: It holds that i.e., one step of refinement leads to a bounded increase of elements. (R2) Parent is union of children: For all T ∈ T • , it holds that i.e., each element T is the union of its successors. (R3) Reduction of children: For all T ∈ T • , it holds that i.e., successors are uniformly smaller than their parent. By induction and the definition of refine(T • ), one easily sees that (R2)-(R3) remain valid if T • is an arbitrary mesh in refine(T • ). In particular, (R2)-(R3) imply that each refined element T ∈ T • \ T • is split into at least two children, which yields that (5.9) Besides (R1)-(R3), we show in Section 5.2.2-5.2.3 the following less trivial requirements (R4)-(R5) with constants C clos , C over > 0: (R4) Closure estimate: Let (T ℓ ) ℓ∈N 0 be an arbitrary sequence in T such that T ℓ+1 = refine(T ℓ , M ℓ ) with some M ℓ ⊆ T ℓ for all ℓ ∈ N 0 . Then, for all ℓ ∈ N 0 , it holds that Summing all components gives the overlay estimate It remains to show that T • is a refinement of T • and T ⋆ , i.e., T • ∈ refine(T • ) ∩ refine(T ⋆ ). Clearly, T • is finer than these meshes. Since Proposition 3.8 states that any finer admissible mesh can be reached via refine, we just have to verify admissibility of T • . As each T •,m is admissible, we only have to show that T ∩ T ′ is a common transformed lower-dimensional hyperrectangle for all T, T ′ ∈ T • with non-empty intersection and T ⊆ Γ m as well as T ′ ⊆ Γ m ′ for some m, m ′ ∈ {1, . . . , M} with m = m ′ . Without loss of generality, we assume that T ∈ T • and T ′ ∈ T ⋆ . Further, we may assume that dim(T ∩ T ′ ) > 0. Then, by definition of T • , there exist T ⋆ ∈ T ⋆ and T ′ • ∈ T • with T ⊆ T ⋆ ⊆ Γ m and T ′ ⊆ T ′ • ⊆ Γ m ′ . Obviously, T and T ′ • have non-empty intersection too. Thus, admissibility of T • shows that T ∩ T ′ • is a common transformed hyperrectangle. We suppose that T ′ is obtained from T ′ • via iterative bisections, i.e., T ′ T ′ • , and lead this to a contradiction. The intersection T ∩ T ′ is only a proper subset of a transformed hyperrectangle of T . Since T ⋆ ⊇ T , the same holds for T ⋆ instead of T , i.e., T ⋆ ∩ T ′ is only a proper subset of a transformed hyperrectangle of T ⋆ . Thus, admissibility of T ⋆ leads to a contradiction, and we see that T ′ = T ′ • shares a common transformed hyperrecangle with T . 5.3. Boundary element space. In the following Sections 5.3.1-5.3.5, we show the existence of constants C inv > 0, q loc , q proj , q supp ∈ N 0 , and 0 < ρ unity < 1, such that the following properties (S1)-(S4) hold for all admissible hierarchical meshes T • ∈ T with associated hierarchical spline space X • : (S1) Inverse inequality: For all Ψ • ∈ X • , it holds that Remark 5.2. Clearly, (S4) is in particular satisfied if the considered ansatz space X • is a product space, i.e., X • = D j=1 (X • ) j , and each component (X • ) j ⊂ L 2 (Γ) satisfies (S4). Besides (S1)-(S4), we show in Section 5.3.6 that there exist constants C sz > 0 as well as q sz ∈ N 0 such that for all admissible hierarchical meshes T • ∈ T with associated hierarchical spline space X • and S ⊆ T • , there exists a linear operator J •,S : L 2 (Γ) D → Ψ • ∈ X • : Ψ • | (T•\S) = 0 with the following properties (S5)-(S6): (S5) Local projection property. Let q loc , q proj ∈ N 0 from (S3). For all ψ ∈ L 2 (Γ) D and • (T )) . 5.3.1. Verification of (S1). For piecewise constants and piecewise affine functions on a triangulation of the boundary of a polyhedral domain Ω, the inverse estimate (S1) is already found in [DFG + 04, Theorem 4.7]. [GHS05, Theorem 3.6] and [Geo08, Theorem 3.9] generalized the result to arbitrary piecewise polynomials on curvilinear triangulations. In the recent own work [FGHP17, Proposition 4.1] and based on the ideas of [DFG + 04], we proved (S1) for non-rational splines on a one-dimensional piecewise smooth boundary Γ. In the proof, we derived the following abstract criterion for the ansatz functions which is sufficient for the inverse inequality (S1). Although, [FGHP17] considered only d = 2, the proof works verbatim for arbitrary dimension d ≥ 2.
Proposition 5.3. Let T • ∈ T be a general mesh as in Section 5.1 which satisfies (M1)-(M5). We assume that the Lipschitz constants of the mappings γ T : T → T are uniformly bounded, i.e., there exists a constant C lip > 0 such that Moreover, let ψ ∈ L 2 (Γ) satisfy the following assumption: There exists a constant ρ inf ∈ (0, 1), such that for all T ∈ T • there exists a hyperrectangular subset R T of the interior T • (i.e., R T has the form for some real numbers a T,i < b T,i ) such that |R T | ≥ ρ inf |T |, ψ does not change its sign on R T , and (5.14) where essinf denotes the essential infimum. We further assume that the shape-regularity constants of the sets R T are uniformly bounded, i.e., there exists a constant C rec > 0 such that Then, there exists a constant C inv > 0 such that The constant C inv depends only on d, C lip , ρ inf , C rec , and (M1)-(M5).
To apply Proposition 5.3 to hierarchical splines, we need the next elementary lemma which was proved in the recent own work [FGHP17, Proposition 4.1].
Lemma 5.4. Let p ∈ N 0 be a fixed polynomial degree. Let I ⊂ R be a compact interval with positive length |I| > 0. Then, there exists a constant ρ ∈ (0, 1) such that for all polynomials P of degree p on I, there exists some interval |P (t)| ≥ ρ P L ∞ (I) . (5.17) The constant ρ depends only on p.
Finally, we come to the proof of the inverse inequality (S1). Let T • ∈ T be an admissible hierarchical mesh on Γ. We recall that X • is a product space of transformed scalar-valued hierarchical splines. Without loss of generality, we assume that we are in the scalar case D = 1. We show that all Ψ • ∈ X • ⊂ L 2 (Γ) satisfy the assumptions of Proposition 5.3 and hence conclude that h 1/2 . We have already seen that (M1)-(M5) are satisfied. Moreover, (5.13) is trivially satisfied since each γ T is just the restriction of some γ m to T = γ −1 m (T ), where m ∈ {1, . . . , M}. For T ∈ T • , we abbreviate Ψ • := Ψ • • γ T . Due to the regularity (3.2) of the parametrizations γ m , it is sufficient to find a uniform constant ρ inf ∈ (0, 1) and a shape-regular hyperrectangle R T ⊂ T • such that | R T | ≥ ρ inf | T |, Ψ • does not change sign on R T , and (5.18) Indeed, one easily shows that | R T | ≥ ρ inf | T | implies that |R T | ≥ ρ inf |T | for some uniform constant ρ inf ∈ (0, 1); see, e.g., [GHP17, Section 5.3]. Recall that Ψ • coincides with a tensor-product polynomial P . Hence, there exist polynomials P i of degree p max such that We define ( R T ) i as the interval of Lemma 5.4 corresponding to the polynomial P i on the interval I = T i . With the constant ρ of Lemma 5.4, we set ρ inf := ρ d−1 . Then, (5.19), and therefore (5.18) is satisfied. Moreover, one sees that | R T | ≥ ρ inf | T |, and that Ψ • doesn't change its sign on R T ⊂ T • . It remains to prove shape-regularity (5.15). Since, the refinement procedure refine only uses uniform bisection of elements, the element T is shape-regular in the sense that | T i | ≃ | T i ′ | for all i, i ′ ∈ {1, . . . , d − 1}. Together with |( R T ) i | ≥ ρ| T i |, this proves (5.15). Altogether, we conclude (S1), where the constant C inv depends only on the dimensions d, D, the (fixed) number M of boundary parts Γ m , the parametrizations γ m and γ z , the initial meshes T 0,m , and the polynomial orders (p 1,m , . . . , p d−1,m ) for m ∈ {1, . . . , M} and z ∈ N γ . 5.3.2. Verification of (S2). Let T • ∈ T and T • ∈ refine(T • ). The nestedness X • ⊆ X • was already stated in (3.37).

5.3.3.
Bases of hierarchical splines on the boundary. In this section, we give two bases for X • . Since X • is a product space of transformed scalar-valued hierarchical splines, it is sufficient to consider D = 1. An alternate basis of X • is given by span Trunc • (β) : β ∈ B • . 5.3.4. Verification of (S3). With Remark 3.4(a) and (5.24), the proof of (S3) with q loc := q proj + 2(p max + 1) works as for IGAFEM in [GHP17, Section 5.8], where the assertion is stated for Γ = Γ 1 ⊂ R d−1 ; see [Gan17, Proposition 5.5.12] for further details. 5.3.5. Verification of (S4). Let T • ∈ T. First, we recall that X • is a product space of transformed scalar-valued hierarchical splines. Without loss of generality, we assume that D = 1; see Remark 5.2. [Fae00, Lemma 2.6] resp. [Fae02, Lemma 3.5] prove a similar version of (S4) for splines on a one-dimensional boundary Γ resp. for certain piecewise polynomials of degree 0, 1, 5, and 6 on curvilinear triangulations of a two-dimensional boundary Γ. There, the proof follows from direct calculations, where [Fae00, Lemma 2.6] actually only provides a proof of (S4) for splines of degree 2. In contrast, we will make use of the following abstract result.
5.4. Proof of Theorem 3.10 for rational hierarchical splines. As mentioned in Remark 3.13, Theorem 3.10 is still valid if one replaces the ansatz space X • for T • ∈ T by rational hierarchical splines, i.e., by the set where W 0,m = W 0 • γ m ∈ S (p 1,m ,...,p d−1,m ) ( K 0,m ) is a fixed positive weight function in the initial space of hierarchical splines for all m ∈ {1, . . . , M}, where we additionally assume the representation (3.52). Indeed, the mesh properties (M1)-(M5) as well as the refinement properties (R1)-(R5) of Section 5 are independent of the discrete spaces. To verify the validity of Theorem 3.10 in the rational setting, it thus only remains to verify the properties (S1)-(S6) for the rational boundary element spaces. To see the inverse estimate (S1), it is again sufficient to consider D = 1. In Section 5.3.1, we proved (S1) for X • by applying Proposition 5.3 for all Ψ • ∈ X • . With the notation from Section 5.3.1, we showed that where Ψ • does not change its sign on R T . With 0 < w min := inf x∈Γ W 0 (x), w max := sup x∈Γ W 0 (x), and ρ inf := ρ inf w min /w max , this yields that, for all Ψ • ∈ X • , In particular, the conditions for Proposition 5.3 are also satisfied for the functions in X W 0 • , which concludes (S1).
The properties (S2)-(S3) depend only on the numerator of the rational hierarchical splines and thus transfer.
For the proof of (S4), we exploit the representation (3.52) to verify the conditions of the abstract Proposition 5.5. Again, we may assume that D = 1. Let T • ∈ T. Note that W 0,m is also an element of the standard tensor-product spline space S (p 1,m ,...,p d−1,m ) ( K uni(k),m ) for all m ∈ {1, . . . , M} and k ∈ N 0 . In particular, it can be written as linear combination of B-splines in B uni(k),m . The representation (3.52) and the two-scale relation with only nonnegative coefficients between bases of consecutive levels of [dB01, Section 11] yields that the corresponding coefficients are non-negative. Therefore, the preservation of coefficients property of [SM16, Theorem 1] or [GJS14,Theorem 12] implies that also the coefficients of the linear combination of W 0,m in Trunc •,m ( β) : β ∈ B •,m are non-negative, i.e., As in Section 5.3.5, one sees that this choice satisfies the assumptions of Proposition 5.5. To see (S5) and (S6), we define the corresponding projection operator The desired properties transfer immediately from the non-rational case.