Adaptive 2D IGA boundary element methods

Engineering Analysis with Boundary Elements 62 (2016) 141–153

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Adaptive 2D IGA boundary element methods Michael Feischl a, Gregor Gantner b,n, Alexander Haberl b, Dirk Praetorius b a b

School of Mathematics and Statistics, University of New South Wales, The Red Centre, Sydney, NSW 2052, Australia Institute for Analysis and Scientific Computing, Vienna University of Technology, Wiedner Hauptstraße 8-10, 1040 Wien, Austria

art ic l e i nf o

a b s t r a c t

Article history: Received 27 March 2015 Received in revised form 7 October 2015 Accepted 7 October 2015 Available online 5 November 2015

We derive and discuss a posteriori error estimators for Galerkin and collocation IGA boundary element methods for weakly singular integral equations of the first-kind in 2D. While recent own work considered the Faermann residual error estimator for Galerkin IGA boundary element methods, the present work focuses more on collocation and weighted-residual error estimators, which provide reliable upper bounds for the energy error. Our analysis allows piecewise smooth parametrizations of the boundary, local mesh-refinement, and related standard piecewise polynomials as well as NURBS. We formulate an adaptive algorithm which steers the local mesh-refinement and the multiplicity of the knots. Numerical experiments show that the proposed adaptive strategy leads to optimal convergence, and related IGA boundary element methods are superior to standard boundary element methods with piecewise polynomials. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Isogeometric analysis Boundary element method Collocation A posteriori error estimate Adaptive mesh-refinement

1. Introduction 1.1. Isogeometric analysis The central idea of isogeometric analysis (IGA) is to use the same ansatz functions for the discretization of the partial differential equation at hand as for the representation of the problem geometry. Usually, Ω is represented in computer aided design (CAD) by means of NURBS, hierarchical splines, or T-splines. This concept, invented in [18] for finite element methods (IGAFEM), has proved very fruitful in applications [18,25]; see also the monograph [6]. Since CAD directly provides a parametrization of the boundary ∂Ω, this makes the boundary element method (BEM) the most attractive numerical scheme, if applicable (i.e., provided that the fundamental solution of the differential operator is explicitly known). Isogeometric BEM (IGABEM) has first been considered in [23] for 2D resp. [27] for 3D. While standard BEM with piecewise polynomials is well-studied in the literature, cf. the monographs [26,28] and the references therein, the numerical analysis of IGABEM in essentially open. We refer to [25,29,24,21] for numerical experiments, to [30] for fast IGABEM with H-matrices, and to [17] for some quadrature analysis. A posteriori error estimation has first been considered for Galerkin IGABEM in our recent work n

Corresponding author. E-mail addresses: [email protected] (M. Feischl), [email protected] (G. Gantner), [email protected] (A. Haberl), [email protected] (D. Praetorius). http://dx.doi.org/10.1016/j.enganabound.2015.10.003 0955-7997/& 2015 Elsevier Ltd. All rights reserved.

[15]. In the present work, we extend the latter result to collocation IGABEM which is preferred in practice for its simpler assembly of the stiffness matrix. 1.2. Model problem Let Ω
R2 be a Lipschitz domain and Γ D∂Ω be a compact, piecewise smooth part of the boundary with finitely many connected components. For a given right-hand side f, we consider the weakly singular boundary integral equation Z 1 log j x yj ϕðyÞ dy ¼ f ðxÞ on Γ ð1:1Þ V ϕðxÞ ≔  2π Γ associated with the 2D Laplacian; see Section 2 below for the mathematical setting and the definition of the problem related energy norm jjjjjj. With some discrete ansatz space X h
L2 ðΓ Þ, the Galerkin BEM computes the unique solution ϕh A X h of the discrete variational formulation Z Z V ϕh ψ h dx ¼ f ψ h dx for all ψ h A X h : ð1:2Þ Γ

Γ

Note that X h
L2 ðΓ Þ ensures V ϕh A CðΓ Þ. The collocation BEM computes ϕh A X h such that V ϕh ðxj Þ ¼ f ðxj Þ

for all xj A fx1 ; …; xNcol g;

ð1:3Þ

where the xj are appropriately chosen collocation points with Ncol ≔ dim X h ; see Section 2.8. In either case (1.2)–(1.3), ϕh is computed by solving a linear system of equations.

142

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

1.3. A posteriori error estimation for Galerkin IGABEM We assume that X h is associated to some partition T h of Γ into a set of connected segments. For each vertex z of T h , let ωh ðzÞ ≔ ⋃fT A T h : z A Tg denote the node patch. If X h is sufficiently rich (e.g., X h contains certain splines or NURBS), it is proved in [15] that Galerkin BEM guarantees reliability and efficiency 0 11=2 X 1 2A @ C rel jjj ϕ  ϕh jjj r ηh ≔ ηh ðzÞ r C eff jjj ϕ  ϕh jjj zANh

ð1:4aÞ with X h -independent constants C eff ; C rel 4 0. Here, r h ≔ f  V ϕh denotes the residual and Z Z j r h ðxÞ  r h ðyÞj 2 ηh ðzÞ2 ≔ dy dx ð1:4bÞ j x  yj 2 ωh ðzÞ ωh ðzÞ is some Sobolev–Slobodeckij seminorm, i.e., the unknown BEM energy error is controlled by some computable a posteriori error estimator ηh. Estimate (1.4) has first been proved by Faermann [12] for closed Γ ¼ ∂Ω and standard spline spaces X h based on the arclength parametrization. Her result is generalized in [15] to a more general setting which also includes isogeometric analysis. We note that [12,15] show that the efficiency estimate ηh r C eff jjj ϕ  ϕh jjj holds even independently of the discretization and, in particular, for collocation. 1.4. A posteriori error estimation for collocation IGABEM In the present manuscript, we focus on the weighted-residual error estimator which has first been proposed in [9,3] for standard BEM with piecewise polynomials and polygonal Γ. We prove that for Galerkin IGABEM (1.2) as well as collocation IGABEM (1.3), there holds the upper bound 0 11=2 X 1 2A @ C rel jjj ϕ  ϕh jjj r μh ≔ μh ðzÞ ; ð1:5aÞ zANh

with an X h -independent constant C rel 4 0. Here, r h ≔ f  V ϕh is again the residual and Z μh ðzÞ2 ≔ j ωh ðzÞj j r 0h ðxÞj 2 dx ð1:5bÞ

IGABEM leads to essentially the same convergence behavior as Galerkin IGABEM, so that an adaptive collocation IGABEM may be favorable in practice. Section 4 recalls the precise statement of (1.4) from [15] and gives a proof of (1.5)–(1.7). The concluding Section 5 comments on our overall findings, open questions, and future research.

2. Preliminaries In this section, we collect the main assumptions on the boundary and its discretization and introduce the BEM ansatz spaces. Further details on Sobolev spaces and the functional analytic setting of weakly singular integral equations are found, e.g., in the monographs [19,20,26] and the references therein. Throughout, j  j denotes the absolute value of scalars, the Euclidean norm of vectors in R2 , the measure of a set in R (e.g., the length of an interval), or the arclength of a curve in R2 . The respective meaning will be clear from the context. We write A ≲ B to abbreviate A r cB with some constant c 4 0 which is clear from the context. Moreover A C B abbreviates A ≲ B ≲ A. 2.1. Function spaces For any measurable subset ω D Γ resp. any interval ω DR, L2 ðωÞ denotes the Lebesgue space of all square integrable functions with corresponding norm Z j uðxÞj 2 dx: ð2:1Þ J u J 2L2 ðωÞ ≔ ω

If u A L2 ðωÞ is differentiable along the arc, u0 denotes the arclength derivative. Define the Sobolev space H 1 ðωÞ ≔ fu A L2 ðωÞ : u0 A L2 ðωÞg with corresponding norm J u J 2H1 ðωÞ ≔ J u J 2L2 ðωÞ þ j uj 2H1 ðωÞ ; Z j uj 2H1 ðωÞ ≔

ω

1 C eff ηh r jjj ϕ  ϕh jjj r C rel μh

ð1:6Þ

which, however, involves different error estimators. In addition to the global relation of the error estimators ηh and μh and independently of the discretization, we prove

ηh ðzÞ rC loc μh ðzÞ for all vertices z of T h ; where C loc 4 0 depends only on

ð1:7Þ

Γ.

1.5. Outline Section 2 recalls the functional analytic framework, provides the assumptions on Γ and its parametrization γ, introduces the ansatz spaces, and presents an adaptive algorithm which is capable to control and adapt the multiplicity of the nodes as well as the local mesh-size (Algorithm 2.2). Section 3 provides the numerical evidence that the proposed adaptive IGABEM is superior to IGABEM with uniform mesh-refinement as well as to adaptive standard BEM with piecewise polynomials. Moreover, we observe that collocation

j u0 ðxÞj 2 dx:

ð2:2bÞ

Furthermore, define the Sobolev space H 1=2 ðωÞ ≔ fu A L2 ðωÞ : J u J H1=2 ðωÞ o 1g with corresponding norm J u J 2H1=2 ðωÞ ≔ J u J 2L2 ðωÞ þ j uj 2H1=2 ðωÞ ;

ωh ðzÞ

is a weighted H1-seminorm, where ðÞ0 denotes the arc-length derivative and j ωh ðzÞj is the length of the node patch. For collocation BEM, we thus control the energy error by

ð2:2aÞ

Z Z j uj 2H1=2 ðωÞ ≔

ω ω

j uðxÞ  uðyÞj 2 dy dx: j x  yj 2

ð2:3aÞ ð2:3bÞ

  1=2

ðωÞ, where duality is understood The dual space of H 1=2 ðωÞ is H with respect to the extended L2 ðωÞ-scalar product, i.e., for u A H 1=2 ðωÞ and ϕ A L2 ðωÞ, it holds Z ð2:4Þ 〈u; ϕ〉ω ¼ uðxÞϕðxÞ dx: ω

  1=2

ðΓ Þ form a Gelfand triple and We note that H 1=2 ðΓ Þ
L2 ðΓ Þ
H all inclusions are dense and compact. Amongst other equivalent definitions of H 1=2 ðωÞ are the characterization as trace space of functions in H 1 ðΩÞ as well as equivalent interpolation techniques. All these definitions provide the same space but different norms, where norm equivalence constants depend only on ω. 2.2. Weakly singular integral equation The operator V from (1.1) extends to a linear and continuous   1=2

ðΓ Þ-H 1=2 ðΓ Þ with additional stability V : L2 operator V : H 1 ðΓ Þ-H ðΓ Þ, see [8] resp. the monographs [20,19]. We additionally

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

suppose that V is even an elliptic isomorphism, which is satisfied, e.g., if diamðΩÞ o 1. In particular, 〈VðÞ; ðÞ〉Γ is thus a scalar product   1=2

on H

ðΓ Þ, and the induced energy norm   1=2

jjj ψ jjj 2 ≔ 〈V ψ ; ψ 〉Γ

for ψ A H

ðΓ Þ

ð2:5Þ

  1=2

ðΓ Þ. is an equivalent norm on H Given f A H 1=2 ðΓ Þ, the weakly singular integral equation (1.1) is   1=2

equivalently reformulated in variational form: Find ϕ A H such that 〈V ϕ; ψ 〉Γ ¼ 〈f ; ψ 〉Γ

  1=2

for all ψ A H

ðΓ Þ:

ðΓ Þ

  1=2

2.4. Discretization of boundary For the discretization, let T h ¼ fT 1 ; …; T n g be a partition of Γ into compact and connected segments Tj. The endpoints of the elements of T h form the set of nodes ( fzj : j ¼ 1; …; ng for Γ ¼ ∂Ω; Nh ¼ fzj : j ¼ 0; …; ng for Γ ⫋ ∂Ω: The arclength of each element T A T h is denoted by hT. Moreover, the shape regularity constant is defined by

κ ðT h Þ ≔ maxfhT =h0T : T; T 0 A T h ; T \ T 0 a ∅g: ð2:6Þ

The Lax–Milgram lemma thus applies and proves existence and uniqueness of the solution ϕ A H

143

ðΓ Þ of (2.6) resp. (1.1).

For Γ ¼ ∂Ω, we extend the nodes, elements, and their length periodically. Moreover, we suppose max hT r j Γ j =4

T AT h

for Γ ¼ ∂Ω:

ð2:9Þ

2.3. Parametrization of boundary

2.5. Discretization of parameter domain

Let Γ ¼ ⋃i Γ i D∂Ω be decomposed into its finitely many connected components Γi. The Γi are compact and piecewise smooth. Note that this yields the existence of some constant c 4 0 such that j x  yj Z c 4 0 for all x A Γ i , y A Γ j , and ia j. Together with j uðxÞ  uðyÞj 2 r 2j uðxÞj 2 þ 2j uðyÞj 2 , we derive the estimate XZ Z j uðxÞ  uðyÞj 2 X X dy dx ≲ J u J 2L2 ðΓ Þ þ J u J 2L2 ðΓ Þ C J u J 2L2 ðΓ Þ : 2 i j j x  yj Γi Γj i;j i j

Given γ : ½a; b-Γ , the partition T h induces a partition T h ¼ fT 1 ; …; T n g of the parameter domain ½a; b. Let a ¼ z 0 o z 1 o ⋯ o z n ¼ b be the endpoints of the elements of T h . We assume T j ¼ ½z j  1 ; z j , γ ðTj Þ ¼ T j , and γ ðz j Þ ¼ zj . We define ( fz j : j ¼ 1; …; ng for Γ ¼ ∂Ω; N h ¼ fz j : j ¼ 0; …; ng for Γ ⫋ ∂Ω:

iaj

This results in norm equivalence X XZ Z j uðxÞ  uðyÞj 2 J u J 2H1=2 ðΓ Þ þ dy dx J u J 2H1=2 ðΓ Þ ¼ i j x  yj 2 Γi Γj i;j i C

X i

0 0 κ ðTh Þ ≔ maxfhT =hT 0 : T ; T A Th ; γ ðT Þ \ γ ðT Þ a ∅g:

Note that κ ðT h Þ C κ ðTh Þ, where the hidden constants depend only on the parametrization γ.

iaj

J u J 2H1=2 ðΓ Þ : i

2.6. B-splines and NURBS in the parameter domain

The usual piecewise polynomial and NURBS basis functions have connected support and are hence supported by some single Γi each. Without loss of generality and to ease the mathematical proofs, we may therefore assume that Γ is connected. All results remain valid for non-connected Γ. We assume that either Γ ¼ ∂Ω is parametrized by a closed continuous and piecewise two times continuously differentiable path γ : ½a; b-Γ such that the restriction γ j ½a;bÞ is even bijective, or that Γ ⫋ ∂Ω is parametrized by a bijective continuous and piecewise two times continuously differentiable path γ : ½a; b-Γ . For Γ ¼ ∂Ω, we denote the ðb  aÞperiodic extension to R also by γ. 0 For the left and right derivative of γ, we assume that γ ℓ ðtÞ a 0 for 0 r t A ða; b and γ ðtÞ a0 for t A ½a; bÞ. Moreover we assume that 0 0 γ ℓ ðtÞ þ cγ r ðtÞ a0 for all c 4 0 and t A ½a; b resp. t A ða; bÞ. By γ L : ½0; L-Γ , we denote the arclength parametrization, i.e., j γ 0Lℓ ðtÞj ¼ 1 ¼ j γ 0Lr ðtÞj , and its periodic extension. Then, elementary differential geometry yields bi-Lipschitz continuity j γ L ðsÞ  γ L ðtÞj rC Γ j s  tj ( j s tj r 34L for all s; t A R; with s a t A ½0; L

C Γ 1 r

C Γ 1 j u○γ L j H1=2 ðIÞ r j uj H1=2 ðγ

L ðIÞÞ

 ≔ ðt i Þ We consider knots K i A Z on R with t i  1 r t i for i A Z and limi- 7 1 t i ¼ 7 1. For the multiplicity of any knot ti, we write #t i . We denote the corresponding set of nodes N ≔ ft i : i A Zg ¼ fz j : jA Zg with z j  1 o z j for all jA Z. For i A Z and p A N0 , the i-th B-Spline of degree p is defined inductively by ( χ ½ti  1 ;ti Þ for p ¼ 0;  K Bi;p ≔ Bi;p ≔ βi  1;p Bi;p  1 þ ð1  βi;p ÞBi þ 1;p  1 for p 4 0; where, for t A R, 8 t  ti > < t i þ p  ti βi;p ðtÞ ≔ > :0

for t i a t i þ p ; for t i ¼ t i þ p :

We collect some basic properties of B-splines from [10]: Lemma 2.1 (de Boor [10, Theorem6, Section 2 and pp. 9–10]). For p A N0 , the following assertions hold: (i) Let I ¼ ½a; bÞ be a finite interval. Then,

for Γ ¼ ∂Ω;

for Γ ⫋ ∂Ω;

ð2:7Þ

see, e.g., [16, Lemma 2.1] for the proof for Γ ¼ ∂Ω which even simplifies for Γ ⫋ ∂Ω. Let I D ½a; b. Suppose j Ij r 34L for Γ ¼ ∂Ω. Then, (2.7) implies

for all u A H 1=2 ðΓ Þ.

The length of each T A T h is denoted by hT . Moreover, we define the shape regularity constant on ½a; b by

r C Γ j u○γ L j H1=2 ðIÞ

ð2:8Þ

fBi;p j I : iA Z; Bi;p j I a 0g

ð2:10Þ

is a basis for the space of all right-continuous N piecewise polynomials of degree lower or equal p on I and which are, at each knot ti, p  #t i times continuously differentiable if p  #t i Z 0. (ii) For i A Z, Bi;p vanishes outside the interval ½t i  1 ; t i þ p Þ. It is positive on the open interval ðt i  1 ; t i þ p Þ. (iii) For iA Z, Bi;p is completely determined by the p þ 2 knots t i  1 ; …; t i þ p .

144

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

(iv) The B-splines of degree p form a locally finite partition of unity, P i.e., i A Z Bi;p ¼ 1 on R.□  ¼ ðt i Þ In addition to the knots K i A Z , we consider weights W ≔ ðwi Þi A Z with wi 4 0. For i A Z and p A N0 , we define the i-th non-uniform rational B-Spline (NURBS) of degree p 

RK;W ≔ Ri;p ≔ P i;p

wi Bi;p

ℓ A Z wℓ Bℓ;p

:

ð2:11Þ

Note that the denominator is positive and locally finite. For any p A N0 , we define the vector spaces ( ) X K  ≔ ai B : ai A R ; S p ðKÞ i;p

p

c ðK  0 ; W 0 Þ
L2 ðΓ Þ
H  1=2 ðΓ Þ: N ð2:12Þ

iAZ

 WÞ ≔ N ðK; p

( ) X K;W  ai Ri;p : ai A R :

ð2:13Þ

iAZ

ð2:20Þ

Choose an error estimator ϱh A fηh ; μh g. The nodal contributions ϱh  h to ðzÞ from (1.4) resp. (1.5) are used to steer knot insertion from K  H . The new weights W H are uniquely chosen the following knots K such that the denominator of the NURBS functions does not change. In particular, this implies nestedness c p ðK c p ðK  h; WhÞ D N  H ; WH Þ N

2.7. NURBS on the boundary For Γ ¼ ∂Ω, each node z A N h has a multiplicity #z r p þ 1. This  h ¼ ðt i ÞN induces a sequence of non-decreasing knots K i ¼ 1 on ða; b. Let W h ¼ ðwi ÞN be a sequence of weights on these knots. We i¼1 extend the knot sequence ðb  aÞperiodically to ðt i Þi A Z and the weight sequence to ðwi Þi A Z by wN þ i ≔ wi for i A Z. For the exten h and W h . We set ded sequences, we also write K c p ðK  h ; W h Þ ≔ N p ðK  h ; W h Þj ½a;bÞ ○γ j  1 : N ½a;bÞ

ð2:14Þ

For Γ ⫋ ∂Ω, each node z A N h has a multiplicity #z rp þ 1 such that #z 0 ¼ #z n ¼ p þ 1. This induces a sequence of nonN p  h ¼ ðt i ÞN decreasing knots K i ¼  p on ½a; b. Let W h ¼ ðwi Þi ¼ 1  p be a h¼ sequence of weights. We extend the sequences arbitrarily to K ðt i Þi A Z with t i r t i þ 1 for i A Z, a 4 t i -  1 for i o p, and b o t i 1 for i 4 N, and W h ¼ ðwi Þi A Z with wi 40. We set p

c ðK  h ; W h Þ ≔ N p ðK  h ; W h Þj ½a;b ○γ  1 : N

ð2:15Þ

Due to Lemma 2.1(ii)–(iii), this definition does not depend on how the sequences are extended. 2.8. Collocation IGABEM In this section, we show how to choose the collocation points xj for j ¼ 1; …; Ncol in (1.3). First, we note that Lemma 2.1(i) implies that   1 Ri;p j ½a;bÞ : i ¼ 1  p; …; N  #b þ 1 ○γ j ½a;bÞ ð2:16Þ

ð2:21Þ

of the related NURBS spaces. Since the weights in W H are just convex combinations of the weights in W 0 , it holds min W 0 r min W H r max W H r max W 0 . For details, we refer to [15, Section 4.2]. Then, the adaptive algorithm reads as follows: Algorithm 2.2. Input: Adaptivity parameter 0 o θ r 1, polynomial 0 ¼K  h , initial order p A N0 , initial partition T 0 ¼ T h with knots K weights W 0 ¼ W h . Adaptive loop: Iterate the following steps (i)–(vi), until ϱh is sufficiently small: p

c ðK  h ; W h Þ from Galerkin BEM (i) Compute approximation ϕh A N (1.2) resp. collocation BEM (1.3). (ii) Compute indicators ϱh ðzÞ for all nodes z A N h . (iii) Determine a set Mh D N h of minimal cardinality such that X θ ϱ2h r ϱh ðzÞ2 : ð2:22Þ z A Mh

(iv) If both nodes of an element T A T h belong to Mh , T will be marked. (v) For all other nodes in Mh , the multiplicity will be increased if it is smaller than p þ 1, otherwise the elements which contain one of these nodes z A Mh will be marked. (vi) Refine all marked elements T A T h by bisection (insertion of a node with multiplicity one) of the corresponding T A T h . Use further bisections to guarantee that the new partition T H satisfies

κ ðT H Þ r2κ ðT 0 Þ:

ð2:23Þ

Update h↦H, i.e., replace T h by T H .

for Γ ¼ ∂Ω resp.   Ri;p j ½a;b : i ¼ 1  p; …; N  p ○γ  1

ð2:17Þ p

c ðK  h ; W h Þ. Recall #b ¼ p þ 1 for for Γ ⫋ ∂Ω forms a basis of N Γ ⫋ ∂Ω. For simplicity, suppose #b ¼ p þ 1 also for Γ ¼ ∂Ω. This gives N col ¼ N:

considered ηh for Galerkin IGABEM, the current focus is on μh and collocation IGABEM. Suppose that Γ is represented by a NURBS curve of degree p A N0 . This induces the initial partition T 0 of Γ with nodes N 0 , related nodes N 0 in the parameter domain, and positive weights W 0 . Each node has a multiplicity lower than or equal to p þ 1, where for Γ ⫋ ∂Ω or collocation IGABEM we suppose #a ¼ #b ¼ p þ 1. For Γ ¼ ∂Ω, we suppose hT r j Γ j =4 for all T A T 0 . As the initial trial space, we consider

ð2:18Þ

For j ¼ 1; …; N, the collocation point xj is defined through the arithmetic mean of p þ 2 knots in the parameter domain Pj t k ¼ jp1 k   : ð2:19Þ xj ¼ γ ðx j Þ with x j ≔ pþ2 2.9. Adaptive algorithm

Output: Adaptively generated partition T h with corresponding solution ϕh and error estimator ϱh.□ Remark 2.3. (i) While θ ¼ 1 leads essentially to uniform refinement, θ⪡1 leads to highly adapted partitions. Note that the smaller the θ, the more iterations of the adaptive loop are required. In our experiments below, θ ¼0.75 appeared to be an appropriate compromise which led to optimal convergence behavior. (ii) The estimate (2.23) in step (iv) of the adaptive algorithm can be achieved by some extended 1D bisection algorithm from [1]. The latter guarantees that the overall number of elements is bounded by the sum of elements in the initial partition plus the number of marked elements.□ 3. Numerical experiments

Finally, we recall an adaptive algorithm from our preceding work [15], which steers the h-refinement of the partition T h as well as the increase of the multiplicity of the nodes N h . While [15]

In this section, we empirically investigate the performance of Algorithm 2.2 for Galerkin as well as collocation IGABEM in three

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

typical situations: in Section 3.2, the boundary Γ ¼ ∂Ω is closed and the solution exhibits a generic (i.e., geometry induced) singularity. In Section 3.3, the solution is piecewise smooth on Γ ¼ ∂Ω with certain jumps which locally require discontinuous ansatz functions. In Section 3.4, we consider a slit problem. In all examples, the exact solution is known. This allows to analyze the reliability and efficiency of the proposed estimators. The boundary part Γ is parametrized by a NURBS curve γ, i.e., the parametrization has the special form X  γ ðtÞ ¼ C i RKi;pγ ;W γ ðtÞ ð3:1Þ

145

0.1

γ(1/6)

0.05

γ(1/3)

γ(1/2)

0

iAZ

 γ and W γ are for all t A ½a; b. Here, p A N is the polynomial degree, K knots and weights as in Section 2.9 and ðC i Þi A Z are control points in R2 which are periodic for closed Γ ¼ ∂Ω. We choose the same polynomial degree p for our ansatz spaces c p ðK  h ; W h Þ. For the initial knots and weights, we choose Xh ¼ N  γ and W h ¼ W γ . As the ansatz spaces are nested, it always h ¼K K holds

γ(1)

γ(2/3) −0.05

γ(5/6) −0.1

γ 1 ; γ 2 A N p ðK h ; W h Þj ½a;b ;

ð3:2Þ

where γ 1 ; γ 2 denote the first resp. second component of γ. Therefore, this approach reflects the main idea of isogeometric analysis, i.e., the same space is used for the geometry and for the approximation. For adaptive Galerkin IGABEM as well as adaptive collocation IGABEM, we compare uniform refinement, where Mh ¼ N h and hence all elements are refined, and adaptive refinement with θ ¼ 0.75. In addition, we also consider discontinuous piecewise polynomials. Note that this is formally only a special case if wj ¼ 1 for all weights wj of W h and #zj ¼ p þ 1 for all nodes zj A N h . As basis for the considered ansatz spaces, we use (2.16) resp. (2.17). To calculate the Galerkin matrix, the collocation matrix, the Faermann error estimator, and the weighted-residual error estimator, we transform the weakly singular integrands into a sum of a smooth part and a logarithmically singular part. Then, we use adapted Gauss quadrature to compute the resulting integrals with appropriate accuracy; see [16, Section 5] for details. For the weighted-residual error estimator (1.5), we replace j ωh ðzÞj by the length j γ  1 ðωh ðzÞÞj , since this eases the calculation. Note that j ωh ðzÞj C j γ  1 ðωh ðzÞÞj , where the hidden constants depend only on the parametrization γ. gal To calculate, the exact error, we proceed as follows: Let ϕh A gal X h be the Galerkin approximation with c h the corresponding coefficient vector. Let ϕh A X h be the collocation approximation the corresponding coefficient vector. Let V gal be the with ccol h h Galerkin matrix of the h-th step. With the Galerkin orthogonality and the energy norm j j j ϕ j j j 2 ¼ 〈V ϕ; ϕ〉, obtained by Aitken's Δ2extrapolation, we can compute the energy error as col

jjj ϕ  ϕh jjj 2 ¼ jjj ϕ jjj 2  jjj ϕh jjj 2 ¼ jjj ϕ jjj 2  〈V gal cgal ; c gal 〉; h h h gal

gal

ð3:3Þ resp. jjj ϕ  ϕh jjj 2 ¼ jjj ϕ  ϕh jjj 2 þ jjj ϕh  ϕh jjj 2 col

gal

¼

jjj ϕ  ϕ

gal h

gal

jjj

2

col

gal col þ〈V gal ðcgal  ccol h Þ; ðc h c h Þ〉: h h

ð3:4Þ

3.1. Laplace–Dirichlet problem In the first two examples, we consider the Laplace–Dirichlet problem  Δu ¼ 0 in Ω

and

u ¼ g on Γ

ð3:5Þ

−0.1

−0.05

0

0.05

0.1

Fig. 1. Geometry and initial nodes for the experiment from Section 3.2.

for given Dirichlet data g A H 1=2 ðΓ Þ and closed boundary Γ ¼ ∂Ω. The problem is equivalent to the integral equation (1.1) with f ¼ ðK þ σ Þg, i.e. V ϕ ¼ ðK þ σ Þg where KgðxÞ ≔ 

on Γ ;

Z 1 gðyÞ∂νðyÞ log j x  yj dy 2π Γ

ð3:6Þ

ð3:7Þ

denotes the double-layer integral operator and σ ðxÞ ¼ 1=2 for all x A Γ except for the corners, where σ ðxÞ ¼ α=ð2π Þ with the corresponding interior angle α. The unique solution of (1.1) is the normal derivative ϕ ¼ ∂u=∂ν of the solution u A H 1 ðΩÞ of (3.5). For more details, see e.g. [28, Sections 6.3 and 6.6]. 3.2. Problem with generic singularity As first example, we consider the Laplace–Dirichlet problem (3.5) on the pacman geometry   π π  1 Ω ≔ rð cos ðβÞ; sin ðβÞÞ : 0 r r o ; β A  ; ; 10 2τ 2τ with τ ¼4/7; see Fig. 1. The geometry is parametrized on ½0; 1 by a NURBS curve of degree p ¼2. We prescribe the exact solution of (3.5) as  uðx; yÞ ¼ r τ cos τβ in polar coordinates ðx; yÞ ¼ rð cos β; sin βÞ. We consider the corresponding integral equation (3.6). The normal derivative ϕ ¼ ∂u=∂ν of u reads   ! cos ðβ Þ cos τβ þ sin ðβ Þ sin τβ    νðx; yÞ  τ  r τ  1 ϕðx; yÞ ¼ sin ðβ Þ cos τβ  cos ðβ Þ sin τβ and has a generic singularity at the origin. In Fig. 2, the solution ϕ is plotted over the parameter domain. The singularity is located at t ¼ 1=2 and two jumps are located at t ¼ 1=3 resp. t ¼ 2=3. In Fig. 3, error and error estimators are plotted. All values are plotted in a double logarithmic scale such that the experimental convergence rates are visible as the slope of the corresponding curves. Since the solution lacks regularity, uniform refinement leads to the suboptimal rate OðN  4=7 Þ for the energy error, whereas

146

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

Fig. 2. Experiment with singular solution on pacman geometry from Section 3.2. The singular solution ϕ○γ is plotted on the parameter interval, where 0.5 corresponds to the origin, where ϕ is singular.

Fig. 4. Experiment with singular solution on pacman geometry from Section 3.2. The plot shows the efficiency indices jjj ϕ ρhϕ jjj for the estimators ρh A fηh ; μh g, h where adaptivity is driven by ρh.

Fig. 3. Experiment with singular solution on pacman geometry from Section 3.2. Error and estimator are plotted versus the number of knots N. Uniform, ηh-driven and μh-driven refinement is considered.

Fig. 5. Experiment with singular solution on pacman geometry from Section 3.2. The errors from all presented adaptive IGABEM strategies are plotted versus the number of knots N.

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

0.5

147

(1)

(3/4)

(1/4)

(1/2)

0.4

0.3

0.2

0.1

0

Fig. 6. Experiment with singular solution on pacman geometry from Section 3.2. Histogram of number of knots over the parameter domain. Knots with maximal multiplicity p þ 1 ¼ 3 are marked.

0

0.1

0.2

0.3

0.4

0.5

Fig. 8. Geometry and initial nodes for the experiments from Section 3.3.

Fig. 7. Experiment with singular solution on pacman geometry from Section 3.2. The errors from uniform/adaptive BEM with discontinuous piecewise polynomials and uniform/adaptive IGABEM are plotted versus the number of knots N.

adaptive refinement leads to the optimal rate OðN  7=2 Þ. In each case, the curves for the two different estimators ηh and μh and the error are parallel. In Fig. 4, we plot the ratios ηh = jjj ϕ  ϕh jjj resp. μh = jjj ϕ  ϕh jjj . Throughout, these ratios stay between 0.5 and 2.7 which underlines an accurate error estimation for both error estimators. Fig. 5 shows the errors of all considered adaptive IGABEM strategies. We observe a very similar behavior. For adaptive refinement, Fig. 6 provides a histogram of the knots in ½a; b of the last refinement step for collocation IGABEM with ρh ¼ μh , for the other adaptive strategies, the output looks similar (not displayed). We see that the algorithm mainly refines the mesh around the singularity at t ¼ 1=2. Additionally, the multiplicity at the jump points t ¼ 1=3 and t ¼ 2=3 appears to be maximal so that the discrete solution ϕh also mimics the discontinuities of the exact solution ϕ. In Fig. 7, we finally compare standard BEM with discontinuous piecewise polynomials against IGABEM. For the error estimation we use the weighted-residual estimator μh. The output looks similar if ηh is used instead (not displayed). All approaches show

Fig. 9. Experiment with jump solution on square from Section 3.3. The solution ϕ○γ is plotted on the parameter interval.

similar convergence rates, however we clearly observe better multiplicative constants for Galerkin IGABEM and collocation IGABEM than for standard BEM. 3.3. Adaptive IGABEM for problem with jump solution As second example, we consider the Laplace–Dirichlet problem (3.5) on the square Ω ¼ ½0; 1=22 ; see Fig. 8. The geometry is parametrized on ½0; 1 by a NURBS curve of degree p ¼1. We prescribe the exact solution of (3.5) as uðx; yÞ ¼ sinhð2π xÞ cos ð2π yÞ: We consider the corresponding integral equation (3.6). The normal derivative ϕ ¼ ∂u=∂ν of u reads ! coshð2π xÞ cos ð2π yÞ ϕðx; yÞ ¼ 2π  νðx; yÞ: sinhð2π xÞ cos ð2π yÞ It is smooth up to four jumps as can be seen in Fig. 9.

148

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

Fig. 10. Experiment with jump solution on square from Section 3.3. Error and estimator are plotted versus the number of knots N. Uniform, ηh-driven and μhdriven refinement is considered.

In Fig. 10 we plot error and error estimators. The solution ϕ○γ has jumps at the points t ¼ 1=4, t ¼ 1=2, t ¼ 3=4 and t ¼1 resp. t ¼0.  γ used for the parametrization of Γ all have mulAs the knots K tiplicity one, the functions of the isogeometric start approximation space are continuous at the points t ¼ 1=4, t ¼ 1=2 and t ¼ 3=4. Uniform refinement, where only h-refinement takes place, leads to the suboptimal rate OðN  1 Þ for the energy error, whereas adaptive refinement increases the knot multiplicity at these problematic points and leads again to the optimal rate OðN  5=2 Þ. In Fig. 11, we plot the efficiency indices ηh = jjj ϕ  ϕh jjj resp. μh = jjj ϕ  ϕh jjj . Throughout, these ratios stay between 0.1 and 2.2. Fig. 12 shows the errors of all considered adaptive IGABEM strategies. We observe that ηh leads to slightly better results than μh, while there appears to be almost no difference between Galerkin IGABEM and collocation IGABEM. In Fig. 13, standard BEM with discontinuous piecewise polynomials is compared against IGABEM. For adaptivity, we use the weighted-residual estimator μh. The output looks similar if the estimator ηh is used (not displayed). We observe that in this example uniform standard BEM is superior to uniform IGABEM. This is of course due to the fact that standard BEM uses ansatz

Fig. 11. Experiment with jump solution on square from Section 3.3. The plot shows the efficiency indices jjj ϕ ρhϕ jjj for the estimators ρh A fηh ; μh g, where adaptivity is h driven by ρh.

Fig. 12. Experiment with jump solution on square from Section 3.3. The errors from all presented adaptive IGABEM strategies are plotted versus the number of knots N.

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

149

Fig. 13. Experiment with jump solution on square from Section 3.3. The errors from uniform BEM with discontinuous piecewise polynomials and uniform/adaptive IGABEM are plotted versus the number of knots N.

1.5

1

0.5

0

γ(0)

γ(1/5)

γ(2/5)

γ(3/5)

γ(1)

γ(4/5)

−0.5

−1

−1.5 −1.5

Fig. 15. Experiment with singular solution on slit from Section 3.4. Error and estimator are plotted versus the number of knots N. Uniform, ηh-driven and μhdriven refinement is considered. −1

−0.5

0

0.5

1

1.5

Fig. 14. Geometry and initial nodes for the experiment from Section 3.4.

spaces which are discontinuous at the jumps of ϕ. However, with the use of adaptive multiplicity increase this is fixed as can be seen in the convergence plot, where we again see that adaptive IGABEM leads to better results than adaptive standard BEM. It is also interesting that adaptive standard BEM converges with a better multiplicative constant than uniform standard BEM. This is due to the fact that the solution is zero on ½1=4; 1=2 and ½3=4; 1, wherefore the adaptive algorithm uses only few elements in this area.

3.4. Adaptive IGABEM for slit problem As last example, we consider a crack problem on the slit

Γ ¼ ½  1; 1  f0g. We parametrize Γ by a NURBS curve of degree p ¼1. For f ðx; 0Þ ≔  x=2 and the single-layer operator V, the exact

Fig. 16. Experiment with singular solution on slit from Section 3.4. The errors from all presented adaptive IGABEM strategies are plotted versus the number of knots N.

150

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

Fig. 18. Experiment with singular solution on slit from Section 3.4. The errors from uniform BEM with discontinuous piecewise polynomials and uniform/adaptive IGABEM are plotted versus the number of knots N.

4. A posteriori error estimation for IGABEM 4.1. Main results For T A T h , we inductively define the patch ωm h ðTÞ D Γ of order m A N0 by

ω0h ðTÞ ≔ T; þ1 ωm ðTÞ ≔ ⋃fT 0 A T h : T 0 \ ωm h ðTÞ a ∅g h

ð4:1Þ

Theorem 4.2 requires the following two assumptions on T h and X h for some fixed integer m A N0 : (A1) For each T A T h , there exists some fixed function ψ T A X h with connected support suppðψ T Þ such that T D suppðψ T Þ D ωm h ðTÞ: Fig. 17. Experiment with singular solution on slit from Section 3.4. The plot shows the efficiency indices j j j ϕ ρhϕ j j j for the estimators ρh A fηh ; μh g, where adaptivity is h driven by ρh.

solution of (1.1) reads x

ϕðx; 0Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffi2ffi: 1x

 ε

Note that ϕ A H ðΓ Þ⧹L2 ðΓ Þ for all ε 4 0 and that ϕ has singularities at the tips x ¼ 7 1 (Fig. 14). In Fig. 15, error and error estimators for the uniform and for the adaptive approach are plotted. The error is obtained via (3.3) resp. (3.4), where jjj ϕ jjj 2 ¼ π =4 is computed analytically. Since the solution lacks regularity, uniform refinement leads to the suboptimal rate OðN  1=2 Þ, whereas adaptive refinement leads to the optimal rate OðN  5=2 Þ. The curves for the two estimators and the error are again parallel. In Fig. 17, we plot the efficiency indices ηh = jjj ϕ  ϕh jjj resp. μh = jjj ϕ  ϕh jjj . Fig. 16 shows the errors of all considered adaptive IGABEM strategies. Here, ηh-adaptive Galerkin IGABEM and μh-adaptive collocation IGABEM lead to the best results. In Fig. 18 we compare standard BEM against IGABEM, where we use ρh ¼ μh . While adaptive Galerkin IGABEM and adaptive standard BEM lead to optimal convergence rates, the best results are achieved with adaptive collocation IGABEM.

ð4:2Þ

(A2) There exists some constant q A ð0; 1 such that J 1  ψ T J 2L2 ðsuppðψ

T ÞÞ

r ð1  qÞj suppðψ T Þj

ð4:3Þ

for all T A T h . The first theorem shows that these assumptions are, in particular, satisfied for NURBS spaces. Theorem 4.1 (Feischl et al. [15, Theorem 4.4]). For p A N0 and c p ðK  h ; W h Þ satisfies the assumptions m ≔ ⌈p=2⌉, the space X h ≔ N (A1)–(A2). The constant 0 oq r1 depends only on κ ðT h Þ, minðW h Þ, maxðW h Þ, p, and γ.□ The main result of [15] reads as follows: Theorem 4.2 (Feischl et al. [15, Theorem 3.1]). For any approximation ϕh A L2 ðΓ Þ, the residual r h ¼ f  V ϕh satisfies the efficiency estimate 0 11=2 X 2A @ ηh ≔ ηh ðzÞ r C eff jjj ϕ  ϕh jjj ð4:4Þ zANh

with ηh ðzÞ ≔ j r h j H1=2 ðωh ðzÞÞ . If the mesh T h and the discrete space X h satisfy assumptions (A1)–(A2), the Galerkin solution ϕh A X h of (1.2) also satisfies the reliability estimate jjj ϕ  ϕh jjj r C rel ηh :

ð4:5Þ

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

Z Z

The constant C eff 4 0 depends only on Γ , while C rel 4 0 additionally depends on m, κ ðT h Þ, and q.□

r

where C Γ 4 0 is the constant from (2.7). If collocation IGABEM as in Section 2.8 is used, the patch ωp þ 1 ðTÞ contains a collocation point and therefore a root of the residual rh, for each T A T h . Hence, the condition of the following theorem is fulfilled with m ¼ p þ 1. Theorem 4.4. Suppose that wither ϕh A X h is the Galerkin solution of (1.2), where X h satisfies (A1)–(A2), or that the residual r h ¼ f  V ϕh has at least one root in each ωm h ðTÞ for all T A T h and some fixed m A N0 . Then, 0 11=2 X 1 jjj ϕ  ϕ jjj r μ ≔ @ μ ðzÞA ð4:7Þ C h

rel

h

h

zANh

The constant C rel 40 depends with μh ðzÞ ≔ j ωh ðzÞj only on Γ, m, κ ðT h Þ, and, in the first case, q. 1=2

J r 0h J L2 ðωh ðzÞÞ .

Lemma 4.5. For any connected ω D Γ , whose length satisfies j ω j r 34L if Γ ¼ ∂Ω, there holds for all u A H 1 ðΓ Þ:

ð4:8Þ 1

Proof. We recall that for a finite interval I
R, H ðIÞ coincides with the space of all absolutely continuous functions on I with L2 derivative; see, e.g., [11, p. 306]. Step 1: First we consider I ¼ ð0; 1Þ and prove ð4:9Þ

We use the transformation theorem, with r ¼ ρðs  tÞ þt and s  t ¼ σ , as well as the Cauchy Schwarz inequality to get  Z Z  uðsÞ  uðtÞ2   j uj 2H1=2 ðIÞ ¼  s  t  ds dt I I Z Z R u0 ðrÞ dr  R u0 ðrÞ dr 2  ð0;sÞ  ð0;tÞ ¼   ds dt  st I I  2 Z Z Z     u0 ρðs tÞ þ t dρ ds dt ¼   I I Z Z ZI  j u0 ρðs  tÞ þ t j 2 dρ ds dt r Z ZI ZI I  j u0 ρσ þ t j 2 dρ dσ dt: ¼ ð  t;1  tÞ I

I

We formally extend u0 by zero to R. This and the Fubini theorem lead to Z Z Z  j uj 2H1=2 ðIÞ r j u0 ρσ þt j 2 dρ dσ dt I

ð  1;1Þ I

I

ð  1;1Þ

J u0 J 2L2 ðRÞ dσ dρ ¼ 2j uj 2H1 ðIÞ :

Step 2: If I DR is an arbitrary finite interval, it holds ð4:10Þ

Without loss of generality, let I ¼ ðc; dÞ be open. We define the function uð0;1Þ : ð0; 1Þ-R : τ↦uðτðd  cÞ þ cÞ. Obviously, it holds uð0;1Þ A H 1 ð0; 1Þ with u0ð0;1Þ ðτÞ ¼ ðd  cÞu0 ðτðd  cÞ þ cÞ. The transformation theorem with s ¼ σ ðd  cÞ þ c, t ¼ τðd  cÞ þ c, and r ¼ ρðd  cÞ þc, and (4.9) yield  Z Z  uðsÞ uðtÞ2   j uj 2H1=2 ðIÞ ¼  s t  ds dt I I   Z Z uðσ ðd  cÞ þ cÞ uðτðd  cÞ þ cÞ2   ds dt ¼   σ τ ð0;1Þ ð0;1Þ ¼ j uð0;1Þ j 2H1=2 ð0;1Þ r 2j uð0;1Þ j 2H1=2 ð0;1Þ Z Z ¼2 j u0ð0;1Þ ðρÞj 2 dρ ¼ 2ðd  cÞ j u0 ðrÞj 2 dr ¼ 2j Ij j uj H1 ðIÞ : ð0;1Þ

I

Step 3: We show (4.8). Let I be a real interval with γ L ðIÞ ¼ ω. Then, (2.8) and (4.10) give j uj H1=2 ðωÞ ¼ j uj 2H1=2 ðγ

L ðIÞÞ

r C 2Γ j u○γ L j 2H1=2 ðIÞ r 2C 2Γ j ω j j u○γ L j 2H1 ðIÞ

¼ 2C 2Γ j ω j J u0 J 2L2 ðωÞ :

We only need the following lemma, whose proof is inspired by [22, Proposition 2.2], where an analogous assertion for norms instead of seminorms is found. The assertion itself is also stated in [4, Lemma 7.4] in a more general way. Indeed a similar version of (4.8) holds even for the Hs-seminorm, 0 o s o 1. However, in [4], the proof is only given for the hardest case 1=2 o s o 1.

j uj 2H1=2 ðIÞ r 2j uj 2H1 ðIÞ :

ð  1;1Þ R

 j u0 ρσ þt j 2 dt dσ dρ

j uj 2H1=2 ðIÞ r2j Ij j uj 2H1 ðIÞ :

4.2. Proof of Theorem 4.3

j uj 2H1=2 ðωÞ r 2 C 2Γ j ω j J u0 J 2L2 ðωÞ

I

Z Z ¼

The following two theorems are the mathematical contributions of this work to the field of IGABEM. They apply to both, Galerkin IGABEM and collocation IGABEM. Theorem 4.3. For any approximation ϕh A L2 ðΓ Þ and r h ≔ f  V ϕh , the indicator ηh ðzÞ ≔ j r h j H1=2 ðωh ðzÞÞ is bounded above by the weighted-residual indicator μh ðzÞ ≔ j ωh ðzÞj 1=2 J r 0h J L2 ðωh ðzÞÞ pffiffiffi ηh ðzÞ r 2 C Γ μh ðzÞ; ð4:6Þ

Z

151

This concludes the proof.□ 4.3. Proof of Theorem 4.4 We use the following estimate from [12, Lemma 2.3]; see [16, Proposition 2.13] for a detailed proof. Lemma 4.6. There exists a constant C 1 4 0 such that, for all u A H 1=2 ðΓ Þ, it holds X 1 X J u J 2H1=2 ðΓ Þ r j uj 2H1=2 ðω ðzÞÞ þ C 1 hT J u J 2L2 ðTÞ : ð4:11Þ h

zANh

T AT h

Γ and κ ðT h Þ.□

The constant only depends on

Proof of Theorem 4.4. If the residual is orthogonal to some X h satisfying (A1)–(A2), the assertion follows at once from Theorem 4.2 in combination with Eq. (4.6). If the residual has local roots, we first note that J ϕ  ϕh J

  1=2

H

ðΓ Þ

C J f  V ϕh J H1=2 ðΓ Þ ¼ J r h J H1=2 ðΓ Þ ;

ð4:12Þ

since V is an isomorphism. The hidden constants only depend on Γ. Taking u ¼ r h in Lemma 4.6, it only remains to estimate the sum P 1 2 m T A T h hT J r h J L2ðTÞ . Note that shape regularity yields j ωh ðTÞj r ð2m þ 1Þκ ðT h Þm hT . Replacing T by ωm ðTÞ, we apply Friedrich's inequality to see X T AT h

1

hT J r h J 2L2 ðTÞ r

X T AT h

1

hT J r h J 2L2 ðωm ðTÞÞ r h

r ð2m þ 1Þ κ ðT h Þ 2

2m

X T AT h

r ð2m þ 1Þ κ ðT h Þ 3

3m

X

T AT h

r ð2m þ 1Þ κ ðT h Þ 3

3m

X

zANh

This concludes the proof.□

X j ωm ðTÞj 2 h J r 0h J 2L2 ðωm ðTÞÞ hT T AT h

hT J r 0h J 2L2 ðωm ðTÞÞ h

hT J r 0h J 2L2 ðTÞ j ωh ðzÞj J r 0h J 2L2 ðωðzÞÞ :

152

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153

5. Conclusion

observed that

ηh C jjj ϕ  ϕh jjj C μh ;

5.1. Analytical results In this work, we considered adaptive BEM for weakly singular integral equations V ϕ ¼ f associated to elliptic PDEs in 2D. As model example served the 2D Laplacian, but the results apply as long as   1=2

ðΓ Þ-H 1=2 ðΓ Þ is an elliptic isomorphism. With the resiV:H dual r h ≔ f  V ϕh , we transferred the weighted-residual error estimator

μh ¼ J h1=2 r0h J L2 ðΓÞ

ð5:2Þ

see Theorem 4.4. In our preceding work [15], we considered the residual error estimator 0 11=2 Z 2 XZ j r ðxÞ  r ðyÞj h h ηh ¼ @ dy dxA ð5:3Þ j x  yj 2 ωh ðzÞ ωh ðzÞ zAN h

proposed in [12]. In [15], we transferred this estimator from standard BEM with piecewise polynomials to IGABEM. Independently of the discretization, we proved the general efficiency estimate

ηh r C eff jjj ϕ  ϕh jjj

ð5:4Þ

while our proof of the converse estimate jjj ϕ  ϕh jjj r C rel ηh is restricted to Galerkin IGABEM. However, the combination of (5.2) and (5.4) provides also full error control 1 ηh r jjj ϕ  ϕh jjj r C rel μh C eff

ð5:5Þ

for collocation IGABEM computations in 2D. Moreover, this estimate implies the global relation ηh ≲ μh , and we even proved

μh ðzÞ r C loc ηh ðzÞ for all z A N h

i.e., both error estimators are efficient and reliable. The efficiency indices ηh = jjj ϕ  ϕh jjj and μh = jjj ϕ  ϕh jjj appeared to be r 3, i.e., the overestimation of the energy error is very moderate. We note that only the equivalence ηh C jjj ϕ  ϕh jjj for Galerkin IGABEM as well as the bounds jjj ϕ  ϕh jjj ≲ μh and ηh ≲ jjj ϕ  ϕh jjj have thoroughly been proved mathematically. 5.3. Open questions and future work

ð5:1Þ

proposed in [9,3] from standard BEM with lowest-order polynomials to IGABEM, where we considered the Galerkin method as well as collocation. For either discretization, we proved that μh is reliable jjj ϕ  ϕh jjj r C rel μh ;

ð5:7Þ

ð5:6Þ

for the respective nodal contributions defined in (1.4) resp. (1.5); see Theorem 4.3 which holds independently of the discretization employed. 5.2. Numerical results We proposed an adaptive algorithm which is capable to steer the mesh-refinement as well as the knot multiplicity in Galerkin and collocation IGABEM computations; see Algorithm 2.2. Numerical experiments in Section 3 underline that generic singularities of the (unknown) exact solutions lead to reduced experimental convergence behavior if the underlying mesh is not appropriately graded. This is a well-known fact for standard BEM with piecewise polynomials, but also applies to IGABEM. Consequently, the gain of adaptive IGABEM (resp. the loss in case of uniform meshes) is huge due to the higherorder ansatz functions of IGABEM, and therefore adaptivity seems to be a must to exploit the full potential of isogeometric analysis. In several numerical experiments, we showed that the proposed algorithm is capable to recover the optimal order of convergence. The gain of IGABEM is that the algorithm chooses smooth NURBS, where the exact solution appears to be smooth, while discontinuities and singularities are well detected and appropriately resolved. Compared to standard BEM with discontinuous piecewise polynomials, this leads to a smaller number of degrees of freedom for comparable accuracies. For collocation IGABEM as well as Galerkin IGABEM and independently of the (uniform or adaptive) mesh-refinement, we

All considered numerical experiments show optimal convergence of the estimator and the error. Understanding this observation mathematically in the spirit of [5] is one of our goals for future research. However, it is questionable if an analogous version of the reduction property on refined element domains [5, (A2)] can be proved for the Faermann estimator ηh. Indeed, this is yet an open problem even for standard BEM with piecewise polynomials; see [14], where at least convergence of an h-adaptive algorithm with ηh is analyzed. For the weighted-residual error estimator μh the axioms of [5] are satisfied for standard Galerkin BEM with piecewise polynomials, see [5, Section 5.4]. For collocation IGABEM there remain two challenging mathematical questions: First, one needs further investigation on the unique solvability of the discrete system. Second, the quasi-orthogonality [5, (A3)] is unclear for collocation methods. As mentioned, we observed in all numerical experiments reliability as well as efficiency of the used error estimators. However, it remains to mathematically verify the reliability estimate jjj ϕ  ϕh jjj ≲ ηh for collocation BEM and the efficiency estimate μh ≲ jjj ϕ  ϕh jjj þosc, at least for some higher-order oscillation terms osc. Again, these estimates are yet open problems even for standard BEM. For lowest-order Galerkin BEM, the efficiency estimate is proved in [1, Theorem 4] under additional regularity assumptions on the Dirichlet data g in (3.5). Finally, the ultimate goal is of course to analyze and apply the estimators ηh and μh in 3D Galerkin IGABEM. For 3D one has to consider, e.g., T-splines [27] or hierarchical B-splines [2], because, in contrast to multivariate NURBS, they naturally allow for local mesh refinement. Ref. [13] shows that ηh is reliable and efficient for standard BEM with piecewise polynomials, whereas [7] proves reliability for μh. In [5, Section 5.4] optimal convergence of adaptive h-refinement for μh is proved. The estimate ηh ≲ μh as well as plain convergence for ηhbased adaptivity is analyzed in [14]. The transfer of the mentioned results from standard BEM to adaptive IGABEM leaves interesting and challenging questions for future research.

Acknowledgment The authors acknowledge support through the Austrian Science Fund (FWF) under Grant P27005 Optimal adaptivity for BEM and FEM–BEM coupling. In addition, M.F., G.G. and D.P. are supported through the FWF doctoral school Nonlinear PDEs funded under Grant W1245.

References [1] Aurada Markus, Feischl Michael, Führer Thomas, Karkulik Michael, Praetorius Dirk. Efficiency and optimality of some weighted-residual error estimator for adaptive 2D boundary element methods. Comput Methods Appl Math 2013;13 (3):305–32. [2] Buffa Annalisa, Giannelli Carlotta. Adaptive isogeometric methods with hierarchical splines: error estimator and convergence. ArXiv preprint arXiv:1502. 00565, 2015.

M. Feischl et al. / Engineering Analysis with Boundary Elements 62 (2016) 141–153 [3] Carstensen Carsten. An a posteriori error estimate for a first-kind integral equation. Math Comput 1997;66(217):139–55. [4] Carstensen Carsten, Faermann Birgit. Mathematical foundation of a posteriori error estimates and adaptive mesh-refining algorithms for boundary integral equations of the first kind. Eng Anal Bound Elem 2001;25(7):497–509. [5] Carstensen Carsten, Feischl Michael, Page Marcus, Praetorius Dirk. Axioms of adaptivity. Comput Math Appl 2014;67(6):1195–253. [6] Cottrell J Austin, Hughes Thomas JR, Bazilevs Yuri. Isogeometric analysis: toward integration of CAD and FEA. Hoboken, New Jersey: John Wiley & Sons; 2009. [7] Carstensen Carsten, Maischak Matthias, Stephan Ernst P. A posteriori error estimate and h-adaptive algorithm on surfaces for Symm's integral equation. Numer Math 2001;90(2):197–213. [8] Costabel Martin. Boundary integral operators on Lipschitz domains: elementary results. SIAM J Math Anal 1988;19(3):613–26. [9] Carstensen Carsten, Stephan Ernst P. Adaptive boundary element methods for some first kind integral equations. SIAM J Numer Anal 1996;33(6):2166–83. [10] de Boor Carl. B (asic)-spline basics. Mathematics Research Center. University of Wisconsin-Madison; 1986. [11] Evans Lawrence C. Partial differential equations. Graduate studies in mathematics, vol. 19, 2nd ed. Providence, RI: American Mathematical Society; 2010. [12] Faermann Birgit. Localization of the Aronszajn–Slobodeckij norm and application to adaptive boundary element methods: I. The two-dimensional case. IMA J Numer Anal 2000;20(2):203–34. [13] Faermann Birgit. Localization of the Aronszajn–Slobodeckij norm and application to adaptive boundary element methods. II. The three-dimensional case. Numer Math 2002;92(3):467–99. [14] Feischl Michael, Führer Thomas, Mitscha-Eibl Gregor, Praetorius Dirk, Stephan Ernst P. Convergence of adaptive BEM and adaptive FEM–BEM coupling for estimators without h-weighting factor. Comput Methods Appl Math 2014;14 (4):485–508. [15] Feischl Michael, Gantner Gregor, Praetorius Dirk. Reliable and efficient a posteriori error estimation for adaptive IGA boundary element methods for weakly-singular integral equations. Comput Methods Appl Mech Eng 2015;290:362–86. [16] Gantner Gregor. Isogeometric adaptive BEM [Master's thesis]. Vienna University of Technology; 2014. [17] Heltai Luca, Arroyo Marino, DeSimone Antonio. Nonsingular isogeometric boundary element method for Stokes flows in 3D. Comput Methods Appl Mech Eng 2014;268:514–39.

153

[18] Hughes Thomas JR, Cottrell J Austin, Bazilevs Yuri. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194(39–41):4135–95. [19] Hsiao George C, Wendland Wolfgang L. Boundary integral equations. Berlin: Springer; 2008. [20] McLean William. Strongly elliptic systems and boundary integral equations. Cambridge: Cambridge University Press; 2000. [21] Marussig Benjamin, Zechner Jürgen, Beer Gernot, Fries Thomas-Peter. Fast isogeometric boundary element method based on independent field approximation. Comput Methods Appl Math 2015;284:458–88. [22] Di Nezza Eleonora, Palatucci Giampiero, Valdinoci Enrico. Hitchhiker's guide to the fractional Sobolev spaces Bulletin des Sciences Mathématiques 2011;136(5):521-573. [23] Politis Costas, Ginnis Alexandros I, Kaklis Panagiotis D, Belibassakis Kostas, Feurer Christian. An isogeometric BEM for exterior potential-flow problems in the plane. In: 2009 SIAM/ACM joint conference on geometric and physical modeling. San Francisco: ACM; 2009. p. 349–54. [24] Peake Michael J, Trevelyan Jon, Coates Graham. Extended isogeometric boundary element method (XIBEM) for two-dimensional Helmholtz problems. Comput Methods Appl Mech Eng 2013;259:93–102. [25] Simpson Robert N, Bordas Stéphane PA, Trevelyan Jon, Rabczuk Timon. A twodimensional isogeometric boundary element method for elastostatic analysis. Comput Methods Appl Mech Eng 2012;209/212:87–100. [26] Sauter Stefan A, Schwab Christoph. Boundary element methods, Springer series in computational mathematics, vol. 39. Berlin: Springer-Verlag; 2011 [Translated and expanded from the 2004 German original]. [27] Scott Michael A, Simpson Robert N, Evans John A, Lipton Scott, Bordas Stephane PA, Hughes Thomas JR, et al. Isogeometric boundary element analysis using unstructured T-splines. Comput Methods Appl Mech Eng 2013;254:197– 221. [28] Steinbach Olaf. Numerical approximation methods for elliptic boundary value problems. New York: Springer; 2008 [Translated from the 2003 German original]. [29] Takahashi Toru, Matsumoto Toshiro. An application of fast multipole method to isogeometric boundary element method for Laplace equation in two dimensions. Eng Anal Bound Elem 2012;36(12):1766–75. [30] Zechner Jürgen, Marussig Benjamin, Beer Gernot, Fries Thomas-Peter. Isogeometric boundary element method with hierarchical matrices. ArXiv preprint arXiv:1406.2817, 2014.