Stable Implementation of Adaptive IGABEM in 2D in MATLAB

1.1 Isogeometric Analysis

The central idea of isogeometric analysis (IGA) is to use the same ansatz functions for the discretization of a partial differential equation (PDE) as is used for the representation of the corresponding problem geometry in computer-aided design (CAD), namely (rational) splines. This concept, originally invented in [33] for finite element methods (IGAFEM), has proved very fruitful in applications; see also the monograph [11].

1.2 Isogeometric Boundary Element Method

The idea of the boundary element method (BEM) is to reformulate the considered PDE on a domain Ω as an equivalent integral equation on the boundary ∂ ⁡ Ω . The solution of the latter is the missing Cauchy data, i.e., the Neumann or the Dirichlet data, which can be plugged into the so-called representation formula to obtain the PDE solution itself. Since CAD usually provides only a parametrization of ∂ ⁡ Ω , this makes BEM the most attractive numerical scheme if applicable (i.e., provided that the fundamental solution of the differential operator is explicitly known); see [40, 41] for the first works on isogeometric BEM (IGABEM) for 2D resp. 3D. Compared to standard BEM, which uses discontinuous or continuous piecewise polynomials as ansatz functions, IGABEM typically saves degrees of freedom by using smooth splines instead.

We refer to [45, 39, 44, 2, 38] for numerical experiments, to [47, 35, 15, 14, 16, 13] for fast IGABEM based on fast multipole ℋ-matrices resp. H 2 -matrices, and to [31, 6, 34, 1, 1, 17] for some quadrature analysis. For further references, we also mention the recent monograph [7].

1.3 Adaptive IGABEM

On the one hand, IGA naturally leads to high-order ansatz functions. On the other hand, however, optimal convergence behavior with higher-order discretizations is only observed in simulations if the (given) data as well as the (unknown) solution are smooth. Therefore, a posteriori error estimation and related adaptive strategies are mandatory to exploit the full potential of IGA. Rate-optimal adaptive strategies for IGAFEM (using hierarchical splines) have been proposed and analyzed independently in [10, 25] for IGAFEM, while the earlier work [9] proves only linear convergence. Meanwhile, these results have been extended to T-splines [27].

As far as IGABEM is concerned, available results focus on weakly singular integral equations with energy space H - 1 / 2 ⁢ ( Γ ) ; see [20, 18] for a posteriori error estimation, as well as [19] for the analysis of a rate-optimal adaptive IGABEM in 2D, and [24, 26, 28] and [8] for corresponding results for IGABEM in 3D with hierarchical splines and T-splines, respectively. Recently, [21] investigated optimal preconditioning for IGABEM in 2D with locally refined meshes. Moreover, [29] proved rate-optimality of an adaptive IGABEM in 2D for hyper-singular integral equations with energy space H 1 / 2 ⁢ ( Γ ) . The recent work [8] provides an overview of adaptive IGAFEM as well as IGABEM.

1.4 Model Problems

The present paper reports on our Matlab implementation of adaptive IGABEM in 2D, which is available online [30]. On a bounded Lipschitz domain Ω ⊂ R 2 with diam ⁢ ( Ω ) < 1 , we consider the Laplace model problem - Δ ⁢ P = 0 with Dirichlet/Neumann boundary conditions. Rewriting this PDE as boundary integral equation involves (up to) four integral operators, namely the weakly singular operator 𝔙, the hyper-singular operator 𝔚, the double-layer operator 𝔎, and the adjoint double-layer operator K ′ . Let Γ ⊆ ∂ ⁡ Ω be a connected part of the boundary. Given boundary densities ϕ , u : Γ → R and x ∈ Γ , the operators are formally defined by

where f : Γ → R is given and the integral density ϕ : Γ → R is sought. Moreover, we consider the hyper-singular integral equation

where 𝜈 denotes the outer unit normal vector on ∂ ⁡ Ω , g : Γ → R is given, and the integral density u : Γ → R is sought. For Γ = ∂ ⁡ Ω , these integral equations are equivalent formulations of the 2D Laplace problem - Δ ⁢ P = 0 in Ω, where (1.3) with f = ( 1 2 + K ) ⁢ ( P | Γ ) is equivalent to the Dirichlet problem and (1.4) with g = ( 1 2 - K ′ ) ⁢ ( ∂ ⁡ P ∂ ⁡ ν ) is equivalent to the Neumann problem. Knowing both Dirichlet as well as Neumann boundary data allows to compute 𝑃 via the representation formula P = V ~ ⁢ ( ∂ ⁡ P ∂ ⁡ ν ) - K ~ ⁢ ( P | Γ ) , where V ~ and K ~ are defined as in (1.1) and (1.2) for x ∈ Ω . While this approach is called direct, an alternative indirect approach is to make the ansatz P = V ~ ⁢ ϕ or P = K ⁢ u and solve V ~ ⁢ ϕ = f := P | Γ or W ⁢ u = g := ∂ ⁡ P ∂ ⁡ ν , respectively. In all cases, the singularities of the involved integral kernels pose a significant challenge, which is overcome in the presented stable IGABEM implementation.

1.5 Matlab Library IGABEM2D and Contributions

The library IGABEM2D, which provides an easily accessible Matlab implementation in the spirit of, e.g., [3, 5, 22], is available online [30]. It is distributed as a single file igabem2d.zip. To install the library, download and unpack the zip archive, start Matlab, and change to the root folder igabem2d/. You can directly run the main files IGABEMWeak.m and IGABEMHyp.m, where the adaptive algorithms, Algorithms 3.2 and 4.1, for the problems of Section 1.4 are implemented. There, you may choose out of a variety of different parameters. These functions both automatically run mexIGABEM.m, which adds all required paths, compiles the C files of source/ and stores the resulting mex-files in MEX-files/. The Matlab functions provided by IGABEM2D are contained in the folder functions/. Further details are also provided by README.txt, which is found in the root folder igabem2d/.

1.6 Outline

Section 2 recalls the definition of B-splines and NURBS and discusses necessary assumptions on parametrized NURBS geometries. Section 3 recalls the weakly singular integral equation along with the required Sobolev spaces on the boundary and formulates an adaptive Galerkin IGABEM, which is applied in numerical experiments at the end of the section. Similarly, Section 4 covers the hyper-singular integral equation. Details on the implementation, i.e., the stable computation of the involved (singular) boundary integrals, are given in Section 5. Finally, Appendix A provides an overview of all functions in IGABEM2D.

1.7 General Notation

Throughout and without any ambiguity, | ⋅ | denotes the absolute value of scalars, the Euclidean norm of vectors in R 2 , the measure of a set in ℝ (e.g., the length of an interval), or the arc length of a curve in R 2 . Throughout, mesh-related quantities have the same index, e.g., N ⋆ is the set of nodes of the partition Q ⋆ , and h ⋆ is the corresponding local mesh-width function etc. The analogous notation is used for partitions Q + , Q ℓ , etc. We sometimes use ⋅ ^ to transform notation on the boundary to the parameter domain or vice versa, e.g., Q ^ ⋆ is the partition of the parameter domain related to the partition Q ⋆ of the boundary.

2.1 B-Splines and NURBS in 1D

We consider an arbitrary but fixed sequence K ^ ⋆ := ( t ⋆ , i ) i ∈ Z on ℝ with multiplicities

which satisfy that t ⋆ , i - 1 ≤ t ⋆ , i for i ∈ Z and lim i → ± ∞ ⁡ t ⋆ , i = ± ∞ . Let N ^ ⋆ := { t ⋆ , i : i ∈ Z } = { x ^ ⋆ , j : j ∈ Z } denote the corresponding set of nodes with x ^ ⋆ , j - 1 < x ^ ⋆ , j for j ∈ Z . For i ∈ Z , the 𝑖-th B-spline [12] of degree p ∈ N 0 is defined for t ∈ R inductively by

where we suppose the convention ( ⋅ ) / 0 := 0 . It is well known that, for arbitrary I = [ a , b ) and p ∈ N 0 , the set { B ^ ⋆ , i , p | I : i ∈ Z ∧ B ^ ⋆ , i , p | I ≠ 0 } is a basis for the space of all right-continuous piecewise polynomials of degree lower than or equal to 𝑝 with breakpoints N ^ ⋆ on 𝐼 which are, at each knot t ⋆ , i , p - # ⋆ ⁢ t ⋆ , i times continuously differentiable if p - # ⋆ ⁢ t ⋆ , i ≥ 0 . Moreover, the B-splines are non-negative and have local support

Indeed, B ^ ⋆ , i , p depends only on the knots t ⋆ , i - 1 , … , t ⋆ , i + p . If t ⋆ , i - 1 < t ⋆ , i + p , then

Moreover, B-splines form a partition of unity, i.e., ∑ i ∈ Z B ⋆ , i , p = 1 for all p ∈ N 0 . For p ≥ 1 and i ∈ Z , the right derivative of the corresponding B-spline can be computed as

Finally, if ∑ i ∈ Z a ⋆ , i ⁢ B ^ ⋆ , i , p is a given spline and K ^ + is obtained from K ^ ⋆ by adding a single knot t ′ with t ′ ∈ ( t ⋆ , ℓ - 1 , t ⋆ , ℓ ] for some ℓ ∈ Z , there exist coefficients ( a + , i ) i ∈ Z such that

With the multiplicity # + ⁢ t ′ of t ′ in the knots K ^ + , the new coefficients can be chosen as convex combinations of the old coefficients

If # ⋆ ⁢ t k ≤ p + 1 for all k ∈ Z , these coefficients are unique. Proofs are found, e.g., in [12].

In addition to the knots K ^ ⋆ = ( t ⋆ , i ) i ∈ Z , we consider fixed positive weights W ⋆ := ( w ⋆ , i ) i ∈ Z with w ⋆ , i > 0 . For i ∈ Z and p ∈ N 0 , we define the 𝑖-th NURBS (non-uniform rational B-spline) by

Note that the denominator is locally finite and positive. For any p ∈ N 0 , we define the spline space as well as the rational spline space

2.2 Boundary Parametrization

Recall that Ω ⊂ R 2 is a bounded Lipschitz domain with diam ⁢ ( Ω ) < 1 and Γ ⊆ ∂ ⁡ Ω is a connected part of its boundary. We further assume that either Γ = ∂ ⁡ Ω is parametrized by a closed continuous and piecewise continuously differentiable path γ : [ a , b ] → Γ such that the restriction γ | [ a , b ) is even bijective, or that Γ ⫋ ∂ ⁡ Ω is parametrized by a bijective continuous and piecewise continuously differentiable path γ : [ a , b ] → Γ . In the first case, we speak of closed Γ = ∂ ⁡ Ω , whereas the second case is referred to as open Γ ⫋ ∂ ⁡ Ω . Throughout and by abuse of notation, we write γ - 1 for the inverse of γ | [ a , b ) and γ | ( a , b ] . The meaning will be clear from the context. For the left- and right-hand derivative of 𝛾, we assume that γ ′ ℓ ⁢ ( t ) ≠ 0 for t ∈ ( a , b ] and γ ′ r ⁢ ( t ) ≠ 0 for t ∈ [ a , b ) . Moreover, we assume that γ ′ ℓ ⁢ ( t ) + c ⁢ γ ′ r ⁢ ( t ) ≠ 0 for all c > 0 and t ∈ [ a , b ] and t ∈ ( a , b ) , respectively. Note that these assumptions provide a pointwise representation of the arc-length derivative

Finally and without loss of generality, we suppose that 𝛾 is positively oriented in the sense that the outer normal vector of Ω has the form

While the given assumptions are sufficient for the abstract analysis of adaptive BEM, the later implementation requires that 𝛾 is a NURBS curve of degree p γ ∈ N in the following sense. Let K ^ γ = ( t γ , i ) i = 1 N γ be a sequence of knots with

and multiplicity # γ ⁢ t γ , i ≤ p γ for i ∈ { 1 , … , N γ - p γ } . Let

denote the corresponding set of nodes with x ^ γ , j - 1 < x ^ γ , j for j ∈ { 2 , … , n γ } . Note that N γ = ∑ j = 1 n γ # γ ⁢ x ^ γ , j . With x ^ γ , 0 := a , we define the induced mesh on [ a , b ] ,

To use the definition of B-splines as in Section 2.1, we extend the knot sequence arbitrarily to ( t γ , i ) i ∈ Z with t γ , - p γ = ⋯ = t γ , 0 = a , t γ , i ≤ t γ , i + 1 , and lim i → ± ∞ ⁡ t γ , i = ± ∞ . For the extended sequence, we also write K ^ γ . Moreover, let W γ = ( w γ , i ) i = 1 - p γ N γ - p γ and C γ = ( C i ) i = 1 - p γ N γ - p γ be given positive weights and control points in R 2 , respectively, which satisfy that w γ , 1 - p γ = w γ , N γ - p γ and C γ , 1 - p γ = C γ , N γ - p γ in the case of Γ = ∂ ⁡ Ω . We extend W γ and C γ arbitrarily to ( w γ , i ) i ∈ Z and ( C γ , i ) i ∈ Z with positive weights w γ , i and control points C γ , i in R 2 , respectively. With the standard NURBS R ^ γ , i , p γ defined in Section 2.1, we suppose that 𝛾 has the form

where the second equality follows from the locality (2.1) of the B-splines resp. NURBS. The locality of NURBS even shows that this definition does not depend on how the knots, the weights, and the control points are precisely extended. We note that the assumptions w γ , 1 - p γ = w γ , N γ - p γ and C γ , 1 - p γ = C γ , N γ - p γ together with (2.2) below ensure that γ ⁢ ( a ) = γ ⁢ ( b - ) in the case of closed Γ = ∂ ⁡ Ω .

2.3 Rational Splines on Γ

Let p ∈ N 0 be an arbitrary but fixed polynomial degree. Let K ^ 0 = ( t 0 , i ) i = 1 N 0 be a sequence of initial knots with

and multiplicity # 0 ⁢ t 0 , i ≤ p + 1 for i ∈ { 1 , … , N 0 - p } . Let

denote the corresponding set of nodes with x ^ 0 , j - 1 < x ^ 0 , j for j ∈ { 2 , … , n 0 } . In order to apply standard quadrature rules, we assume that N ^ γ ⊆ N ^ 0 . Moreover, let W 0 = ( w 0 , i ) i = 1 - p N 0 - p be given positive weights. We extend the knots and weights as in Section 2.2 and define the weight function

where B ^ 0 , k , p are the standard B-splines defined in Section 2.1. We define W ^ 0 : [ a , b ] → R by continuously extending this function at 𝑏.

Now, let K ^ ⋆ = ( t 0 , i ) i = 1 N ⋆ be a finer knot vector, i.e., K ^ 0 is a subsequence of K ^ ⋆ which satisfies the same properties as K ^ 0 . Outside the interval ( a , b ] , we extend K ^ ⋆ as K ^ 0 . Let

denote the corresponding set of nodes on [ a , b ] with x ^ ⋆ , j - 1 < x ^ ⋆ , j for j ∈ { 1 , … , n ⋆ } , and let

denote the corresponding nodes on Γ. Note that N ^ ⋆ does not contain the node x ^ ⋆ , 0 = a , which is natural for the case Γ = ∂ ⁡ Ω , since then γ ⁢ ( x ^ ⋆ , 0 ) = γ ⁢ ( x ^ ⋆ , n ⋆ ) . We define the induced mesh on [ a , b ] ,

where we set x ^ ⋆ , 0 := a and the induced mesh on Γ, Q ⋆ := { γ ⁢ ( Q ^ ) : Q ^ ∈ Q ^ ⋆ } . Recall that the non-vanishing B-splines on [ a , b ) form a basis; see Section 2.1. This proves the existence and uniqueness of weights W ⋆ = ( w ⋆ , i ) i = 1 - p N ⋆ - p such that

Choosing these weights, we ensure that the denominator of the rational splines does not change. These weights are just convex combinations of the initial weights W 0 and can be computed via the knot insertion procedure (2.4). Finally, we extend W ⋆ arbitrarily to ( w ⋆ , i ) i ∈ Z with w ⋆ , i > 0 and define the space of (transformed) rational splines on Γ as

where S ^ p ⁢ ( K ^ ⋆ ) denotes the space of all right-continuous piecewise polynomials of degree lower than or equal to 𝑝 with breakpoints N ^ ⋆ on [ a , b ) which are, at each knot t ⋆ , i , p - # ⋆ ⁢ t ⋆ , i times continuously differentiable if p - # ⋆ ⁢ t ⋆ , i ≥ 0 ; see also Section 2.1. Here, we extend each quotient S ^ ⋆ / W ^ ⋆ left-continuously at 𝑏. The locality (2.1) of B-splines shows that this definition does not depend on how the knots and the weights are precisely extended. With the standard NURBS R ^ ⋆ , i , p from Section 2.1, a basis of S p ⁢ ( K ⋆ , W ⋆ ) is given by

where the functions R ^ ⋆ , i , p are again left-continuously extended at 𝑏.

3.1 Sobolev Spaces

The usual Lebesgue and Sobolev spaces on Γ are denoted by L 2 ⁢ ( Γ ) = H 0 ⁢ ( Γ ) and H 1 ⁢ ( Γ ) . With the weak arc-length derivative ∂ Γ , the H 1 -norm reads

We stress that H 1 ⁢ ( Γ ) is continuously embedded in C 0 ⁢ ( Γ ) . Moreover, we equip the space

with the same norm. We define the Sobolev space H 1 / 2 ⁢ ( Γ ) as the space of all functions u ∈ L 2 ⁢ ( Γ ) with finite Sobolev–Slobodeckij norm

We will also need the seminorm | u | H 1 / 2 ⁢ ( ω ) for subsets ω ⊆ Γ , which is defined analogously. Moreover, H ~ 1 / 2 ⁢ ( Γ ) := { v ∈ H 1 / 2 ⁢ ( ∂ ⁡ Ω ) : supp ⁢ ( v ) ⊆ Γ } is equipped with the same norm.

Sobolev spaces of negative order are defined by duality H - 1 / 2 ⁢ ( Γ ) := H ~ 1 / 2 ⁢ ( Γ ) * and H ~ - 1 / 2 ⁢ ( Γ ) := H 1 / 2 ⁢ ( Γ ) * , where duality is understood with respect to the extended L 2 ⁢ ( Γ ) -scalar product ⟨ ⋅ , ⋅ ⟩ Γ . Note that we have H ~ ± σ ⁢ ( ∂ ⁡ Ω ) = H ± σ ⁢ ( ∂ ⁡ Ω ) in case of Γ = ∂ ⁡ Ω even with equal norms for σ ∈ { 0 , 1 2 , 1 } .

All details and equivalent definitions of the Sobolev spaces are, for instance, found in the monographs [37, 32, 46, 42].

3.2 Weakly Singular Integral Equation

For σ ∈ { 0 , 1 2 , 1 } , the single- and double-layer operators V : H ~ σ - 1 ⁢ ( Γ ) → H σ ⁢ ( Γ ) and K : H ~ σ ⁢ ( Γ ) → H σ ⁢ ( Γ ) , respectively, are well-defined, linear, and continuous. As H 1 ⁢ ( Γ ) ⊂ C 0 ⁢ ( Γ ) , the case σ = 1 yields continuous functions, which is essential for the implementation described in Section 5.

For σ = 1 2 , V : H ~ - 1 / 2 ⁢ ( Γ ) → H 1 / 2 ⁢ ( Γ ) is symmetric and elliptic under the assumption that diam ⁢ ( Ω ) < 1 . In particular, ⟨ ϕ , ψ ⟩ V := ⟨ V ⁢ ϕ , ψ ⟩ Γ defines an equivalent scalar product on H ~ - 1 / 2 ⁢ ( Γ ) with corresponding norm ∥ ⋅ ∥ V . With this notation, the strong form (1.3) with data f = u for some u ∈ H 1 / 2 ⁢ ( Γ ) or f = ( 1 2 + K ) ⁢ u for some u ∈ H ~ 1 / 2 ⁢ ( Γ ) is equivalently stated as

The Lax–Milgram lemma applies, and hence (1.3) admits a unique solution ϕ ∈ H ~ - 1 / 2 ⁢ ( Γ ) . If f = u , the approach is called indirect; otherwise, if f = ( 1 2 + K ) ⁢ u , it is called direct. More details and proofs are found, e.g., in [37, 32, 46, 42].

3.3 Galerkin IGABEM

Let p ∈ N 0 be a fixed polynomial degree. Moreover, let K ^ ⋆ and W ⋆ be knots and weights as in Section 2.3. We introduce the ansatz space

Recall that

where the set even forms a basis. We define the Galerkin approximation Φ ⋆ ∈ X ⋆ of 𝜙 by

Note that (3.2) is equivalent to solving the finite-dimensional linear system

where Φ ⋆ = ∑ i ′ = 1 - p N ⋆ - p c ⋆ , i ′ ⁢ R ⋆ , i ′ , p . The computation of the matrix and the right-hand side vector in (3.3) is realized in VMatrix.c and RHSVectorWeak.c and can be called via buildVMatrix and buildRHSVectorWeak in Matlab; see Appendix A. For details on the implementation, we refer to Sections 5.1 and 5.2.

3.4 A Posteriori Error Estimation

Let K ^ ⋆ be a knot vector as in Section 3.3 with corresponding ansatz space X ⋆ . We define the mesh-size function h ⋆ ∈ L ∞ ⁢ ( Γ ) by h ⋆ | Q := | Q | for all Q ∈ Q ⋆ . Moreover, we abbreviate the node patch

We consider the following three different node-based estimators: the ( h - h 2 ) -estimator

where K ^ + are the uniformly refined knots obtained from K ^ ⋆ by adding the midpoint of each element Q ^ ∈ Q ^ ⋆ with multiplicity one (see Algorithm 3.1); the Faermann estimator

and the weighted-residual estimator

which requires the additional regularity f ∈ H 1 ⁢ ( Γ ) to ensure that f - V ⁢ Φ ⋆ ∈ H 1 ⁢ ( Γ ) and hence (3.6) is well-defined. If f = ( K + 1 2 ) ⁢ u stems from a Dirichlet problem as in Section 1.4, Section 3.2 shows that this is satisfied for u ∈ H ~ 1 ⁢ ( Γ ) .

The computation of the estimators is implemented in HHEstWeak.c, FaerEstWeak.c, and ResEstWeak.c and can be called in Matlab via buildHHEstWeak, buildFaerEstWeak, and buildResEstWeak; see Appendix A. The residual estimators require the evaluation of the boundary integral operators 𝔙 and 𝔎 applied to some function, which is implemented in ResidualWeak.c. The evaluations can be used in Matlab via evalV and evalRHSWeak; see Appendix A. The implementation of the estimators and the evaluations are discussed in Section 5.

3.5 Adaptive IGABEM Algorithm

To refine a given ansatz space X ⋆ , we compute one of the corresponding error estimators η ⋆ and determine a set of marked nodes with large error indicator η ⋆ ⁢ ( x ) . By default, we then apply the following refinement algorithm, which uses ℎ-refinement controlling the maximal mesh ratio in the parameter domain

as well as multiplicity increase.

The marked elements M ℓ ′ and the nodes whose multiplicity should be increased are determined in the function markNodesElements.m. The multiplicity increase and the computation of the new weights are realized in increaseMult.m. An optimal 1D bisection algorithm as in (iv) is discussed and analyzed in [4]. Together with the computation of the new weights, it is realized in refineBoundaryMesh.m. We stress that we have also implemented two further relevant strategies that rely on ℎ-refinement only: replace (iii) by adding all elements Q ∈ Q ⋆ containing one of the other nodes in M ⋆ to M ⋆ ′ ; or replace (iii) by adding all elements Q ∈ Q ⋆ containing one of the other nodes in M ⋆ to M ⋆ ′ and insert the midpoints in (iv) with multiplicity p + 1 . The first strategy leads to refined splines of full regularity, whereas the second one leads to lowest regularity.

We fix the considered error estimator η ⋆ ∈ { η V , hh2 , ⋆ , η V , fae , ⋆ , η V , res , ⋆ } . The corresponding error indicators η ⋆ ⁢ ( x ) are defined accordingly.

The adaptive algorithm is realized in the main function IGABEMWeak.m. The considered problem as well as the used parameters can be changed there by the user. In Section 5, we discuss how the arising (singular) boundary integrals are computed. We also mention that Algorithm 3.2 (iii) is realized in markNodesElements.m, which also determines the corresponding set of marked elements M ℓ ′ and the nodes whose multiplicity should be increased from Algorithm 3.1 (ii)–(iii).

3.6 Numerical Experiments

In this section, we empirically investigate the performance of Algorithm 3.2 for different geometries, ansatz spaces, and error estimators. Figure 1 shows the different geometries, namely a pacman geometry and a slit. The boundary of the pacman geometry can be parametrized via rational splines of degree 2, while the slit can be parametrized by splines of degree 1. As initial ansatz space, we either consider the same (rational) splines, i.e., K ^ 0 = K ^ γ and W 0 = W γ , splines of degree 𝑝 and smoothness C p - 1 , i.e.,

and W 0 = ( 1 , … , 1 ) , or piecewise polynomials of degree 𝑝, i.e.,

and W 0 = ( 1 , … , 1 ) . In the latter case, we always consider ℎ-refinement with new knots having multiplicity p + 1 as explained after Algorithm 3.1.

Figure 1

Geometries and initial nodes for examples from Section 3.6.