Optimal additive Schwarz methods for the hp -BEM: the hypersingular integral operator in 3D on locally refined meshes
Thomas Führer, Jens Markus Melenk, Dirk Praetorius, Alexander Rieder
OOptimal additive Schwarz methods for the hp -BEM: the hypersingular integraloperator in 3D on locally refined meshes T. F¨uhrer a , J. M. Melenk b, ∗ , D. Praetorius b , A. Rieder b a Pontificia Universidad Cat´olica de Chile, Facultad de Matem´aticas, Vicu˜na Mackenna 4860, Santiago, Chile. b Technische Universit¨at Wien, Institut f¨ur Analysis und Scientific Computing, Wiedner Hauptstraße 8-10, A-1040 Vienna
Abstract
We propose and analyze an overlapping Schwarz preconditioner for the p and hp boundary element method for thehypersingular integral equation in 3D. We consider surface triangulations consisting of triangles. The conditionnumber is bounded uniformly in the mesh size h and the polynomial order p . The preconditioner handles adaptivelyrefined meshes and is based on a local multilevel preconditioner for the lowest order space. Numerical experimentson different geometries illustrate its robustness. Keywords: hp -BEM, hypersingular integral equation, preconditioning, additive Schwarz method
1. Introduction
Many elliptic boundary value problems that are solved in practice are linear and have constant (or at leastpiecewise constant) coefficients. In this setting, the boundary element method (BEM, [21, 35, 39, 31]) has estab-lished itself as an effective alternative to the finite element method (FEM). Just as in the FEM applied to thisparticular problem class, high order methods are very attractive since they can produce rapidly convergent schemeson suitably chosen adaptive meshes. The discretization leads to large systems of equations, and a use of iterativesolvers brings the question of preconditioning to the fore.In the present work, we study high order Galerkin discretizations of the hypersingular operator. This is anoperator of order 1, and we therefore have to expect the condition number of the system matrix to increase asthe mesh size h decreases and the approximation order p increases. We present an additive overlapping Schwarzpreconditioner that offsets this degradation and results in condition numbers that are bounded independently ofthe mesh size and the approximation order. This is achieved by combining the recent H / -stable decompositionof spaces of piecewise polynomials of degree p of [23] and the multilevel diagonal scaling preconditioner of [14, 16]for the hypersingular operator discretized by piecewise linears.Our additive Schwarz preconditioner is based on stably decomposing the approximation space of piecewisepolynomials into the lowest order space (i.e., piecewise linears) and spaces of higher order polynomials supportedby the vertex patches. Such stable localization procedures were first developed for the hp -FEM in [34] for meshesconsisting of quadrilaterals (or, more generally, tensor product elements). The restriction to tensor product elementsstems from the fact that the localization is achieved by exploiting stability properties of the 1D-Gauß-Lobattointerpolation operator, which, when applied to polynomials, is simultaneously stable in L and H (see, e.g., [7,eqns. (13.27), (13.28)]). This simultaneous stability raises the hope for H / -stable localizations and was pioneeredin [18] for the hp -BEM for the hypersingular operator on meshes consisting of quadrilaterals. Returning to the ∗ Corresponding author
Preprint submitted to Elsevier August 29, 2018 a r X i v : . [ m a t h . NA ] A ug p -FEM, H -stable localizations on triangular/tetrahedral meshes were not developed until [36]. The techniquesdeveloped there were subsequently used in [23] to design H / -stable decompositions on triangular meshes andthus paved the way for condition number estimates that are uniform in the approximation p for overlappingSchwarz methods for the hp -version BEM applied to the hypersingular operator. Non-overlapping additive Schwarzpreconditioners for high order discretizations of the hypersingular operator are also available in the literature, [1];as it is typical of this class of preconditioners, the condition number still grows polylogarithmically in p .Our preconditioner is based on decomposing the approximation space into the space of piecewise linears andspaces associated with the vertex patches. It is highly desirable to decompose the space of piecewise linears furtherin a multilevel fashion. For sequences of uniformly refined meshes, the first such multilevel space decompositionappears to be [43] (see also [33]). For adaptive meshes, local multilevel diagonal scaling was first analyzed in [2],where for a sequence T (cid:96) of successively refined adaptive meshes a uniformly bounded condition number for thepreconditioned system is established. Formally, however, [2] requires that T (cid:96) ∩ T (cid:96) +1 ⊂ T (cid:96) + k for all (cid:96), k ∈ N , i.e.,as soon as an element K ∈ T (cid:96) is not refined, it remains non-refined in all succeeding triangulations. While thiscan be achieved implementationally, the recent works [14, 16] avoid such a restriction by considering sequences ofmeshes that are obtained in typical h -adaptive environments with the aid of newest vertex bisection (NVB). Wefinally note that the additive Schwarz decomposition on adaptively refined meshes is a subtle issue. Hierarchicalbasis preconditioners (which are based on the new nodes only) lead to a growth of the condition number with O ( | log h min | ); see [44]. Global multilevel diagonal preconditioning (which is based on all nodes) leads to a growth O ( | log h min | ); see [28, 14].The paper is organized as follows: In Section 2 we introduce the hypersingular equation and the discretizationby high order piecewise polynomial spaces. Section 3 collects properties of the fractional Sobolev spaces includingthe scaling properties. Section 4 studies in detail the p -dependence of the condition number of the unpreconditionedsystem. The polynomial basis on the reference triangle chosen by us is a hierarchical basis of the form first proposedby Karniadakis & Sherwin, [25, Appendix D.1.1.2]; the precise form is the one from [50, Section 5.2.3]. We provebounds for the condition number of the stiffness matrix not only in the H / -norm but also in the norms of L and H . This is also of interest for hp -FEM and could not be found in the literature. Section 5 develops severalpreconditioners. The first one (Theorem 5.3) is based on decomposing the high order approximation space intothe global space of piecewise linears and local high order spaces of functions associated with the vertex patches.The second one (Theorem 5.6) is based on a further multilevel decomposition of the global space of piecewiselinears. The third one (Theorem 5.10) exploits the observation that topologically, only a finite number of vertexpatches can occur. Hence, significant memory savings for the preconditioner are possible if the exact bilinear formsfor the vertex patches are replaced with scaled versions of simpler ones defined on a finite number of referenceconfigurations. Numerical experiments in Section 6 illustrate that the proposed preconditioners are indeed robustwith respect to both h and p .We close with a remark on notation: The expression a (cid:46) b signifies the existence of a constant C > a ≤ C b . The constant C does not depend on the mesh size h and the approximation order p , but may depend onthe geometry and the shape regularity of the triangulation. We also write a ∼ b to abbreviate a (cid:46) b (cid:46) a . hp -discretization of the hypersingular integral equation Let Ω ⊂ R be a bounded Lipschitz polyhedron with a connected boundary ∂ Ω, and let Γ ⊆ ∂ Ω be an open,connected subset of ∂ Ω. If Γ (cid:54) = ∂ Ω, we assume it to be a Lipschitz hypograph, [31]; the key property needed is thatΓ is such that the ellipticity condition (2.2) holds. Furthermore, we will use affine, shape regular triangulationsof Γ, which further imposes conditions on Γ. In this work, we are concerned with preconditioning high orderdiscretizations of the hypersingular integral operator, which is given by( Du ) ( x ) := − ∂ intn x (cid:90) Γ ∂ intn y G ( x, y ) u ( y ) ds y for x ∈ Γ , (2.1)2here G ( x, y ) := π | x − y | is the fundamental solution of the 3D-Laplacian and ∂ intn y denotes the (interior) normalderivative with respect to y ∈ Γ.We will need some results from the theory of Sobolev and interpolation spaces, see [31, Appendix B]. For anopen subset ω ⊂ ∂ Ω, let L ( ω ) and H ( ω ) denote the usual Sobolev spaces. The space (cid:101) H ( ω ) consists of thosefunctions whose zero extension to ∂ Ω is in H ( ∂ Ω). (In particular, for ω = ∂ Ω, H ( ∂ Ω) = (cid:101) H ( ∂ Ω).) When thesurface measure of the set ∂ Ω \ ω is positive, we use the equivalent norm (cid:107) u (cid:107) (cid:101) H ( ω ) := (cid:107)∇ Γ u (cid:107) L ( ω ) .We will define fractional Sobolev norms by interpolation. The following Proposition 2.1 collects key propertiesof interpolation spaces that we will need; we refer to [41, 45] for a comprehensive treatment. For two Banach spaces( X , (cid:107)·(cid:107) ) and ( X , (cid:107)·(cid:107) ), with continuous inclusion X ⊆ X and a parameter s ∈ (0 ,
1) the interpolation norm isdefined as (cid:107) u (cid:107) X ,X ] s := (cid:90) ∞ t =0 t − s (cid:18) inf v ∈ X (cid:107) u − v (cid:107) + t (cid:107) v (cid:107) (cid:19) dtt . The interpolation space is given by [ X , X ] s := (cid:110) u ∈ X : (cid:107) u (cid:107) [ X ,X ] s < ∞ (cid:111) .An important result, which we use in this paper, is the following interpolation theorem: Proposition 2.1.
Let X i , Y i , i ∈ { , } , be two pairs of Banach spaces with continuous inclusions X ⊆ X and Y ⊆ Y . Let s ∈ (0 , .(i) If a linear operator T is bounded as an operator X → Y and X → Y , then it is also bounded as an operator [ X , X ] s → [ Y , Y ] s with (cid:107) T (cid:107) [ X ,X ] s → [ Y ,Y ] s ≤ (cid:107) T (cid:107) − sX → Y (cid:107) T (cid:107) sX → Y . (ii) There exists a constant C > such that for all x ∈ X : (cid:107) x (cid:107) [ X ,X ] s ≤ C (cid:107) x (cid:107) − sX (cid:107) x (cid:107) sX . We define the fractional Sobolev spaces by interpolation. For s ∈ (0 , H s ( ω ) := (cid:2) L ( ω ) , H ( ω ) (cid:3) s , (cid:101) H s ( ω ) := (cid:104) L ( ω ) , (cid:101) H ( ω ) (cid:105) s . Here, we will only consider the case s = 1 /
2. We define H − / (Γ) as the dual space of (cid:101) H / (Γ), where duality isunderstood with respect to the (continuously) extended L (Γ)-scalar product and denoted by (cid:104)· , ·(cid:105) Γ . An equivalentnorm on H / (Γ) is given by (cid:107) u (cid:107) H / (Γ) ∼ (cid:107) u (cid:107) L (Γ) + | u | H / (Γ) , where |·| H / (Γ) is given by the Sobolev-Slobodeckijseminorm (see [35] for the exact definition).We now state some important properties of the hypersingular operator D from (2.1), see, e.g., [35, 31, 21, 39].First, the operator D : (cid:101) H / (Γ) → H − / (Γ) is a bounded linear operator.For open surfaces Γ (cid:36) ∂ Ω the operator is elliptic (cid:104)
Du, u (cid:105) Γ ≥ c ell (cid:107) u (cid:107) (cid:101) H / (Γ) ∀ u ∈ (cid:101) H / (Γ) , (2.2)with some constant c ell > ∂ Ω we note that (cid:101) H / (Γ) = H / (Γ) and the operator D is still semi-elliptic, i.e. (cid:104) Du, u (cid:105) Γ ≥ c ell | u | H / (Γ) ∀ u ∈ (cid:101) H / (Γ) . Moreover, the kernel of D then consists of the constant functions only: ker( D ) = span(1).3o get unique solvability and strong ellipticity for the case of a closed surface, it is customary to introduce astabilized operator (cid:101) D given by the bilinear form (cid:68) (cid:101) Du, v (cid:69) Γ := (cid:104) Du, v (cid:105) Γ + α (cid:104) u, (cid:105) Γ (cid:104) v, (cid:105) Γ , α > . (2.3)In order to avoid having to distinguish the two cases Γ = ∂ Ω and Γ (cid:36) ∂ Ω, we will only work with the stabilizedform on (cid:101) H / (Γ) and just set α = 0 in the case of Γ (cid:36) ∂ Ω. The basic integral equation involving the hypersingularoperator D then reads: For given g ∈ H − / (Γ), find u ∈ (cid:101) H / (Γ) such that (cid:68) (cid:101) Du, v (cid:69) Γ = (cid:104) g, v (cid:105) Γ ∀ v ∈ (cid:101) H / (Γ) . (2.4)We note that in the case of the closed surface Γ = ∂ Ω, the solution of the stabilized system above is equivalentto the solution of (cid:104)
Du, v (cid:105) Γ = (cid:104) g, v (cid:105) Γ under the side constraint (cid:104) u, (cid:105) Γ = (cid:104) g, (cid:105) Γ α | Γ | . Moreover, it is well known that (cid:68) (cid:101) D · , · (cid:69) Γ is symmetric, elliptic and induces an equivalent norm on (cid:101) H / (Γ), i.e., (cid:68) (cid:101) Du, u (cid:69) Γ ∼ (cid:107) u (cid:107) (cid:101) H / (Γ) ∀ u ∈ (cid:101) H / (Γ) . Let T = { K , . . . , K N } denote a regular (in the sense of Ciarlet) triangulation of the two-dimensional manifoldΓ ⊆ ∂ Ω into compact, non-degenerate planar surface triangles. We say that a triangulation is γ -shape regular, ifthere exists a constant γ > K ∈T diam( K ) | K | ≤ γ. (2.5)Let (cid:98) K := conv { (0 , , (1 , , (0 , } be the reference triangle. With each element K we associate an affine,bijective element map F K : (cid:98) K → K . We will write P p ( (cid:98) K ) for the space of polynomials of degree p on (cid:98) K . Thespace of piecewise polynomials on T is given by P p ( T ) := (cid:110) u ∈ L (Γ) : u ◦ F K ∈ P p ( (cid:98) K ) for all K ∈ T (cid:111) . (2.6)The elementwise constant mesh width function h := h T ∈ P ( T ) is defined by ( h T ) | K := diam( K ) for all K ∈ T .Let V = { z , . . . , z M } denote the set of all vertices of the triangulation T that are not on the boundary of Γ.We define the (vertex) patch ω z for a vertex z ∈ V by ω z := interior (cid:91) { K ∈T : z ∈ K } K , (2.7)where the interior is understood with respect to the topology of Γ. For p ≥
1, define (cid:101) S p ( T ) := P p ( T ) ∩ (cid:101) H / (Γ) . (2.8)Then, the Galerkin discretization of (2.4) consists in replacing (cid:101) H / (Γ) with the discrete subspace (cid:101) S p ( T ), i.e.: Find u h ∈ (cid:101) S p ( T ) such that (cid:68) (cid:101) Du h , v h (cid:69) Γ = (cid:104) g, v h (cid:105) Γ for all v h ∈ (cid:101) S p ( T ) . (2.9)4 emark 2.2. We employ the same polynomial degree for all elements. This is not essential and done for simplicityof presentation. For details on the more general case, see [23].
After choosing a basis of (cid:101) S p ( T ), the problem (2.9) can be written as a linear system of equations, and we write (cid:101) D ph for the resulting system matrix. Our goal is to construct a preconditioner for (cid:101) D ph . It is well-known that thecondition number of (cid:101) D ph depends on the choice of the basis of (cid:101) S p ( T ), which we fix in Definition 2.4 below. Weremark in passing that the preconditioned system of Section 5.2 will no longer depend on the basis. For the matrix representation of the Galerkin formulation (2.9) we have to specify a polynomial basis on thereference triangle (cid:98) K . We use a basis that relies on a collapsed tensor product representation of the triangle andis given in [50, Section 5.2.3]. This kind of basis was first proposed for the hp -FEM by Karniadakis & Sherwin,[25, Appendix D.1.1.2]; closely related earlier works on polynomial bases that rely on a collapsed tensor productrepresentation of the triangle are [26, 9]. Definition 2.3 (Jacobi polynomials) . For coefficients α , β > − the family of Jacobi polynomials on the interval ( − , is denoted by P ( α,β ) n , n ∈ N . They are orthogonal with respect to the L ( − , inner product with weight (1 − x ) α (1 + x ) β . (See for example [25, Appendix A] or [50, Appendix A.3] for the exact definitions and a list ofimportant properties). The Legendre polynomials are a special case of the Jacobi polynomials for α = β = 0 anddenoted by (cid:96) n ( s ) := P (0 , n ( s ) . The integrated Legendre polynomials L n and the scaled polynomials are defined by L n ( s ) := (cid:90) s − (cid:96) n − ( t ) dt for n ∈ N , P S , ( α,β ) n ( s, t ) := t n P ( α,β ) n ( s/t ) , L S n ( s, t ) := t n L n ( s/t ) . (2.10)On the reference triangle, our basis reads as follows: Definition 2.4 (polynomial basis on the reference triangle) . Let p ∈ N and let λ , λ , λ be the barycentriccoordinates on the reference triangle (cid:98) K . Then the basis functions for the reference triangle consist of three vertexfunctions, p − edge functions per edge, and ( p − p − / cell-based functions:(a) for i = 1 , , the vertex functions are: ϕ V i := λ i ; (b) for m = 1 , , and an edge E m with edge vertices e , e , the edge functions are given by: ϕ E m i := (cid:114) i + 32 L S i +2 ( λ e − λ e , λ e + λ e ) , ≤ i ≤ p − (c) for ≤ i + j ≤ p − the cell based functions are: ϕ I ( i,j ) := c ij λ λ λ P S , (2 , i ( λ − λ , λ + λ ) P (2 i +5 , j (2 λ − . with c ij such that (cid:13)(cid:13)(cid:13) ϕ I ( i.j ) (cid:13)(cid:13)(cid:13) L ( (cid:98) K ) = 1 . Remark 2.5.
In order to get a basis of (cid:101) S p ( T ) we take the composition with the element mappings ϕ ◦ F K . Toensure continuity along edges we take an arbitrary orientation of the edges and observe that the edge basis functions ϕ E i are symmetric under permutation of λ e and λ e up to a sign change ( − i . . Properties of (cid:102) H / (Γ) (cid:101) H / (Γ)Several results of the present paper depend on results in [23]. Therefore we present a short summary of the mainresults of that paper in this section. In [23] the authors propose an H / -stable space decomposition on meshesconsisting of triangles. It is based on quasi-interpolation operators constructed by local averaging on elements.We introduce the following product spaces: X := (cid:89) z ∈V L ( ω z ) , X := (cid:89) z ∈V (cid:101) H ( ω z ) . The spaces L ( ω z ) and (cid:101) H ( ω z ) are endowed with the L - and (cid:101) H -norm, respectively. Proposition 3.1 (localization, [23]) . There exists an operator J : L (Γ) → (cid:16) (cid:101) S ( T ) , (cid:107)·(cid:107) L (Γ) (cid:17) × X with thefollowing properties:(i) J is linear and bounded.(ii) J | (cid:101) H (Γ) is also bounded as an operator (cid:101) H (Γ) → (cid:16) (cid:101) S ( T ) , (cid:107)·(cid:107) (cid:101) H (Γ) (cid:17) × X .(iii) If u ∈ (cid:101) S p ( T ) then each component of Ju is in (cid:101) S p ( T ) .(iv) If we write Ju =: ( u , U ) , and furthermore U z for the component of U in X corresponding to the space L ( ω z ) , then Ju represents an (cid:101) H / (Γ) − stable decomposition of u , i.e., u = u + (cid:88) z ∈V U z and (cid:107) u (cid:107) (cid:101) H / (Γ) + (cid:88) z ∈V (cid:107) U z (cid:107) (cid:101) H / ( ω z ) ≤ C (cid:107) u (cid:107) (cid:101) H / (Γ) . (3.1) The norms of J in (i)—(ii) and the constant C > in (iv) depend only on Γ and the shape regularity constant γ .Sketch of proof: The first component of J (i.e., the mapping u (cid:55)→ u ) consists of the Scott-Zhang projectionoperator, as modified in [3, Section 3.2]. The local components (i.e., the functions U z , z ∈ V ) then are based on asuccessive decomposition into vertex, edge and interior parts, similar to what is done in [36]. We give a flavor of theprocedure. Set u := u − u . In order to define the vertex parts for a vertex z , we select an element K ⊂ ω z of thepatch ω z and perform a suitable local averaging of u on that element; this averaged function u loc,K is defined on K in terms of u | K and vanishes on the edge opposite z . In order to extend u loc,K to the patch ω z and thus obtainthe function u z , we define u z by “rotating” u loc,K around the vertex z . The averaging process can be done in sucha way that for continuous functions u one has u ( z ) = u z ( z ) and that one has appropriate stability properties in L and H . The edge contributions are constructed from the function u := u − (cid:80) z ∈V u z . Let E ( T ) denote theset of interior edges of T . For an edge E ∈ E ( T ) one selects an element (of which E is an edge), averages there,and extends the obtained averaged function to the edge patch by symmetry across the edge E . In this way, thefunction u E is constructed for each edge E . It again holds for sufficiently smooth u that u ( x ) = u E ( x ) ∀ x ∈ E .For u ∈ (cid:101) H (Γ) and in turn u ∈ (cid:101) H (Γ), we have that u := u − (cid:80) E∈E ( T ) u E − (cid:80) z ∈V u z vanishes on all edges; hence u | K ∈ (cid:101) H ( K ) for all K ∈ T . The terms u z , u E , u | K can be rearranged to take the form of patch contributions U z as given in the statement of the proposition. (The decomposition is not unique.) The (cid:101) H / stability is adirect consequence of the L and H stability and interpolation properties given in Proposition 2.1, (i). We finallymention that assertion (iii) follows from the fact that the averaging operators employed at the various stages ofthe decomposition are polynomial preserving. Remark 3.2.
Independently, a decomposition similar to Proposition 3.1 was presented in [10]. (cid:101) H / ( ω z )-norm to the (cid:101) H / ( O )-normif O ⊃ ω z : Corollary 3.3.
Let z ∈ V and O be the union of some triangles of T with ω z ⊆ O ⊆ Γ . Then there exist constants c , c that depend only on O , Γ , and the γ -shape regularity of T such that for all u ∈ (cid:101) H / ( O ) with supp( u ) ⊆ ω z ,we can estimate: c (cid:107) u (cid:107) (cid:101) H / ( ω z ) ≤ (cid:107) u (cid:107) (cid:101) H / ( O ) ≤ c (cid:107) u (cid:107) (cid:101) H / ( ω z ) . Proof.
To see the second inequality, consider the extension operator E that extends the function u by 0 outsideof ω z . This operator is continuous L ( ω z ) → L ( O ) and (cid:101) H ( ω z ) → (cid:101) H ( O ), both with constant 1. ApplyingProposition 2.1, (i) to this extension operator E gives the second inequality with c = 1. The first inequality ismore involved. We start by noting that the stability assertion (3.1) of Proposition 3.1 gives u = u + (cid:80) z (cid:48) ∈V U z (cid:48) and (cid:107) u (cid:107) (cid:101) H / ( O ) + (cid:88) z (cid:48) ∈V (cid:107) U z (cid:48) (cid:107) (cid:101) H / ( ω z (cid:48) ) ≤ C (cid:107) u (cid:107) (cid:101) H / ( O ) . (3.2)The constant C depends only on the set O and the shape regularity of the triangulation, when Proposition 3.1 isapplied with Γ replaced by O . The decomposition in Proposition 3.1 is not unique, and we will now exploit thisby requiring more. Specifically, we assert that the operator J , which effects the decomposition, can be chosen suchthat, for given ω z , we have supp u ⊂ ω z and U z (cid:48) = 0 for z (cid:48) (cid:54) = z . If this can be achieved, we get u = u + U z .Since (3.2) contains a term (cid:107) u (cid:107) (cid:101) H / ( O ) , we also need to reinvestigate the stability proof. The decomposition is L - and H -stable, and maps to functions with supp u ⊂ ω z . Therefore, we can interpret the first component of J as an operator mapping to (cid:101) H ( ω z ) and apply Proposition 2.1 (i) to get: (cid:107) u (cid:107) (cid:101) H / ( ω z ) + (cid:107) U z (cid:107) (cid:101) H / ( ω z ) ≤ C (cid:107) u (cid:107) (cid:101) H / ( O ) . The triangle inequality (cid:107) u (cid:107) (cid:101) H / ( ω z ) ≤ (cid:107) u (cid:107) (cid:101) H / ( ω z ) + (cid:107) U z (cid:107) (cid:101) H / ( ω z ) then concludes the proof.It therefore remains to see that we can construct the operator J of Proposition 3.1 with the additional propertythat u = u + U z if u is such that supp u ⊂ ω z . This follows by carefully selecting the elements on which the localaveraging is done, namely, whenever one has to choose an element on which to average, one selects, if possible, anelement that is not contained in ω z . For example for vertex contributions z (cid:48) ∈ ∂ω z we make sure to use elements K (cid:48) which are not in ω z . This implies that for supp( u ) ⊆ ω z we get u z (cid:48) = 0. A similar choice is made when definingthe edge contributions.The stable space decomposition of Proposition 3.1 is one of several ingredients of the proof that the interpolationspace obtained by interpolating the space (cid:101) S p ( T ) endowed with the L -norm and the H -norm yields the space (cid:101) S p ( T ) endowed with the appropriate fractional Sobolev norm: Proposition 3.4 ([23]) . Let s ∈ (0 , and let T be a shape regular triangulation of Γ . Let p ≥ . Then: (cid:104)(cid:16) (cid:101) S p ( T ) , (cid:107)·(cid:107) L (Γ) (cid:17) , (cid:16) (cid:101) S p ( T ) , (cid:107)·(cid:107) (cid:101) H (Γ) (cid:17)(cid:105) s = (cid:16) (cid:101) S p ( T ) , (cid:107)·(cid:107) (cid:101) H s (Γ) (cid:17) . The constants implied in the norm equivalence depend only on Γ , s , and the γ -shape regularity of T . We note that such a result is clearly valid for fixed p ≥ p .7 z b ω z ∂ b ω z Figure 3.1: Example of an interior reference patch (cid:98) ω z We recall that V is the set of all inner vertices, i.e., z / ∈ ∂ Γ for all z ∈ V . We define the patch size for a vertex z ∈ V as h z := diam( ω z ) and stress that γ -shape regularity implies h z ∼ h | K for all elements T ⊆ ω z .Due to shape regularity, the number of elements meeting at a vertex is bounded. The following definition allowsus to transform the vertex patches ω z to a finite number of reference configurations. Definition 3.5 (reference patch) . Let ω z be an interior patch consisting of n triangles. We may define a Lipschitzcontinuous bijective map F z : (cid:98) ω z → ω z , where (cid:98) ω z ⊆ R is a regular polygon with n edges, see Fig. 3.1. The map F z is piecewise defined as a concatenation of affine maps from triangles comprising the regular polygon to the referenceelement (cid:98) K with the element maps F K . We note that F z ( ∂ (cid:98) ω z ) = ∂ω z . The following lemma tells us how the hypersingular integral operator behaves under the patch transformation.
Lemma 3.6.
Let u ∈ (cid:101) H / (Γ) with supp( u ) ⊆ ω z . Define (cid:98) u := u ◦ F z and the integral operator (cid:98) D as (cid:98) D (cid:98) u ( x ) := − ∂ intn x (cid:82) (cid:99) ω z ∂ intn y G ( x, y ) (cid:98) u ( y ) ds ( y ) for x ∈ (cid:98) ω z , where we treat (cid:99) ω z ⊆ R as a screen embedded in R and ∂ n x is thederivative in direction of the vector (0 , , . Then the hypersingular operator scales like (cid:104) Du, u (cid:105) ω z ∼ h z (cid:68) (cid:98) D (cid:98) u, (cid:98) u (cid:69) (cid:98) ω z , where the implied constants depend only on Γ and the γ -shape regularity of T .Proof. We will prove this in three steps:(i) (cid:104)
Du, u (cid:105) Γ ∼ (cid:107) u (cid:107) (cid:101) H / ( ω z ) ,(ii) (cid:107) u (cid:107) (cid:101) H / ( ω z ) ∼ h z (cid:107) (cid:98) u (cid:107) (cid:101) H / ( (cid:98) ω z ) ,(iii) (cid:107) (cid:98) u (cid:107) (cid:101) H / ( (cid:98) ω z ) ∼ (cid:68) (cid:98) D (cid:98) u, (cid:98) u (cid:69) (cid:98) ω z . Proof of (i):
It is well-known that D is continuous and elliptic on (cid:101) H / ( ω z ). In our case, the ellipticity constantscan be chosen independently of the patch ω z and instead only depend on Γ. To do so we embed the spaces (cid:101) H / ( ω z )into finitely many larger spaces (cid:101) H / ( O j ), where the sub-surfaces O j are open and for each z there is O j such that ω z ⊆ O j ⊆ Γ. (For the screen problem we may use the single sub-surface O := Γ, for the case of closed surfaceswe can, for example, use O j := Γ \ F j , where F j is the j -th face of the polyhedron Ω such that ω z ∩ F j = ∅ ). It8s important that these surfaces are open, since for closed surfaces Γ we do not have full ellipticity of D but onlyfor (cid:101) D , and the stabilization term has a different scaling behavior. Since u vanishes outside of ω z we can use theellipticity on (cid:101) H / ( O j ) to see (cid:104) Du, u (cid:105) Γ ∼ (cid:107) u (cid:107) (cid:101) H / ( O j ) . By Corollary 3.3 the norms on ω z and O j are equivalent,which implies the statement (i). Proof of (ii):
The scalings of the L -norm and H -seminorm (we can use the seminorm, since we are workingon (cid:101) H / of an open surface) is well-known to be (cid:107) u (cid:107) L ( ω z ) ∼ h z (cid:107) (cid:98) u (cid:107) L ( (cid:98) ω z ) , (cid:107)∇ u (cid:107) L ( ω z ) ∼ (cid:107)∇ (cid:98) u (cid:107) L ( (cid:98) ω z ) . The interpolation theorem (Proposition 2.1, (i)) then proves part (ii).
Proof of (iii):
We again use ellipticity and continuity of (cid:98) D . Since there are only finitely many reference patches,the constants can be chosen independently of the individual patches.
4. Condition number of the hp -Galerkin matrix In this section we investigate the condition number of the unpreconditioned Galerkin matrix to motivate theneed for good preconditioning. We will work on the reference triangle (cid:98) K = conv { (0 , , (1 , , (0 , } and will needthe following well-known inverse inequalities for polynomials on (cid:98) K : Proposition 4.1 (Inverse inequalities, [37, Theorem 4.76]) . Let (cid:98) K denote the reference triangle and let E be oneof its edges. There exists a constant C such that for all p ∈ N and for all v ∈ P p ( (cid:98) K ) the following estimates hold: (cid:107) v (cid:107) L ∞ ( (cid:98) K ) ≤ Cp (cid:107) v (cid:107) L ( (cid:98) K ) , (4.1) (cid:107) v (cid:107) L ∞ ( (cid:98) K ) ≤ C (cid:112) log( p + 1) (cid:107) v (cid:107) H ( (cid:98) K ) , (4.2) (cid:107) v (cid:107) H ( (cid:98) K ) ≤ Cp (cid:107) v (cid:107) L ( (cid:98) K ) , (4.3) (cid:107) v (cid:107) L ( E ) ≤ Cp (cid:107) v (cid:107) L ( (cid:98) K ) . (4.4)First we investigate the L and H conditioning of our basis on the reference triangle. Lemma 4.2.
Let u ∈ P p ( (cid:98) K ) and let α V j , α E m j , α I ( ij ) be the coefficients with respect to the basis in Definition 2.4,i.e., we decompose u = u V + u E + u E + u E + u I with u V = (cid:88) j =1 α V j ϕ V j , u E m = p − (cid:88) j =0 α E m j ϕ E m j , u I = (cid:88) i + j ≤ p − α I ( ij ) ϕ I ( ij ) . (4.5) Then for a constant
C > that does not depend on u or p : (cid:107) u V (cid:107) L ( (cid:98) K ) ≤ C (cid:88) j =1 (cid:12)(cid:12) α V j (cid:12)(cid:12) , (cid:107) u E m (cid:107) L ( (cid:98) K ) ≤ C p − (cid:88) j =0 (cid:12)(cid:12)(cid:12) α E m j (cid:12)(cid:12)(cid:12) , (cid:107) u I (cid:107) L ( (cid:98) K ) = (cid:88) i + j ≤ p − (cid:12)(cid:12)(cid:12) α I ( ij ) (cid:12)(cid:12)(cid:12) . (4.6) Combined this gives: (cid:107) u (cid:107) L ( (cid:98) K ) ≤ C (cid:16) (cid:88) j =1 (cid:12)(cid:12) α V j (cid:12)(cid:12) + (cid:88) m =1 p − (cid:88) j =0 (cid:12)(cid:12)(cid:12) α E m j (cid:12)(cid:12)(cid:12) + (cid:88) i + j ≤ p − (cid:12)(cid:12)(cid:12) α I ( ij ) (cid:12)(cid:12)(cid:12) (cid:17) . (4.7)9 roof. The estimate for u V in (4.6) is clear. For the edge contributions in (4.6), we restrict ourselves to the edge(0 , − (1 , m = 1 and drop the index m in the notation. The other edges can be treated analogously. (cid:107) u E (cid:107) L ( (cid:98) K ) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) p − (cid:88) j =0 α E j (cid:114) i + 32 L S i +2 ( λ − λ , λ + λ ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( (cid:98) K ) = p − (cid:88) i,j =0 α E j α E i (cid:90) (cid:98) K (cid:114) i + 32 L S i +2 ( λ − λ , λ + λ ) (cid:114) j + 32 L S j +2 ( λ − λ , λ + λ ) dx = 14 p − (cid:88) i,j =0 α E j α E i (cid:90) − (cid:90) − (cid:114) i + 32 L i +2 ( ξ ) (cid:114) j + 32 L j +2 ( ξ ) (cid:18) − η (cid:19) i + j +5 dξdη. In the last step we transformed the reference triangle to ( − , × ( − ,
1) via the map ( ξ, η ) (cid:55)→ ( (1 + ξ )(1 − η ) , (1 + η )).It is well-known (see [37, p. 65]) that the 1D-mass matrix M of the integrated Legendre polynomials is pen-tadiagonal, and the non-zero entries satisfy (cid:12)(cid:12) M ( ij ) (cid:12)(cid:12) ∼ i +1) ( j +1) . It is easy to check that (cid:12)(cid:12)(cid:12)(cid:82) − (cid:0) − η (cid:1) i + j +5 dη (cid:12)(cid:12)(cid:12) ≤ C ( i + j + 6) − . Together with a Cauchy-Schwarz estimate, we obtain: (cid:107) u E (cid:107) L ( (cid:98) K ) (cid:46) p − (cid:88) j =0 (cid:12)(cid:12) α E j (cid:12)(cid:12) j + 1) (cid:46) p − (cid:88) j =0 (cid:12)(cid:12) α E j (cid:12)(cid:12) . The bubble basis functions are chosen L -orthogonal. Thus, using our scaling of the bubble basis functions, (cid:107) u I (cid:107) L ( (cid:98) K ) = (cid:88) i + j ≤ p − (cid:12)(cid:12)(cid:12) α I ( ij ) (cid:12)(cid:12)(cid:12) (cid:13)(cid:13) ϕ ( ij ) (cid:13)(cid:13) L ( (cid:98) K ) = (cid:88) i + j ≤ p − (cid:12)(cid:12)(cid:12) α I ( ij ) (cid:12)(cid:12)(cid:12) . Finally, we split the function u into vertex, edge and inner components, apply the triangle inequality and get (cid:107) u (cid:107) L ( (cid:98) K ) ≤ (cid:107) u V (cid:107) L ( (cid:98) K ) + 5 (cid:88) m =1 (cid:13)(cid:13) u E j (cid:13)(cid:13) L ( (cid:98) K ) + 5 (cid:107) u I (cid:107) L ( (cid:98) K ) , which is (4.7).More interesting are the reverse estimates. Lemma 4.3.
There is a constant
C > independent of p such for every u ∈ P p ( (cid:98) K ) the coefficients of itsrepresentation in the basis of Definition 2.4 as in Lemma 4.2 satisfy:vertex parts: (cid:88) j =1 (cid:12)(cid:12) α V j (cid:12)(cid:12) ≤ Cp (cid:107) u (cid:107) L ( (cid:98) K ) as well as (cid:88) j =1 (cid:12)(cid:12) α V j (cid:12)(cid:12) ≤ C log( p + 1) (cid:107) u (cid:107) H ( (cid:98) K ) , edge parts: p − (cid:88) j =0 (cid:12)(cid:12)(cid:12) α E m j (cid:12)(cid:12)(cid:12) ≤ Cp (cid:107) u (cid:107) L ( (cid:98) K ) and p − (cid:88) j =0 (cid:12)(cid:12)(cid:12) α E m j (cid:12)(cid:12)(cid:12) ≤ Cp (cid:107) u (cid:107) H ( (cid:98) K ) . Moreover, if u vanishes on ∂ (cid:98) K , then (cid:88) i + j ≤ p − (cid:12)(cid:12)(cid:12) α I ( ij ) (cid:12)(cid:12)(cid:12) = (cid:107) u (cid:107) L ( (cid:98) K ) . (4.8)10 roof. Since α V j = u ( z j ) where z j denotes the j − th vertex, we can use the L ∞ -inverse estimates (4.1) and (4.2)to get estimates for the vertex part.For the edge parts, we again only consider the bottom edge, E = E m with m = 1. First we assume that u vanishes in all vertices. If we consider the restriction of u to the edge E we only have contributions by the edgebasis, i.e., we can write u ( x,
0) = p − (cid:88) i =0 α E i (cid:114) i + 32 L i +2 (2 x − , x ∈ (0 , ,∂∂x u ( x,
0) = p − (cid:88) i =0 α E i (cid:114) i + 32 (cid:96) i +1 (2 x − , x ∈ (0 , . The factor was chosen to get an L -normalized basis, since we have (cid:107) (cid:96) i +1 (cid:107) L ( − , = i +3 . The Legendrepolynomials are orthogonal on ( − , (cid:13)(cid:13)(cid:13)(cid:13) ∂u∂x (cid:13)(cid:13)(cid:13)(cid:13) L ( E ) = 2 p − (cid:88) i =0 (cid:12)(cid:12) α E i (cid:12)(cid:12) . If we consider a general u ∈ P p ( (cid:98) K ), we apply the previous estimate to u := u − I u where I denotes thenodal interpolation operator to the linears. Then we get from the triangle inequality p − (cid:88) i =0 (cid:12)(cid:12) α E i (cid:12)(cid:12) ≤ (cid:13)(cid:13)(cid:13)(cid:13) ∂u∂x (cid:13)(cid:13)(cid:13)(cid:13) L ( E ) + 2 (cid:13)(cid:13)(cid:13)(cid:13) ∂I u∂x (cid:13)(cid:13)(cid:13)(cid:13) L ( E ) . We apply the trace estimate (4.4) to the first part and the trace and norm equivalence for the second. We obtain: p − (cid:88) i =0 (cid:12)(cid:12) α E i (cid:12)(cid:12) (cid:46) p (cid:13)(cid:13)(cid:13)(cid:13) ∂u∂x (cid:13)(cid:13)(cid:13)(cid:13) L ( (cid:98) K ) + (cid:13)(cid:13) I u (cid:13)(cid:13) L ( (cid:98) K ) . The H estimate then follows from the L ∞ estimate for the nodal interpolant (4.2). For the L estimate we thensimply use the inverse estimate (4.3).For the equality (4.8) we note that if u | ∂T = 0 then u = u I and thus we can use the equality in (4.6). Lemma 4.4.
There exist constants c , C , c , C > independent of p such that for every u ∈ P p ( (cid:98) K ) itscoefficients in the basis of Definition 2.4 as in Lemma 4.2 satisfy: c (cid:107) u (cid:107) L ( (cid:98) K ) ≤ (cid:88) j =1 (cid:12)(cid:12) α V j (cid:12)(cid:12) + (cid:88) m =1 p − (cid:88) j =0 (cid:12)(cid:12)(cid:12) α E m j (cid:12)(cid:12)(cid:12) + (cid:88) i + j ≤ p − (cid:12)(cid:12)(cid:12) α I ( ij ) (cid:12)(cid:12)(cid:12) ≤ C p (cid:107) u (cid:107) L ( (cid:98) K ) , (4.9) c p − (cid:107) u (cid:107) H ( (cid:98) K ) ≤ (cid:88) j =1 (cid:12)(cid:12) α V j (cid:12)(cid:12) + (cid:88) m =1 p − (cid:88) j =0 (cid:12)(cid:12)(cid:12) α E m j (cid:12)(cid:12)(cid:12) + (cid:88) i + j ≤ p − (cid:12)(cid:12)(cid:12) α I ( ij ) (cid:12)(cid:12)(cid:12) ≤ C p (cid:107) u (cid:107) H ( (cid:98) K ) . (4.10) Proof of (4.9) : The lower bound was already shown in Lemma 4.2. For the upper bound we apply the precedingLemma 4.3 to u and see (cid:80) j =1 (cid:12)(cid:12) α V j (cid:12)(cid:12) (cid:46) p (cid:107) u (cid:107) L ( (cid:98) K ) . For any edge E m we get p − (cid:88) j =0 (cid:12)(cid:12)(cid:12) α E m j (cid:12)(cid:12)(cid:12) (cid:46) p (cid:107) u (cid:107) L ( (cid:98) K ) . − − − − O ( p − ) O ( p − ) O ( p − )polynomial degree λ m i n M E M E + S E MM + S − O ( p ) O (1)polynomial degree λ m a x M I M I + S I MM + S Figure 4.1: Numerical computation of the extremal eigenvalues of the mass matrices and the sum of mass and H -stiffness matrix M + S for the full system and different sub-blocks. Next we set u := u − u V − u E , where u E is the sum of the edge contributions u E m . This function vanishes on theboundary of (cid:98) K and we can apply Lemma 4.3 to get: (cid:88) i + j ≤ p − (cid:12)(cid:12)(cid:12) α I ( ij ) (cid:12)(cid:12)(cid:12) (cid:46) (cid:107) u (cid:107) L ( (cid:98) K ) . Since we can always estimate the L norms by the (cid:96) norms of the coefficients (Lemma 4.2), we obtain (cid:107) u (cid:107) L ( (cid:98) K ) ≤ (cid:107) u (cid:107) L ( (cid:98) K ) + 3 (cid:107) u V (cid:107) L ( (cid:98) K ) + 3 (cid:107) u E (cid:107) L ( (cid:98) K ) (cid:46) (cid:107) u (cid:107) L ( (cid:98) K ) + p (cid:107) u (cid:107) L ( (cid:98) K ) + p (cid:107) u (cid:107) L ( (cid:98) K ) . Proof of (4.10) : The proof for the upper H estimate works along the same lines, but using the sharper H -estimates from Lemma 4.3 for the vertex and edge parts. For the lower estimate, we just make use of the inverseestimate (4.3) and the fact that the L -norm is uniformly bounded by the coefficients, to conclude the proof. Example 4.5.
In Figure 4.1 we compare our theoretical bounds on the reference element from Lemma 4.4 with anumerical experiment that studies the maximal and minimal eigenvalues of the mass matrix M and the stiffnessmatrix S (corresponding to the bilinear form ( ∇· , ∇· ) L ( (cid:98) K ) ). We focus on the full system and the subblocks thatcontributed the highest order in our theoretical investigations, i.e. the edge blocks M E , S E , and the block of innerbasis functions M I , S I . We see that the estimates on the full condition numbers are not overly pessimistic: thenumerics show a behavior of the minimal eigenvalue of O ( p . ) instead of O ( p ) . If we focus solely on the edgecontributions, we see that the bound we used for the lower eigenvalue is not sharp there. This can partly be explainedby the fact that if no inner basis functions are present it is possible to improve the estimate (4.4) by a factor of p .But since we also need to include the coupling of inner and edge basis functions this improvement in order is lostagain when looking at the full systems. The estimates on the reference triangle can now be transferred to the global space (cid:101) S p ( T ) on quasiuniformmeshes. Theorem 4.6.
Let T be a quasiuniform triangulation with mesh size h . With the polynomial basis on the referencetriangle (cid:98) K given by Definition 2.4, let { ϕ i | i = 1 , . . . , N } be the basis of (cid:101) S p ( T ) . Then there exist constants c , c / , c , C , C / , C > that depend only on Γ and the γ -shape regularity of T , such that, for every u ∈ R N and = (cid:80) Nj =1 u j ϕ j ∈ (cid:101) S p ( T ) : c h (cid:107) u (cid:107) L (Γ) ≤ (cid:107) u (cid:107) (cid:96) ≤ C p h (cid:107) u (cid:107) L (Γ) , (4.11) c p − (cid:107) u (cid:107) H (Γ) ≤ (cid:107) u (cid:107) (cid:96) ≤ C (cid:0) p + h − (cid:1) (cid:107) u (cid:107) H (Γ) , (4.12) c / h − p − (cid:107) u (cid:107) (cid:101) H / (Γ) ≤ (cid:107) u (cid:107) (cid:96) ≤ C / (cid:18) p h + h − (cid:19) (cid:107) u (cid:107) (cid:101) H / (Γ) . (4.13) Proof.
The L -estimate (4.11) can easily be shown by transforming to the reference element and applying (4.9).To prove the other estimates (4.12), (4.13), we need the Scott-Zhang projection operator J h : L (Γ) → (cid:101) S ( T )as modified in [3, Section 3.2]. It has the following important properties:1. J h is a bounded linear operator from L (Γ) to ( (cid:101) S ( T ) , (cid:107) · (cid:107) L (Γ) ).2. For every s ∈ [0 ,
1] there holds (cid:107) J h v (cid:107) (cid:101) H s (Γ) ≤ C stab ( s ) (cid:107) v (cid:107) (cid:101) H s (Γ) ∀ v ∈ (cid:101) H s (Γ).3. For every K ∈ T let ω K := (cid:83) { K (cid:48) ∈ T : K ∩ K (cid:48) (cid:54) = ∅} denote the element patch, i.e., the union of all elementsthat touch K . Then, for all v ∈ (cid:101) H (Γ) (cid:107) (1 − J h ) v (cid:107) L ( K ) ≤ C sz h K (cid:107)∇ v (cid:107) L ( ω K ) , (4.14) (cid:107)∇ (1 − J h ) v (cid:107) L ( K ) ≤ C sz (cid:107)∇ v (cid:107) L ( ω K ) . (4.15)The constant C sz depends only on the γ -shape regularity of T , and C stab ( s ) additionally depends on Γ and s .We will use the following notation: For a function u ∈ (cid:101) S p ( T ) we will write u ∈ R N for its representation in thebasis { ϕ i | i = 1 , . . . , N } . For an element K ∈ T we write u | K for the part of the coefficient vector that belongs tobasis functions whose support intersects the interior of K . In addition to the function u ∈ (cid:101) S p ( T ), we will employthe function ˜ u := u − J h u . Its vector representation will be denoted ˜ u ∈ R N . Finally, the vector representation of J h u (again u ∈ (cid:101) S p ( T )) will be J h u ∈ R N .
1. step:
We claim the following stability estimates: (cid:107) J h u (cid:107) (cid:96) (cid:46) h − (cid:107) J h u (cid:107) L (Γ) (cid:46) h − (cid:107) u (cid:107) L (Γ) (cid:46) (cid:107) u (cid:107) (cid:96) , (4.16) (cid:107) J h u (cid:107) (cid:101) H / (Γ) (cid:46) h − (cid:107) J h u (cid:107) L (Γ) , (4.17) (cid:107) J h u (cid:107) (cid:101) H / (Γ) (cid:46) h (cid:107) u (cid:107) (cid:96) . (4.18)The inequalities (4.16) are just a simple scaling argument combined with the L stability of the Scott-Zhangprojection and (4.11). The inequality (4.17) follows from the inverse inequality (note that J h u has degree 1).Finally, (4.18) follows from combining (4.17) and (4.16).
2. step:
Next, we investigate the function ˜ u = u − J h u . We claim the following estimates: (cid:107) ˜ u (cid:107) L (Γ) (cid:46) h (cid:107) u (cid:107) (cid:96) , (4.19) (cid:107) ˜ u (cid:107) H (Γ) (cid:46) p (cid:107) u (cid:107) (cid:96) , (4.20) (cid:107) ˜ u (cid:107) (cid:101) H / (Γ) (cid:46) p h (cid:107) u (cid:107) (cid:96) . (4.21)The estimate (4.19) is a simple consequence of (4.11) and the L -stability of the Scott-Zhang operator J h . For theproof of (4.20), we combine a simple scaling argument with (4.10) and the stability estimate (4.16) to get (cid:107) ˜ u (cid:107) H (Γ) (4.10) (cid:46) p (cid:107) ˜ u (cid:107) (cid:96) (4.16) (cid:46) p (cid:107) u (cid:107) (cid:96) .
3. step:
We assert: (cid:107) ˜ u (cid:107) (cid:96) (cid:46) p h (cid:107) u (cid:107) L (Γ) , (4.22) (cid:107) ˜ u (cid:107) (cid:96) (cid:46) p (cid:107) u (cid:107) H (Γ) , (4.23) (cid:107) ˜ u (cid:107) (cid:96) (cid:46) p h (cid:107) u (cid:107) (cid:101) H / (Γ) . (4.24)Again, (4.22) is a simple consequence of (4.11) and the L -stability of the Scott-Zhang operator J h . For the bound(4.23) we calculate, using the equivalence (4.10) of the coefficient vector and the H -norm on the reference triangle,together with the scaling properties of the H - and L -norms, (cid:107) ˜ u (cid:107) (cid:96) ≤ (cid:88) K ∈T (cid:107) ˜ u | K (cid:107) (cid:96) (4.10) , scaling (cid:46) p (cid:88) K ∈T (cid:16) h − (cid:107) ˜ u (cid:107) L ( K ) + | ˜ u | H ( K ) (cid:17) . By applying the local L -interpolation estimate (4.14) and H -stability (4.15) we get (cid:107) ˜ u (cid:107) (cid:96) (cid:46) p (cid:88) K ∈T (cid:107) u (cid:107) H ( ω K ) (cid:46) p (cid:107) u (cid:107) H (Γ) , where in the last step we used the fact for shape regular meshes each element is contained in at most M differentpatches, where M depends solely on the shape regularity constant γ .We next prove (4.24). We apply Proposition 2.1, (i) to the map (cid:96) → (cid:101) S p ( T ) : u (cid:55)→ u , where the space (cid:101) S p ( T ) isonce equipped with the L - and once with the H -norm. By Proposition 3.4 interpolating between (4.19) and (4.20)yields (cid:107) ˜ u (cid:107) (cid:96) (cid:46) (cid:0) p h − p (cid:1) / (cid:107) u (cid:107) (cid:101) H / (Γ) (cid:46) p h − (cid:107) u (cid:107) (cid:101) H / (Γ) .
4. step:
The above steps allow us to obtain the H and H / estimates (4.12)—(4.13) of the theorem. Wedecompose u = ˜ u + J h u and correspondingly u = ˜ u + J h u . Then: (cid:107) u (cid:107) (cid:96) (cid:46) (cid:107) ˜ u (cid:107) (cid:96) + (cid:107) J h u (cid:107) (cid:96) , (4 . , (4 . (cid:46) p h (cid:107) u (cid:107) (cid:101) H / (Γ) + 1 h (cid:107) u (cid:107) L (Γ) (cid:46) p h + 1 h (cid:107) u (cid:107) (cid:101) H / (Γ) , (cid:107) u (cid:107) (cid:101) H / (Γ) (cid:46) (cid:107) ˜ u (cid:107) (cid:101) H / (Γ) + (cid:107) J h u (cid:107) (cid:101) H / (Γ) (4 . , (4 . (cid:46) p h (cid:107) u (cid:107) (cid:96) + h (cid:107) u (cid:107) (cid:96) (cid:46) h p (cid:107) u (cid:107) (cid:96) . This shows (4.13). The H estimate (4.12) follows along the same lines: An elementwise inverse estimate gives (cid:107) u (cid:107) H (Γ) (cid:46) p h (cid:107) u (cid:107) L (Γ) (4.11) (cid:46) p (cid:107) u (cid:107) (cid:96) , and the splitting u = ˜ u + J h u produces (cid:107) u (cid:107) (cid:96) (cid:46) (cid:107) ˜ u (cid:107) (cid:96) + (cid:107) J h u (cid:107) (cid:96) (4 . , (4 . (cid:46) (cid:0) p + h − (cid:1) (cid:107) u (cid:107) H (Γ) . Corollary 4.7.
The spectral condition number of the unpreconditioned Galerkin matrix (cid:101) D ph can be bounded by κ (cid:16) (cid:101) D ph (cid:17) ≤ C (cid:18) p h + p (cid:19) with a constant C > that depends only on Γ and the γ -shape regularity of T . roof. The bilinear form induced by the stabilized hypersingular operator is elliptic and continuous with respectto the (cid:101) H / -norm. By applying the estimates (4.13) to the Rayleigh quotients we get the stated result. Remark 4.8.
In this section we did not consider the effect of diagonal scaling. The numerical results in Section 6suggest that it improves the p -dependence of the condition number significantly.4.1. Quadrilateral meshes The present paper focuses on meshes consisting of triangles. Nevertheless, in order to put the results of Section 4in perspective, we include a short section on quadrilateral meshes. In [29], estimates similar to those of Lemma 4.4have been derived for the case of the Babuˇska-Szab´o basis on the reference square (cid:98) S := [ − , . We give a briefsummary of the definitions and results. Definition 4.9 (Babuˇska-Szab´o basis) . (i) On the reference interval (cid:98) I := [ − , , the basis functions are basedon the integrated Legendre polynomials L j as defined in (2.10) : ϕ ( x ) := 12 (1 − x ) , ϕ ( x ) := 12 (1 + x ) , ϕ j ( x ) := 1 (cid:107) (cid:96) j − (cid:107) L ( I ) L j ( x ) ∀ ≤ j ≤ p. (ii) On the reference square (cid:98) S , the basis of the “tensor product space” Q p ( (cid:98) S ) is given by the tensor product of the1D basis functions: { ϕ i ⊗ ϕ j : 0 ≤ i, j ≤ p } . For this basis, the following estimates hold:
Proposition 4.10 ([29, Theorem 1], [22, Theorem 4.1]) . Let u ∈ Q p ( (cid:98) S ) and let u denote its coefficient vector withrespect to the basis of Definition 4.9. Then the following estimates hold: (cid:107) u (cid:107) L ( (cid:98) S ) (cid:46) (cid:107) u (cid:107) (cid:96) (cid:46) p (cid:107) u (cid:107) L ( (cid:98) S ) , (cid:107) u (cid:107) H ( (cid:98) S ) (cid:46) (cid:107) u (cid:107) (cid:96) (cid:46) p (cid:107) u (cid:107) H ( (cid:98) S ) . (4.25) Remark 4.11.
In [29, Thm. 1], the estimates were only shown for the inner degrees of freedom, i.e., if u | ∂ (cid:98) S = 0 .This restriction is removed in [22, Thm. 4.1]. Theorem 4.12.
Let T be a quasi-uniform, shape-regular affine mesh of quadrilaterals of size h , and let F K : (cid:98) S → K be the affine element map for K ∈ T . Let u ∈ (cid:101) Q p ( T ) := (cid:110) u ∈ (cid:101) H / (Γ) : u | K ◦ F K ∈ Q p ( (cid:98) S ) ∀ K ∈ T (cid:111) , and let u denote its coefficient vector with respect to the basis of Definition 4.9. Then there exist constants c , c , C , C > that depend only on Γ and the γ -shape regularity of T such that: c h − (cid:107) u (cid:107) L (Γ) ≤ (cid:107) u (cid:107) (cid:96) ≤ C h − p (cid:107) u (cid:107) L (Γ) , c (cid:107) u (cid:107) H (Γ) ≤ (cid:107) u (cid:107) (cid:96) ≤ C (cid:0) p + h − (cid:1) (cid:107) u (cid:107) H (Γ) . Proof.
The proof is completely analogous to that of Theorem 4.6. The only additional ingredient to Proposition 4.10is an operator J h : L → (cid:101) Q ( T ) that is bounded with respect to the L and the H -norm, reproduces homogeneousDirichlet boundary conditions for the case of open surfaces, and has the approximation property (cid:107) u − J h u (cid:107) L (Γ) (cid:46) h | u | H (Γ) . Such an operator was proposed, e.g., in [6]. The important estimates that need to be shown are (weagain write (cid:101) u := u − J h u ): (cid:107) u (cid:107) L (Γ) (cid:46) h (cid:107) u (cid:107) (cid:96) , (cid:107) u (cid:107) H (Γ) (cid:46) (cid:107) u (cid:107) (cid:96) (cid:107) u (cid:107) (cid:96) (cid:46) p h (cid:107) u (cid:107) L (Γ) , (cid:107) u (cid:107) (cid:96) (cid:46) (cid:107) (cid:101) u (cid:107) (cid:96) + (cid:107) J h u (cid:107) (cid:96) (cid:46) p (cid:107) u (cid:107) H (Γ) + h − (cid:107) u (cid:107) L (Γ) . emark 4.13. In the case of triangular meshes, Proposition 3.4 allowed us to infer (cid:101) H / -condition numberestimates in Corollary 4.7 by interpolating the discrete norm equivalences in L and H of Theorem 4.6. For meshesconsisting of quadrilaterals, the result corresponding to Proposition 3.4 is currently not available in the literature.If we conjecture (cid:104)(cid:16) (cid:101) Q p ( T ) , (cid:107)·(cid:107) L (Γ) (cid:17) , (cid:16) (cid:101) Q p ( T ) , (cid:107)·(cid:107) H (Γ) (cid:17)(cid:105) / = (cid:16) (cid:101) Q p ( T ) , (cid:107)·(cid:107) (cid:101) H / (Γ) (cid:17) with equivalent norms, thenTheorem 4.12 implies the following estimates: h − (cid:107) u (cid:107) (cid:101) H / (Γ) (cid:46) (cid:107) u (cid:107) (cid:96) (cid:46) (cid:0) h − p + h − (cid:1) (cid:107) u (cid:107) (cid:101) H / (Γ) ,κ (cid:16) (cid:101) D ph (cid:17) (cid:46) h − + p . Remark 4.14. [29] also analyzes the influence of diagonal preconditioning and shows that the condition numberis improved by a factor of two in the exponents of p . Although we did not make any theoretical investigations inthis direction for the H / -case, our numerical experiments in Examples 6.2 and 6.3 for triangular meshes showthat diagonal scaling improves the p -dependence of the condition number from O ( p . ) to O ( p . ) . Remark 4.15.
The Babuˇska-Szab´o basis of Definition 4.9 is not the only one used on quadrilaterals or hexahedra.An important representative of other bases are the Lagrange interpolation polynomials associated with the Gauß-Lobatto points. This basis has the better O ( p ) conditioning for the stiffness matrix and O ( p ) for the mass matrix(see [32],[29, Sect. 6]). Using the same arguments as in the proof of Theorem 4.12, we get for the global H -problemthat the condition number behaves like O ( p h − + p ) . If the conjecture of Remark 4.13 is valid, then we obtain forthis basis for the hypersingular integral operator the condition number estimate κ ( (cid:101) D ph ) (cid:46) p / + p − / h − . (See theAppendix for details.) hp -preconditioning Additive Schwarz preconditioners are based on decomposing a vector space V into smaller subspaces V i , i =0 , . . . , J , on which a local problem is solved. We recall some of the basic definitions and important results. Detailscan be found in [42, chapter 2].Let a ( · , · ) : V × V → R be a symmetric, positive definite bilinear form on the finite dimensional vector space V . For a given f ∈ V (cid:48) consider the problem of finding u ∈ V such that a ( u, v ) = f ( v ) ∀ v ∈ V . We will write A for the corresponding Galerkin matrix.Let V i ⊂ V , i = 0 , . . . , J , be finite dimensional vector spaces with corresponding prolongation operators R Ti : V i → V . We will commit a slight abuse of notation and also denote the matrix representation of the operatorby R Ti , and R i is its transposed matrix. We also assume that V permits a (in general not direct) decompositioninto V = R T V + J (cid:88) i =1 R Ti V i . We assume that for each subspace V i a symmetric and positive definite bilinear form (cid:101) a i ( · , · ) : V i × V i → R , i = 0 , . . . , J, is given. We write (cid:101) A i for the matrix representation of (cid:101) a i ( · , · ). Sometimes these bilinear forms are referred to asthe “local solvers”; in the simplest case of “exact local solvers” they are just restrictions of a ( · , · ), i.e., (cid:101) a i ( u i , v i ) :=16 (cid:0) R Ti u i , R Ti v i (cid:1) for all u i , v i ∈ V i . Then, the corresponding additive Schwarz preconditioner is given by B − := J (cid:88) i =0 R Ti (cid:101) A − i R i . The following proposition allows us to bound the condition number of the preconditioned system B − A . Thefirst part is often referred to as the Lemma of Lions (see [51, 27, 30]). Proposition 5.1. (a) Assume that there exists a constant C > such that every u ∈ V admits a decomposition u = (cid:80) Ji =0 R Ti u i with u i ∈ V i such that J (cid:88) i =0 (cid:101) a i ( u i , u i ) ≤ C a ( u, u ) . Then, the minimal eigenvalue λ min ( B − A ) of B − A satisfies λ min (cid:0) B − A (cid:1) ≥ C − . (b) Assume that there exists C > such that for every decomposition u = (cid:80) Ji =0 R Ti v i with v i ∈ V i the followingestimate holds: a ( u, u ) ≤ C J (cid:88) i =0 (cid:101) a i ( v i , v i ) . Then, the maximal eigenvalue λ max ( B − A ) of B − A satisfies λ max (cid:0) B − A (cid:1) ≤ C .(c) These two estimates together give an estimate for the condition number of the preconditioned linear system: κ (cid:0) B − A (cid:1) := λ max λ min ≤ C C . hp -stable preconditioner In order to define an additive Schwarz preconditioner, we decompose the boundary element space V := (cid:101) S p ( T )into several overlapping subspaces. We define V h := (cid:101) S ( T ) as the space of globally continuous and piecewise linearfunctions on T that vanish on ∂ Γ and denote the corresponding canonical embedding operator by R Th : V h → V .We also define for each vertex z ∈ V the local space V p z := { u ∈ (cid:101) S p ( T ) | supp( u ) ⊂ ω z } and denote the canonical embedding operators by R T z : V p z → V . The space decomposition then reads V = V h + (cid:88) z ∈V V p z . (5.1)We will denote the restriction of the Galerkin matrix (cid:101) D ph to the subspaces V h and V p z as (cid:101) D h and (cid:101) D ph, z , respectively. Lemma 5.2.
There exist constants c , c > , which depend only on Γ and the γ -shape regularity of T , such thatthe following holds:(a) For every u ∈ (cid:101) S p ( T ) there exists a decomposition u = u + (cid:80) z ∈V u z with u ∈ V h and u z ∈ V p z and (cid:107) u (cid:107) (cid:101) H / (Γ) + (cid:88) z ∈V (cid:107) u z (cid:107) (cid:101) H / (Γ) ≤ c (cid:107) u (cid:107) (cid:101) H / (Γ) . b) Any decomposition u = v + (cid:80) z ∈V v z with v ∈ V h and v z ∈ V p z satisfies (cid:107) u (cid:107) (cid:101) H / (Γ) ≤ c (cid:32) (cid:107) v (cid:107) (cid:101) H / (Γ) + (cid:88) z ∈V (cid:107) v z (cid:107) (cid:101) H / (Γ) (cid:33) . Proof.
The first estimate is the assertion of Proposition 3.1, (iv). The second estimate can be shown by a so-calledcoloring argument, along the same lines as in [18, Lemma 2]. It is based on the following estimate (see [35, Lemma4.1.49] or [46, Lemma 3.2]): Let w j , j = 1 , . . . , n be functions in (cid:101) H s (Γ) for s ≥ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) n (cid:88) i =1 w i (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H s (Γ) ≤ C n (cid:88) i =1 (cid:107) w i (cid:107) (cid:101) H s (Γ) , (5.2)where C > γ -shape regularity, the number of elements in any vertex patch, and thereforealso the number of vertices in a patch, is uniformly bounded by some constant N c which depends solely on γ . Thus,we can divide the vertices into sets J , . . . , J N c such that (cid:83) N c i =1 J i = V and | ω z ∩ ω z (cid:48) | = 0 for all z , z (cid:48) in the sameindex set J i . Repeated application of the triangle inequality together with (5.2) then gives: (cid:107) u (cid:107) (cid:101) H / (Γ) ≤ (cid:107) v (cid:107) (cid:101) H / (Γ) + 2 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:88) z ∈V v z (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H / (Γ) ≤ (cid:107) v (cid:107) (cid:101) H / (Γ) + 2 N c N c (cid:88) i =1 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:88) z ∈ J i v z (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H / (Γ) ≤ (cid:107) v (cid:107) (cid:101) H / (Γ) + 2 N c C (cid:88) z ∈V (cid:107) v z (cid:107) (cid:101) H / (Γ) . The previous lemma only made statements about the (cid:101) H / (Γ)-norm. Theorem 5.3.
Let T be a γ -shape regular triangulation of Γ . Then there is a constant C > that depends solelyon Γ and the γ -shape regularity of T such that the following is true: The preconditioner B − := R Th (cid:16) (cid:101) D h (cid:17) − R h + (cid:88) z ∈V R T z (cid:16) (cid:101) D ph, z (cid:17) − R z , which is implied by the space decomposition (5.1) , leads to the spectral condition number estimate κ ( B − (cid:101) D ph ) ≤ C. Proof.
The bilinear form (cid:68) (cid:101) D · , · (cid:69) Γ is equivalent to (cid:107)·(cid:107) (cid:101) H / (Γ) . Hence, the combination of Lemma 5.2 and Proposi-tion 5.1 give the boundedness of the condition number. The preconditioner of Theorem 5.3 relies on the space decomposition (5.1). In this section, we discuss howthe space (cid:101) S ( T ) of piecewise linear function can be further decomposed in a multilevel fashion. Our setting willbe one where T is the finest mesh of a sequence ( T (cid:96) ) L(cid:96) =0 of nested meshes that are generated by newest vertexbisection (NVB); see Figure 5.1 for a description. We point the reader to [40, 24] for a detailed discussion of NVB.A key feature of NVB is that it creates sequences of meshes that are uniformly shape regular. We mention inpassing that further properties of NVB were instrumental in proving optimality of h -adaptive algorithms in bothFEM [8, 13] and BEM [11, 12, 15, 17]. Before discussing the details of the multilevel space decomposition, westress that the preconditioner described in Section 5.2 is independent of the chosen refinement strategy (such as18 igure 5.1: For each element K ∈ T (cid:96) there exists exactly one reference edge indicated by the red line (upper left plot). The element K is refined by bisecting its reference edge. This leads to a new node (red dot) and two son elements. The reference edges of the sonelements are opposite to the newest vertex (lower left plot). Hanging nodes are avoided as follows: Assume that some of the edges ofthe triangle, but at least the reference edge, are marked for refinement (upper plots). The triangle will be split into two, three or fourson elements by iterative application of the newest vertex bisection (NVB) . NVB) as long as it satisfies the assumptions in Section 5.2, whereas the condition number estimates for the localmultilevel preconditioner discussed in the present Section 5.3 depend on the fact that the underlying refinementstrategy is based on NVB.Adaptive algorithms create sequences of meshes ( T (cid:96) ) L(cid:96) =0 . Typically, the procedure starts with an initial tri-angulation T and the further members of the sequences are created inductively. That is, mesh T (cid:96) is obtainedfrom T (cid:96) − by refining some elements of T (cid:96) − . In an adaptive environment, these elements are determined by amarking criterion (“marked elements”) and a mesh closure condition. Usually, the following assumptions on themesh refinement are made: • T (cid:96) is regular for all (cid:96) ∈ N , i.e., there exist no hanging nodes; • The meshes T , T , . . . are uniformly γ -shape-regular, i.e., with | K | denoting the surface area of an element K ∈ T (cid:96) and diam( K ) the Euclidean diameter, we havesup (cid:96) ∈ N max K ∈T (cid:96) diam( K ) | K | ≤ γ. (5.3)We consider a sequence of triangulations T , . . . , T L , which is created by iteratively applying NVB. The corre-sponding sets of vertices are denoted V , . . . , V L . For a vertex z ∈ V (cid:96) , the associated patch is denoted by ω (cid:96), z .In the construction of the p -preconditioner in Section 5.2 we only considered a single mesh T . For the remainderof the paper, the p part will always be constructed with respect to the finest mesh T L . For a simpler presentationwe set T := T L and V := V L . The space decomposition from (5.1) involves the global lowest-order space V h = (cid:101) S ( T L ). Therefore, thecomputation of the corresponding additive Schwarz operator needs the inversion of a global problem, which is, inpractice, very costly, and often even infeasible. To overcome this disadvantage, we consider a refined splitting ofthe space V h that relies on the hierarchy of the adaptively refined meshes T , . . . , T L . The corresponding localmultilevel preconditioner was introduced and analyzed in [14, 16]. See also [20, 47, 49, 48] for local multilevelpreconditioners for (adaptive) FEM , [43, 19] for (uniform) BEM, and [2, 28, 44] for (restricted) approaches foradaptive BEM.Set (cid:101) V := V and define the local subsets (cid:101) V (cid:96) := V (cid:96) \V (cid:96) − ∪ { z ∈ V (cid:96) − : ω (cid:96), z (cid:40) ω (cid:96) − , z } for (cid:96) ≥ igure 5.2: Visualization of the definition (5.4) of the local subsets (cid:101) V (cid:96) : Starting with a triangulation T (cid:96) − (left), we mark two elements,indicated by green triangles (middle), for refinement. Using iterated NVB refinement we obtain the mesh T (cid:96) (right). The set (cid:101) V (cid:96) consistsof the new vertices (red) and neighboring vertices, where the corresponding vertex patches have changed (blue). of newly created vertices plus some of their neighbors, see Figure 5.2 for a visualization. Based on these sets, weconsider the space decomposition V h = L (cid:88) (cid:96) =0 (cid:88) z ∈ (cid:101) V (cid:96) V (cid:96), z with V (cid:96), z := span { ϕ (cid:96), z } , (5.5)where ϕ (cid:96), z ∈ (cid:101) S ( T (cid:96) ) is the nodal hat function with ϕ (cid:96), z ( z ) = 1 and ϕ (cid:96), z ( z (cid:48) ) = 0 for all z (cid:48) ∈ V (cid:96) \{ z } . The basicidea of this splitting is that we do a diagonal scaling only in the regions where the meshes have been refined. Wewill use local exact solvers, i.e., (cid:101) a (cid:96), z ( u (cid:96), z , v (cid:96), z ) := (cid:68) (cid:101) D ( R (cid:96), z ) T u (cid:96), z , ( R (cid:96), z ) T v (cid:96), z (cid:69) Γ for all u (cid:96), z , v (cid:96), z ∈ V (cid:96), z , where ( R (cid:96), z ) T : V (cid:96), z → V h denotes the canonical embedding operator. Let (cid:101) D h denote the Galerkin matrix of (cid:101) D with respect to the basis ( ϕ L, z ) z ∈V L of V h and define (cid:101) D (cid:96), z := (cid:101) a (cid:96), z ( ϕ (cid:96), z , ϕ (cid:96), z ). Then, the local multilevel diagonal(LMLD) preconditioner associated to the splitting (5.5) reads( B h ) − := L (cid:88) (cid:96) =0 (cid:88) z ∈ (cid:101) V (cid:96) ( R (cid:96), z ) T (cid:16) (cid:101) D (cid:96), z (cid:17) − R (cid:96), z . (5.6)We stress that this preconditioner corresponds to a diagonal scaling with respect to the local subset of vertices (cid:101) V (cid:96) on each level (cid:96) = 0 , . . . , L . Further details and the proof of the following result are found in [14, 16]. Proposition 5.4.
The splitting (5.5) together with (cid:101) a (cid:96), z ( · , · ) and the operators R T(cid:96), z satisfies the requirements ofProposition 5.1 with constants depending only on Γ and the initial triangulation T . For the additive Schwarzoperator P h := ( B h ) − (cid:101) D h , there holds in particular c (cid:68) (cid:101) Du h , u h (cid:69) Γ ≤ (cid:68) (cid:101) DP h u h , u h (cid:69) Γ ≤ C (cid:68) (cid:101) Du h , u h (cid:69) Γ for all u h ∈ V h . (5.7) The constants c , C > depend only on Γ , the initial triangulation T , and the use of NVB for refinement, i.e., T (cid:96) +1 = refine( T (cid:96) , M (cid:96) ) with arbitrary set M (cid:96) ⊆ T (cid:96) of marked elements.
20e replace the space V h in (5.1) by the refined splitting (5.5) and end up with the space decomposition V = L (cid:88) (cid:96) =0 (cid:88) z ∈ (cid:101) V (cid:96) V (cid:96), z + (cid:88) z ∈V L V pL, z . (5.8)The following Lemma 5.5 shows that the preconditioner resulting from the decomposition (5.8) is hp -stable. Theresult formalizes the observation that the combination of stable subspace decompositions leads again to a stablesubspace decomposition. It is a simple consequence of the well-known theory for additive Schwarz methods; seeSection 5.1. Therefore, details are left to the reader. Lemma 5.5.
Let V be a finite dimensional vector space, and let V j , R T V ,j , and (cid:101) a V ,j ( · , · ) for j = 0 , . . . , J be adecomposition of V in the sense of Section 5.1 that satisfies the assumptions of Proposition 5.1 with constants C , V and C , V . Consider an additional decomposition W (cid:96) , R T W ,(cid:96) and (cid:101) a W ,(cid:96) ( · , · ) with (cid:96) = 0 , . . . , L of V that also satisfiesthe requirements of Proposition 5.1 for the bilinear form (cid:101) a V , ( · , · ) with constants C , W and C , W . Define a newadditive Schwarz preconditioner as: (cid:101) B − := R T V , (cid:32) L (cid:88) (cid:96) =0 R T W ,(cid:96) (cid:101) A − W ,(cid:96) R W ,(cid:96) (cid:33) R V , + J (cid:88) j =1 R T V ,j (cid:101) A − V ,j R V ,j . This new preconditioner satisfies the assumptions of Proposition 5.1 with C = max (1 , C , W ) C , V and C =max (1 , C , W ) C , V . Theorem 5.6.
Assume that T is generated from a regular and shape-regular initial triangulation T by successiveapplication of NVB. Based on the space decomposition (5.8) define the preconditioner B − := R Th ( B h ) − R h + (cid:88) z ∈V L R T z (cid:16) (cid:101) D ph, z (cid:17) − R z . Then, for constants c , C > that depend only on Γ , T , and the use of NVB refinement, the extremal eigenvaluesof B − (cid:101) D ph satisfy c ≤ λ min ( B − (cid:101) D ph ) ≤ λ max ( B − (cid:101) D ph ) ≤ C. In particular, the condition number κ ( B − (cid:101) D ph ) is bounded independently of h and p .Proof. The proof follows from Lemma 5.2, Proposition 5.4, and Lemma 5.5.
For each vertex patch, we need to store the dense matrix (cid:16) (cid:101) D ph, z (cid:17) − . For higher polynomial orders, storingthese blocks is a significant part of the memory consumption of the preconditioner. To reduce these costs, wecan make use of the fact that the abstract additive Schwarz theory allows us to replace the local bilinear forms a ( R Ti u i , R Ti v i ) with spectrally equivalent forms, as long as they satisfy the conditions stated in Proposition 5.1.This is for example the case, if the decomposition is stable for the exact local solvers and if there exist constants c , c > c (cid:101) a i ( u i , u i ) ≤ a ( R Ti u i , R Ti u i ) ≤ c (cid:101) a i ( u i , u i ) ∀ u i ∈ V i . The new preconditioner will be based on a finite number of reference patches, for which the Galerkin matrix hasto be inverted.First we prove the simple fact that we can drop the stabilization term from (2.4) when assembling the localbilinear forms: 21 emma 5.7.
There exists a constant c > that depends only on Γ and the γ -shape regularity of T such that forany vertex patch ω z the following estimates hold: (cid:104) Du, u (cid:105) Γ ≤ (cid:68) (cid:101) Du, u (cid:69) Γ ≤ c (cid:104) Du, u (cid:105) Γ ∀ u ∈ V p z . Proof.
The first estimate is trivial, as (cid:101) D only adds an additional non-negative term. For the second inequality,we note that the functions in V p z all vanish outside of ω z and therefore V p z ∩ ker( D ) = { } . We transform to thereference patch, use the fact that (cid:98) D is elliptic on (cid:101) H / ( (cid:98) ω z ), and transform back by applying Lemma 3.6: (cid:107) u (cid:107) L ( ω z ) (cid:46) h z (cid:107) (cid:98) u (cid:107) L ( (cid:98) ω z ) (cid:46) h z (cid:107) (cid:98) u (cid:107) (cid:101) H / ( (cid:98) ω z ) (cid:46) h z (cid:68) (cid:98) D (cid:98) u, (cid:98) u (cid:69) (cid:98) ω z (cid:46) h z (cid:104) Du, u (cid:105) ω z . Thus, we can simply estimate the stabilization: α (cid:104) u, (cid:105) ω z ≤ α (cid:107) u (cid:107) L ( ω z ) (cid:107) (cid:107) L ( ω z ) ≤ Cα (cid:107) (cid:107) L ( ω z ) h z (cid:104) Du, u (cid:105) Γ . This gives the full estimate with the constant c := max (cid:16) , α (cid:107) (cid:107) L ( ω z ) Ch z (cid:17) ≤ max (cid:0) , Cα h z (cid:1) . Remark 5.8.
The proof of the previous lemma shows that this modification does not significantly affect the stabilityof the preconditioner and its effect will even vanish with h -refinement. We are now able to define the new local bilinear forms as:
Definition 5.9.
Take z ∈ V and let F z : (cid:98) ω z → ω z be the pullback mapping to the reference patch as in Defini-tion 3.5. Set (cid:101) a z ( u, v ) := h z (cid:68) (cid:98) D ( u ◦ F z ) , v ◦ F z (cid:69) (cid:98) ω z ∀ u, v ∈ V p z . (see Lemma 3.6 for the definition of (cid:98) D ). We denote the Galerkin matrix corresponding the bilinear form (cid:101) a z on thereference patch by (cid:98) D ph, ref( z ) . The above definition only needs to evaluate (cid:68) (cid:98) D (cid:98) u, (cid:98) v (cid:69) Γ on the reference patch. Since the reference patch dependsonly on the number of elements belonging to the patch, the number of blocks that need to be stored, depends onlyon the shape regularity and is independent of the number of vertices in the triangulation T . Theorem 5.10.
Assume that T is generated from a regular and shape-regular initial triangulation T by successiveapplication of NVB. The preconditioner using the local solvers from Definition 5.9 is optimal, i.e., for B − := R Th ( B h ) − R h + (cid:88) z ∈V L h − z R T z (cid:16) (cid:98) D ph, ref( z ) (cid:17) − R z , the condition number of the preconditioned system satisfies κ ( B − ˜ D ph ) ≤ C, where C > depends only on Γ , T and the use of NVB refinement. It is in particular independent of h and p .Proof. The scaling properties of (cid:104)
Du, u (cid:105) Γ were stated in Lemma 3.6. Therefore, we can conclude the argument byusing the standard additive Schwarz theory. 22
23 4 5 60 F − K (0) F − K (1) F − K (2) 1 2345 6 0 F − K (0) F − K (1) F − K (2) F z QF − K F K Figure 5.3: The mapping to the reference patch described as a combination of element maps.
When implementing the preconditioner as defined above, it is important to note that for a basis function ϕ i on ω z the transformed function ϕ i ◦ F − z does not necessarily correspond to the i -th basis function on (cid:98) ω z . Dependingon the chosen basis we may run into orientation difficulties. This can be resolved in the following way:Let z ∈ V be fixed. Choose a numbering for the vertices z i and elements K i of ω z such that adjacentelements have adjacent numbers (for example, enumerate clockwise or counter-clockwise). We also choose a similarenumeration on the reference patch and denote it as (cid:98) z i and (cid:98) K i . The enumeration is such that the reference map F z maps z i to (cid:98) z i and K i to (cid:98) K i . Let N z be the number of vertices in the patch.For elements K ⊂ ω z and K (cid:48) ⊂ (cid:98) ω z , the bases on ω z and on (cid:98) ω z are locally defined by the pullback of polynomialson the reference triangle (cid:98) K . We denote the element maps as F K : (cid:98) K → K and F (cid:48) K : (cid:98) K → K (cid:48) , respectively. Thebasis functions are then given as ϕ j := (cid:98) ϕ j ◦ F K on ω z and ψ j := (cid:98) ψ j ◦ F K (cid:48) on (cid:98) ω z . Corresponding local elementmaps do not necessarily map the same vertices of the reference element (cid:98) K to vertices with the same numbers inthe local ordering. Hence, we need to introduce another map Q : (cid:98) K → (cid:98) K that represents a vertex permutation.Then, we can write the patch-pullback restricted to K (cid:48) as F z | K (cid:48) = F K ◦ Q ◦ F − K (cid:48) (see Figure 5.3). We observe:i) For the hat function the mapping is trivial: ϕ z ◦ F z = ψ (cid:98) z .ii) For the edge basis, permuting the vertices on the reference element only changes the sign of the correspondingedge functions. Thus, we have ϕ E m j ◦ F z = ( − j ψ E m j , if the orientation of the edge in the global triangulationdoes not match the orientation of the reference patch.iii) The inner basis functions transformation under Q is not so simple. Since the basis functions all have supporton a single element we can restrict our consideration to this element and assemble the necessary basis trans-formations for all 5 permutations of vertices on the reference triangle without losing the memory advantage ofusing the reference patch. Remark 5.11.
One could also exploit the symmetry (up to a sign change) of the permutation of λ and λ in thedefinition of the inner basis functions to reduce the number of basis transformation matrices needed from 5 to 2. Figure 6.1: Adaptive meshes on the Fichera cube and for a screen problem.
6. Numerical results
The following numerical experiments confirm that the proposed preconditioners (Theorem 5.3, Theorem 5.6, andTheorem 5.10) do indeed yield a system with a condition number that is bounded uniformly in h and p , whereasthe condition number of the unpreconditioned system grows in p with a rate slightly smaller than predicted inCorollary 4.7: We observe numerically κ ∼ O ( p . ). Diagonal preconditioning appears to reduce the conditionnumber to O ( p . ). All of the following experiments were performed using the BEM++ software library ([38]; ) with the AHMED software library for H -matrix compression, [4], [5]. We used the polynomialbasis described in Section 2.2. Example 6.1 (unpreconditioned p -dependence) . We consider a quadratic screen in R (see Figure 6.1, right).We study the p -dependence of the unpreconditioned system on different uniformly refined meshes. In accordancewith the estimates of Corollary 4.7, Figure 6.2 shows that one has, depending on the mesh size h , a preasymptoticphase in which the O ( h − p ) term dominates, and an h -independent asymptotic O ( p . ) behavior. The latter isslightly better than the prediction of O ( p ) of Corollary 4.7. Example 6.2 (Fichera’s cube) . We compare the preconditioner that uses the local multilevel preconditioner for the h -part and the inexact local solvers based on the reference patches to the unpreconditioned system and to simplediagonal scaling. We consider the problem on a closed surface, namely, the surface of the Fichera cube with sidelength , and employ a stabilization parameter α = 0 . . To generate the adaptive meshes, we used NVB, wherein each step, the set of marked elements originated from a lowest order adaptive algorithm with a ZZ-type errorestimator (as described in [3]). The left part of Figure 6.1 shows an example of one of the meshes used.Figure 6.3 confirms that the condition number of the preconditioned system does not depend on the polynomialdegree of the discretization. Figure 6.4 confirms the robustness of the preconditioner with respect to the adaptiverefinement level. The unpreconditioned and the diagonally preconditioned system do not show a bad behavior withrespect to h , probably due to the already large condition number for p > . Example 6.3 (screen problem) . We consider the screen problem in R with a quadratic screen of side length 1(see Figure 6.1, right), which represents the case Γ (cid:54) = ∂ Ω and α = 0 in (2.3) , and perform the same experiments aswe did for Fichera’s cube in Example 6.2. In Figure 6.5 we again observe that the condition number is independentof the polynomial degree. Figures 6.6–6.9 demonstrate the independence of the mesh size h . O ( p . )polynomial degree λ m a x / λ m i n Mesh with 16 ElementsMesh with 256 ElementsMesh with of 4096 ElementsMesh with 16384 Elements
Figure 6.2: Comparison of the condition number of (cid:101) D ph for the screen problem on different uniform meshes (Example 6.1) O ( p . ) O ( p . )polynomial degree λ m a x / λ m i n no preconditionerdiagonal preconditioneradditive Schwarz preconditioner (Theorem 5.10) Figure 6.3: Fichera cube, condition numbers for fixed uniform mesh with 70 elements (Example 6.2).
00 1 ,
000 2 ,
000 5 , number of elements λ m a x / λ m i n no preconditionerdiagonal preconditioneradditive Schwarz preconditioner (Theorem 5.10) Figure 6.4: Fichera cube, adaptive h -refinement for p = 3 (Example 6.2). O ( p . ) O ( p )polynomial degree λ m a x / λ m i n no preconditionerdiagonal preconditioneradditive Schwarz preconditioner (Theorem 5.10) Figure 6.5: Screen problem, condition numbers for uniform mesh with 45 elements (Example 6.3). number of elements λ m a x / λ m i n no preconditionerdiagonal preconditioneradditive Schwarz preconditioner (Theorem 5.10) Figure 6.6: Screen problem, uniform h -refinement for p = 4 (Example 6.3). number of elements λ m a x / λ m i n no preconditionerdiagonal preconditioneradditive Schwarz preconditioner (Theorem 5.10) Figure 6.7: Screen problem, adaptive h -refinement for p = 1 (Example 6.3). number of elements λ m a x / λ m i n no preconditionerdiagonal preconditioneradditive Schwarz preconditioner (Theorem 5.10) Figure 6.8: Screen problem, adaptive h -refinement for p = 2 (Example 6.3). number of elements λ m a x / λ m i n no preconditionerdiagonal preconditioneradditive Schwarz preconditioner (Theorem 5.10) Figure 6.9: Screen problem, adaptive h -refinement for p = 3 (Example 6.3). λ m a x / λ m i n full vertex blocks and global inverse for e S ( T ) (Theorem 5.3)full vertex blocks and multilevel preconditioner (Theorem 5.6)multilevel for h -part and reference patches (Theorem 5.10) Figure 6.10: Comparison of the different proposed preconditioners for a fixed uniform mesh with 70 elements on the Fichera cube(Example 6.4). xample 6.4 (inexact local solvers) . We compare the different preconditioners proposed in this paper. While thenumerical experiments all show that the preconditioner is indeed robust in h and p , the constant differs if we usethe different simplifications described in the Sections 5.3 and 5.4 to the preconditioner. In Figures 6.10 and 6.11,we can observe the different constants for the geometry given by Fichera’s cube of Example 6.2. ,
000 1 ,
500 2 ,
000 2 ,
500 3 ,
000 3 ,
500 4 ,
000 4 , λ m a x / λ m i n full vertex blocks and global inverse for e S ( T ) (Theorem 5.3)full vertex blocks and multilevel preconditioner (Theorem 5.6)multilevel for h -part and reference patches (Theorem 5.10) Figure 6.11: Comparison of the different proposed preconditioners for adaptive mesh refinement on the Fichera cube with p = 2. xample 6.5 (inexact local solvers) . We continue with the geometry of Example 6.4, i.e., Fichera’s cube. Wemotivated Section 5.4 by stating the large memory requirement of the preconditioner when storing the dense localblock inverses. It can be seen in Table 6.1 that the reference patch based preconditioner resolves this issue: we presentthe memory requirements for the various approaches when excluding the memory requirement for the treatment ofthe lowest order space V h . For comparison, we included the storage requirements for the full matrix (cid:101) D ph and the H -matrix approximation with accuracy − which is denoted as D p, H h . While we still get linear growth in the numberof elements, due to some bookkeeping requirements, such as element orientation etc., which could theoreticallyalso be avoided, we observe a much reduced storage cost. For p = 3 and , degrees of freedom, the memoryrequirement is less than . of the full block storage. For p = 4 and , degrees of freedom the memoryrequirement is just . and for higher polynomial orders, this ratio would become even smaller. Comparing onlythe number of blocks that need to be stored, we see that in this particular geometry we only need to store the inversefor reference blocks. p N dof mem (cid:16) (cid:101) D ph (cid:17) /N dof mem (cid:16) (cid:101) D p, H h (cid:17) /N dof mem (cid:0) B − (cid:1) /N dof mem (cid:0) B − (cid:1) /N dof [KB] [KB] [KB] [KB]2 98 0 .
765 62 0 .
769 45 0 .
094 547 0 .
039 2222 298 2 . . .
098 259 0 .
027 3182 986 7 . . .
100 25 0 .
020 5222 3558 27 .
797 15 .
006 0 .
100 68 0 .
018 3942 8950 69 .
922 19 .
817 0 .
100 78 0 .
017 9023 218 1 . . .
313 14 0 .
083 7873 668 5 . . .
325 08 0 .
041 2613 2216 17 .
312 10 .
921 0 .
332 0 .
017 8953 5969 46 .
633 16 .
452 0 . .
011 5563 16 310 127 .
42 22 .
738 0 .
332 74 0 .
009 182 43 20 135 157 . .
04 0 .
333 72 0 .
008 922 24 386 3 . . .
670 86 0 .
169 734 1954 15 .
266 10 .
879 0 .
702 33 0 .
048 294 5634 44 .
016 17 .
39 0 .
710 85 0 .
019 6194 14 226 111 .
14 21 .
209 0 .
713 64 0 .
010 4244 35 794 279 .
64 27 .
596 0 .
714 28 0 .
006 790 85 602 4 . . . .
293 245 1852 14 .
469 14 .
476 1 . .
129 455 8802 68 .
766 68 .
772 1 . .
029 4595 16 577 129 .
51 129 . . .
016 9615 45 302 353 .
92 353 .
79 1 . .
007 989 85 55 927 436 .
93 436 .
78 1 . .
007 006 2
Table 6.1: Comparison of the memory requirement relative to the number of degrees of freedom N dof between storing the full blockstructure and the reference block based preconditioner from Section 5.4 (Example 6.5). cknowledgments: The research was supported by the Austrian Science Fund (FWF) through the doctoralschool “Dissipation and Dispersion in Nonlinear PDEs” (project W1245, A.R.) and “Optimal Adaptivity for BEMand FEM-BEM Coupling” (project P27005, T.F, D.P.). T.F. furthermore acknowledges funding through theInnovative Projects Initiative of Vienna University of Technology and the CONICYT project “Preconditionedlinear solvers for nonconforming boundary elements” (grant FONDECYT 3150012).
Appendix
From [32, Prop. 2.8] for the stiffness matrix and from classical estimates for the quadrature weights of theGauß-Lobatto quadrature, we have on the reference square (cid:98) Sp − (cid:107) u (cid:107) (cid:96) (cid:46) (cid:107) u (cid:107) H ( (cid:98) S ) (cid:46) p (cid:107) u (cid:107) (cid:96) , p − (cid:107) u (cid:107) (cid:96) (cid:46) (cid:107) u (cid:107) L ( (cid:98) S ) (cid:46) p − (cid:107) u (cid:107) (cid:96) ∀ u ∈ Q p ( (cid:98) S )For quasi-uniform meshes, we therefore obtain h − p (cid:107) u (cid:107) L (Γ) (cid:46) (cid:107) u (cid:107) (cid:96) (cid:46) h − p (cid:107) u (cid:107) L (Γ) . Furthermore, we have (cid:107) J h u (cid:107) (cid:96) (cid:46) h − (cid:107) J h u (cid:107) L (Γ) (cid:46) h − (cid:107) u (cid:107) L (Γ) (cid:46) p − (cid:107) u (cid:107) (cid:96) , (cid:107) u − J h u (cid:107) (cid:96) (cid:46) (cid:88) K ∈T (cid:107) ( u − J h u ) | K ) (cid:107) (cid:96) (cid:46) p (cid:88) K ∈T (cid:107) (cid:98) u − (cid:100) J h u (cid:107) H ( (cid:98) S ) (cid:46) p (cid:88) K ∈T | u − J h u | H ( K ) + h − (cid:107) u − J h u (cid:107) L ( K ) (cid:46) p (cid:107) u (cid:107) H (Γ) . We obtain (cid:107) u (cid:107) L (Γ) (cid:46) h p − (cid:107) u (cid:107) (cid:96) , (cid:107) u (cid:107) H (Γ) = (cid:107) u (cid:107) L (Γ) + | u | H (Γ) (cid:46) h p − (cid:107) u (cid:107) (cid:96) + p (cid:107) u (cid:107) (cid:96) (cid:46) p (cid:107) u (cid:107) (cid:96) , so that interpolation yields (cid:107) u (cid:107) (cid:101) H / (Γ) (cid:46) hp − / (cid:107) u (cid:107) (cid:96) . For the converse estimate, we observe (cid:107) u − J h u (cid:107) (cid:96) (cid:46) p (cid:107) u (cid:107) H (Γ) , (cid:107) u − J h u (cid:107) (cid:96) (cid:46) h − p (cid:107) u (cid:107) L (Γ) , so that an interpolation argument (which we assume to be admissible!) produces (cid:107) u − J h u (cid:107) (cid:96) (cid:46) p h − (cid:107) u (cid:107) (cid:101) H / (Γ) . Hence, (cid:107) u (cid:107) (cid:96) (cid:46) (cid:107) u − J h u (cid:107) (cid:96) + (cid:107) J h u (cid:107) (cid:96) (cid:46) p h − (cid:107) u (cid:107) (cid:101) H / (Γ) + h − (cid:107) u (cid:107) L (Γ) (cid:46) (cid:0) p h − + h − (cid:1) (cid:107) u (cid:107) (cid:101) H / (Γ) . Putting things together, we get p / h − (cid:107) u (cid:107) (cid:101) H / (Γ) (cid:46) (cid:107) u (cid:107) (cid:96) (cid:46) (cid:0) p h − + h − (cid:1) (cid:107) u (cid:107) (cid:101) H / (Γ) , which in turn gives the condition number estimate κ ( (cid:101) D ph ) (cid:46) p / + h − p − / . For the H -condition number, we note the estimates (cid:107) u (cid:107) H (Γ) = | u | H (Γ) + (cid:107) u (cid:107) L (Γ) (cid:46) p (cid:107) u (cid:107) (cid:96) + h p − (cid:107) u (cid:107) (cid:96) (cid:46) p (cid:107) u (cid:107) (cid:96) , (cid:107) u (cid:107) (cid:96) (cid:46) (cid:107) u − J h u (cid:107) (cid:96) + (cid:107) J h u (cid:107) (cid:96) (cid:46) p (cid:107) u (cid:107) H (Γ) + h − (cid:107) u (cid:107) L (Γ) (cid:46) (cid:0) p + h − (cid:1) (cid:107) u (cid:107) H (Γ) , so that we get p − (cid:107) u (cid:107) H (Γ) (cid:46) (cid:107) u (cid:107) (cid:96) (cid:46) (cid:0) p + h − (cid:1) (cid:107) u (cid:107) H (Γ) eferencesReferences [1] M. Ainsworth, B. Guo, An additive Schwarz preconditioner for p -version boundary element approximation of the hypersingularoperator in three dimensions, Numer. Math. 85 (2000) 343–366.[2] M. Ainsworth, W. McLean, Multilevel diagonal scaling preconditioners for boundary element equations on locally refined meshes,Numer. Math. 93 (3) (2003) 387–413.[3] M. Aurada, M. Feischl, T. F¨uhrer, M. Karkulik, D. Praetorius, Energy norm based error estimators for adaptive BEM forhypersingular integral equations, Appl. Numer. Math. (2015) published online first.[4] M. Bebendorf, Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems, vol. 63 of Lect. NotesComput. Sci. Eng., Springer-Verlag, 2008, iSBN 978-3-540-77146-3.[5] M. Bebendorf, Another software library on hierarchical matrices for elliptic differential equations (AHMED), http://bebendorf.ins.uni-bonn.de/AHMED.html (Jun. 2014).[6] C. Bernardi, V. Girault, A local regularization operator for triangular and quadrilateral finite elements, SIAM J. Numer. Anal.35 (5) (1998) 1893–1916.[7] C. Bernardi, Y. Maday, Spectral methods, in: P. Ciarlet, J. Lions (eds.), Handbook of Numerical Analysis, Vol. 5, North Holland,Amsterdam, 1997.[8] J. M. Cascon, C. Kreuzer, R. H. Nochetto, K. G. Siebert, Quasi-optimal convergence rate for an adaptive finite element method,SIAM J. Numer. Anal. 46 (5) (2008) 2524–2550.[9] M. Dubiner, Spectral methods on triangles and other domains, J. Sci. Comp. 6 (1991) 345–390.[10] R. Falk, R. Winther, The bubble transform: A new tool for analysis of finite element methods, Found. Comput. Math. (2015)1–32.[11] M. Feischl, T. F¨uhrer, M. Karkulik, J. M. Melenk, D. Praetorius, Quasi-optimal convergence rates for adaptive boundary elementmethods with data approximation, part I: weakly-singular integral equation, Calcolo 51 (2014) 531–562.[12] M. Feischl, T. F¨uhrer, M. Karkulik, J. M. Melenk, D. Praetorius, Quasi-optimal convergence rates for adaptive boundary elementmethods with data approximation. Part II: Hyper-singular integral equation, Electron. Trans. Numer. Anal. 44 (2015) 153–176.[13] M. Feischl, T. F¨uhrer, D. Praetorius, Adaptive FEM with optimal convergence rates for a certain class of non-symmetric andpossibly non-linear problems, SIAM J. Numer. Anal. 52(2) (2014) 601 – 625.[14] M. Feischl, T. F¨uhrer, D. Praetorius, E. P. Stephan, Efficient additive Schwarz preconditioning for hypersingular integral equationson locally refined triangulations, ASC Report , Vienna University of Technology.[15] M. Feischl, M. Karkulik, J. M. Melenk, D. Praetorius, Quasi-optimal convergence rate for an adaptive boundary element method,SIAM J. Numer. Anal. 51 (2) (2013) 1327–1348.[16] T. F¨uhrer, Zur Kopplung von finiten Elementen und Randelementen, Ph.D. thesis, Vienna University of Technology, in German(2014).[17] T. Gantumur, Adaptive boundary element methods with convergence rates, Numer. Math. 124 (3) (2013) 471–516.[18] N. Heuer, Additive Schwarz methods for indefinite hypersingular integral equations in R —the p -version, Appl. Anal. 72 (3-4)(1999) 411–437.[19] R. Hiptmair, S. Mao, Stable multilevel splittings of boundary edge element spaces, BIT 52 (3) (2012) 661–685.[20] R. Hiptmair, H. Wu, W. Zheng, Uniform convergence of adaptive multigrid methods for elliptic problems and Maxwell’s equations,Numer. Math. Theory Methods Appl. 5 (3) (2012) 297–332.[21] G. C. Hsiao, W. L. Wendland, Boundary integral equations, vol. 164 of Applied Mathematical Sciences, Springer-Verlag, Berlin,2008.[22] N. Hu, X.-Z. Guo, I. N. Katz, Bounds for eigenvalues and condition numbers in the p -version of the finite element method, Math.Comp. 67 (224) (1998) 1423–1450.[23] M. Karkulik, J. M. Melenk, A. Rieder, Optimal additive Schwarz methods for the p -BEM: the hypersingular integral operator (inpreparation).[24] M. Karkulik, D. Pavlicek, D. Praetorius, On 2D newest vertex bisection: Optimality of mesh-closure and H -stability of L -projection, Constr. Approx. 38 (2013) 213–234.[25] G. E. Karniadakis, S. J. Sherwin, Spectral/ hp element methods for CFD, Numerical Mathematics and Scientific Computation,Oxford University Press, New York, 1999.[26] T. Koornwinder, Two-variable analogues of the classical orthogonal polynomials, in: Theory and application of special functions(Proc. Advanced Sem., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1975), Academic Press, New York, 1975, pp. 435–495.Math. Res. Center, Univ. Wisconsin, Publ. No. 35.[27] P.-L. Lions, On the Schwarz alternating method. I, in: First International Symposium on Domain Decomposition Methods forPartial Differential Equations (Paris, 1987), SIAM, Philadelphia, PA, 1988, pp. 1–42.[28] M. Maischak, A multilevel additive Schwarz method for a hypersingular integral equation on an open curve with graded meshes,Appl. Numer. Math. 59 (9) (2009) 2195–2202.[29] J.-F. Maitre, O. Pourquier, Condition number and diagonal preconditioning: comparison of the p -version and the spectral elementmethods, Numer. Math. 74 (1) (1996) 69–84.[30] A. M. Matsokin, S. V. Nepomnyaschikh, A Schwarz alternating method in a subspace, Soviet Math. 29(10) (1985) 78–84.[31] W. McLean, Strongly elliptic systems and boundary integral equations, Cambridge University Press, Cambridge, 2000.
32] J. M. Melenk, On condition numbers in hp -FEM with Gauss-Lobatto-based shape functions, J. Comput. Appl. Math. 139 (1)(2002) 21–48.[33] P. Oswald, Interface preconditioners and multilevel extension operators, in: Eleventh International Conference on Domain De-composition Methods (London, 1998), DDM.org, Augsburg, 1999, pp. 97–104.[34] L. F. Pavarino, Additive Schwarz methods for the p -version finite element method, Numer. Math. 66 (4) (1994) 493–515.[35] S. A. Sauter, C. Schwab, Boundary element methods, vol. 39 of Springer Series in Computational Mathematics, Springer-Verlag,Berlin, 2011, translated and expanded from the 2004 German original.[36] J. Sch¨oberl, J. M. Melenk, C. Pechstein, S. Zaglmayr, Additive Schwarz preconditioning for p -version triangular and tetrahedralfinite elements, IMA J. Numer. Anal. 28 (1) (2008) 1–24.[37] C. Schwab, p - and hp -finite element methods, Numerical Mathematics and Scientific Computation, The Clarendon Press, OxfordUniversity Press, New York, 1998, theory and applications in solid and fluid mechanics.[38] W. ´Smigaj, T. Betcke, S. Arridge, J. Phillips, M. Schweiger, Solving boundary integral problems with BEM++, ACM Trans.Math. Softw. 41 (2) (2015) 6:1–6:40.[39] O. Steinbach, Numerical approximation methods for elliptic boundary value problems, Springer, New York, 2008, finite andboundary elements, Translated from the 2003 German original.[40] R. Stevenson, The completion of locally refined simplicial partitions created by bisection, Math. Comp. 77 (261) (2008) 227–241.[41] L. Tartar, An introduction to Sobolev spaces and interpolation spaces, vol. 3 of Lecture Notes of the Unione Matematica Italiana,Springer, Berlin, 2007.[42] A. Toselli, O. Widlund, Domain decomposition methods—algorithms and theory, vol. 34 of Springer Series in ComputationalMathematics, Springer-Verlag, Berlin, 2005.[43] T. Tran, E. P. Stephan, Additive Schwarz methods for the h -version boundary element method, Appl. Anal. 60 (1-2) (1996) 63–84.[44] T. Tran, E. P. Stephan, P. Mund, Hierarchical basis preconditioners for first kind integral equations, Appl. Anal. 65 (3-4) (1997)353–372.[45] H. Triebel, Interpolation theory, function spaces, differential operators, 2nd ed., Johann Ambrosius Barth, Heidelberg, 1995.[46] T. von Petersdorff, Randwertprobleme der Elastizit¨atstheorie f¨ur Polyeder - Singularit¨aten und Approximation mit Randelement-methoden, Ph.D. thesis, Technische Hochschule Darmstadt, in German (1989).[47] H. Wu, Z. Chen, Uniform convergence of multigrid V-cycle on adaptively refined finite element meshes for second order ellipticproblems, Sci. China Ser. A 49 (10) (2006) 1405–1429.[48] J. Xu, L. Chen, R. H. Nochetto, Optimal multilevel methods for H (grad), H (curl), and H (div) systems on graded and unstructuredgrids, in: Multiscale, nonlinear and adaptive approximation, Springer, Berlin, 2009, pp. 599–659.[49] X. Xu, H. Chen, R. H. W. Hoppe, Optimality of local multilevel methods on adaptively refined meshes for elliptic boundary valueproblems, J. Numer. Math. 18 (1) (2010) 59–90.[50] S. Zaglmayr, High order finite element methods for electromagnetic field computation, Ph.D. thesis, Johannes Kepler UniversityLinz (2006).[51] X. Zhang, Multilevel Schwarz methods, Numer. Math. 63 (4) (1992) 521–539.(div) systems on graded and unstructuredgrids, in: Multiscale, nonlinear and adaptive approximation, Springer, Berlin, 2009, pp. 599–659.[49] X. Xu, H. Chen, R. H. W. Hoppe, Optimality of local multilevel methods on adaptively refined meshes for elliptic boundary valueproblems, J. Numer. Math. 18 (1) (2010) 59–90.[50] S. Zaglmayr, High order finite element methods for electromagnetic field computation, Ph.D. thesis, Johannes Kepler UniversityLinz (2006).[51] X. Zhang, Multilevel Schwarz methods, Numer. Math. 63 (4) (1992) 521–539.