ADAPTIVE ISOGEOMETRIC METHODS WITH HIERARCHICAL SPLINES: AN OVERVIEW

. We consider an adaptive isogeometric method (AIGM) based on (truncated) hierarchical B-splines and present the study of its numerical prop- erties. By following [10, 12, 11], optimal convergence rates of the AIGM can be proved when suitable approximation classes are considered. This is in line with the theory of adaptive methods developed for ﬁnite elements, recently well reviewed in [43]. The important output of our analysis is the deﬁnition of classes of admissibility for meshes underlying hierarchical splines and the design of an optimal adaptive strategy based on these classes of meshes. The adaptiv- ity analysis is validated on a selection of numerical tests. We also compare the results obtained with suitably graded meshes related to diﬀerent classes of admissibility for 2D and 3D conﬁgurations.


1.
Introduction. The design of efficient numerical schemes for the approximation of partial differential equations (PDEs) naturally demands suitable adaptive techniques, that automatically enable the refinement in well localized regions of the computational domain. The use of adaptive schemes is even more important in presence of singularities, when standard methods that do not support local mesh refinements do not achieve optimal convergence rates. In the context of isogeometric analysis [27], the design and development of adaptive methods necessarily require suitable spline spaces with local refinement capabilities.
The hierarchical spline model is one of the most elegant solutions to easily define adaptive spline constructions, see e.g., [31,49], and it has been successfully applied in challenging engineering applications [34,46,36]. In particular, truncated hierarchical B-splines (THB-splines) [24] properly meet the requirements needed for developing not only effective geometric modeling tools, but also efficient computational schemes [30,23]. Indeed, THB-splines have been succesfully used in the design and analysis of adaptive techniques for elliptic problems in [10,12,11], on which this review is based upon, while analogous schemes for hierarchical splines have been studied in [20] or in [8] as well as in the very recent contributions [1,41]. Adaptive techniques for boundary integro-differential equations have been introduced in [18,19]. The main difference between the approach presented in this review and the ones presented in the papers cited above is that our refinement strategy is based on THB-splines, while other approaches are based on the standard hierarchical basis. For this reason, our adaptive strategy may produce less refined meshes for the same error thresholds [6].
Besides hierarchical splines, other types of splines have been successfully used to design adaptive isogeometric methods. For example, the first use of adaptivity in the context of isogeometric methods was based on T-splines [47,2,15], see also [17,44], while their analysis suitable [35] or dual compatible [3,4] formulations have been used to design refinement strategies with linear complexity, as in [40,38]. Another example is the one of LR-splines, first defined in [14] and used in an adaptive setting in e.g., [28,32,33]. For a numerical comparison of all these approaches, we refer to [29,26]. We remark that, unlike for hierarchical splines, the mathematical theory of adaptive methods based on T-splines or LR-splines has not been developed so far.
For our analysis, we consider the elliptic model problem: − div(A∇u) = f in Ω, where Ω ⊂ R d , d ≥ 1, is a bounded domain with Lipschitz boundary ∂Ω, f is any square integrable function and the tensor A satisfies ∀x ∈ Ω, ξ ∈ R d η 1 |ξ| 2 ≤ A(x)ξ · ξ and |A(x)ξ| ≤ η 2 |ξ| with 0 < η 1 ≤ η 2 . Analogously to adaptive finite elements, the solution of the model problem with an adaptive isogeometric method (AIGM) is obtained by an iterative procedure, where each iteration of the adaptivity loop consists of the following four key steps: Following the guidelines from the theory of adaptive finite elements, see e.g. [42,43], the theory of convergence for adaptive isogeometric methods, using THB-splines with residual error estimators, has been recently developed in [10,12,11]. One of the key ingredients of the theoretical analysis is the definition of (strictly) admissible classes of hierarchical meshes. This class governs the maximum number of levels of THB-splines that do not vanish on any mesh element, and in practice it requires the mesh to be sufficiently graded.
In the present paper we review the convergence results developed in [10,12,11], considering and analyzing each of the four steps of the adaptive loop. We also complement the previous works by studying the performance of the method in a selection of numerical tests.
The content of the next sections is as follows. Some preliminary aspects of hierarchical B-spline constructions are reviewed in Section 2 together with the definition of THB-splines and related properties, before introducing the notion of (strictly) admissible meshes. The module SOLVE and ESTIMATE of the adaptive isogeometric method are discussed in Sections 3 and 4, respectively. To define the module SOLVE and compute the numerical solution of problem (1), we consider the Galerkin method on hierarchical spline spaces in connection with admissible mesh configurations. An a posteriori error analysis in terms of both local upper and lower bound for the energy error is also presented. Section 5 deals with the MARK and REFINE modules. We first recall a well-known marking strategy, and then introduce a refinement strategy that preserves the class of admissibility during the adaptive loop. Section 6 is devoted to the result of optimal convergence, while Section 7 presents some numerical examples to analyze how the class of admissibility influences convergence. We end with some concluding remarks.

2.
Hierarchical splines and admissible meshes. We consider a nested se- Let B be the tensor-product B-spline basis of level , = 0, 1, . . . , N − 1, and degree p = (p 1 , . . . , p d ), defined on the grid G . In each coordinate direction, the knot vectors at level zero are assumed to be open so that the first and the last knots are repeated p i + 1 times, for i = 1, . . . , d. B-splines are the standard form for spline representation in computer aided (geometric) design, see e.g., [13,45]. They are locally linearly independent and non-negative. In addition they have local support and form a partition of unity. When nested spline spaces are considered, there exists a two-scale relation so that any function s ∈ V ⊂ V +1 can be expressed as in terms of the coefficients c +1 β .
Let Ω 0 ⊇ Ω 1 ⊇ . . . ⊇ Ω N −1 be a nested sequence of closed subsets of D, such that Ω = (3) be the collection of active (i.e., no more refined) elements at level and the hierarchical mesh, respectively. For any element Q of the hierarchical mesh Q, we assume that where h Q := | Q| 1/d . Note that the symbol is used here and it what follows to denote an inequality which does not depend on the number of hierarchical levels. A mesh Q * is a refinement of Q, usually indicated as Q * Q, if each element Q * ∈ Q * either also belongs to Q or is obtained by refinement of an element of Q. The hierarchical B-spline basis can be defined as follows, see also [31,49].
Definition 1. The hierarchical B-spline (HB-spline) basis H with respect to the mesh Q is defined as Hierarchical B-splines are non-negative, linearly independent, and define nested hierarchical spline spaces, see e.g., [49]. The study of the hierarchical spline space was presented in [22] and [37] for the bivariate case with single knots and the general multivariate setting with arbitrary knot multiplicities, respectively. By considering suitable mesh configurations, it was there proved that span H( Q) contains all piecewise polynomial functions defined over the hierarchical mesh. From a practical point of view, uniform knot vectors and a fixed degree are usually considered at all levels, however the hierarchical B-spline model can be also applied with non-uniform mesh configurations and increasing degrees, as long as the tensor-product spaces in the sequence remain nested.
The following definition introduces the truncated basis for hierarchical splines [24] and hinges on the notion of the so-called truncation, namely where c +1 β (s) is the coefficient of the function s ∈ V with respect to β introduced in equation (2). By starting from the two-scale relation (2), the truncation of the function s with respect to level + 1 defined by (5) considers only the B-splines of level + 1 which do not belong to H( Q). Figure 1 shows an example of truncation for a univariate quadratic B-spline. When iteratively applied to hierarchical basis functions, the truncation eliminates from coarser B-splines the contribution of B-splines introduced at subsequent refinement levels, according to the following definition. where Trunc +1 β := trunc N −1 (trunc N −2 (. . . (trunc +1 ( β)) . . . )), for any β ∈ B ∩ H( Q).
The level of a THB-spline τ = Trunc +1 β in T ( Q) is the level of its corresponding mother function β. We will denote H = H( Q), T = T ( Q) when there will be no ambiguity in the text. THB-splines span the same space of HB-splines, are nonnegative, linearly independent, and form a partition of unity [24]. In addition, they are a strongly stable basis, see [25] for the details. The two properties of nonnegativity and partition of unity imply the convex hull property, which facilitates adaptive modeling operations [30]. Note that the reduced support of THB-splines with respect to HB-splines decreases the overlapping of basis functions introduced at different levels of the spline hierarchy.
In order to develop an adaptivity theory that exploits the reduced support of THB-splines with respect to standard hierarchical B-splines, the notion of admissible meshes was introduced in [7].

Definition 3.
A mesh Q is admissible of class m if the truncated basis functions in T ( Q) which take non-zero values over any element Q ∈ Q belong to at most m successive levels.
When an admissible mesh Q of class m is considered, the number of THB-splines in T ( Q) which are non-zero on each mesh element is bounded. In particular, it is less than m Note that the hidden constants in the above inequalities depend on m but not on τ , Q or N .
To easily link the local support of THB-splines to a certain subset of admissible meshes, we extend the definition of the support extension of a mesh element to the hierarchical setting as follows.
Definition 4. The support extension S( Q, k) of an element Q ∈ G with respect to level k, with 0 ≤ k ≤ , is defined as To keep the notation as simple as possible, we will also denote by S( Q, k) the region occupied by the closure of elements in S( Q, k). By considering the auxiliary subdomains for = 0, . . . , N − 1, we arrive at the following proposition.
Proposition 5. Let Q be the mesh of active elements defined according to (3) with respect to the domain hierarchy If for = m, m + 1, . . . , N − 1, then the mesh Q is admissible of class m.
The proof of this result can be found in [10]. As we will see in Section 5, it is possible to define refinement algorithms such that the constructed meshes satisfy the assumptions in the proposition. This relevant set of admissible meshes, called strictly admissible meshes, is defined as follows. Remark 7. In [20] admissible meshes for hierarchical B-splines, instead of THBsplines, were considered by limiting the interaction between the basis functions acting on any mesh element to 2 consecutive values (analogous to consider m = 2 with respect to HB-splines in Definition 3). This kind of admissible refinements for m ≥ 2 were also considered in [39]. A general framework for the design and implementation of (strictly) admissible refinements for HB-and THB-splines was recently presented in [6]. The properties of the hierarchical mesh configurations obtained with these algorithms were there analyzed and compared.
3. The module SOLVE. Let T 0 be the truncated basis defined on an initial strictly admissible mesh Q 0 . We assume that the computational domain is defined as x We also assume that the mapping F : Ω 0 → Ω is a bi-Lipschitz homeomorphism: where c F and and C F are independent constants bounded away from infinity. To define the variational formulation of problem (1), we consider the space of functions in H 1 (Ω) with vanishing trace on ∂Ω where a : and f ∈ L 2 (Ω). Due to the assumptions on A, the bilinear form a(u, v) is coercive and continuous: with constants α 1 and α 2 , respectively. In addition, it induces the energy norm |||v||| Ω := a(v, v) 1/2 , ∀v ∈ V. The coercivity and continuity properties of a(u, v) imply the equivalence between the energy and the H 1 (Ω) norms on V. The Lax-Milgram theorem ensures the existence and uniqueness of the weak solution (7). All the notation introduced in the previous section for the parametric domain is now used in the computational domain by simply removing the ·. In particular, for any admissible mesh Q Q 0 and truncated basis T ( Q), the corresponding mesh and functions of the physical domain are defined as The basis T (Q) collects the mapped THB-splines with respect to the hierarchical mesh Q on the domain Ω, and S(Q) := span T (Q). We denote by Q the mesh element Q = F( Q), and by h Q = |Q| 1/d , where |Q| represents the volume of Q, its size. In view of (4) and (6), we have: h Q diam(Q) h Q . We also set: Analogously, the support extension of Q ∈ G with respect to level k is given by Finally, when Q is a refinement of Q, we will write Q Q, when their preimages Q and Q satisfy Q Q. We define the module SOLVE of the AIGM as the Galerkin discretization of (7) in terms of hierarchical splines on Ω, that corresponds to where S D (Q) = S(Q) ∩ H 1 0 (Ω). For simplicity, even if not strictly needed, we assume S D (Q) ⊂ C 1 (Ω). The general case could be treated on the line of the classical theory of adaptive finite element methods. 4. The module ESTIMATE. Let u be the exact weak solution of the model problem (7). The residual r, v := f, v − a(U, v) associated to U ∈ S D is the functional in the dual space to V that satisfies The module ESTIMATE of the AIGM computes the error indicator which is defined in terms of the element residual r = f + div(A∇U ). Note that the residual does not contain any edge contribution as in a typical finite element setting, due to the assumption S D (Ω) ⊂ C 1 (Ω). The a posteriori error analysis of the adaptive isogeometric methods was presented in [10] leading to the upper and lower bound for the Galerkin error: where the oscillations are defined as and Π n : L 2 (Q) → Q n , n = (n 1 , . . . , n d ), denotes the L 2 projector onto the space of polynomials of degree n j in the space direction j. Note that the hidden constants in (11) do not depend on the mesh size and the hierarchical level. While (11) directly leads to a local version of the lower bound, namely a separate study is necessary to derive a local upper bound. Let Q and Q * be two admissible meshes so that Q * Q. The corresponding Galerkin solutions U ∈ S D (Q) and U * ∈ S D (Q * ) of problem (9) satisfy where R := R Q→Q * is the refined set of elements, namely the elements of Q that do not belong to Q * , and Ω R := Q : Q ∈ R .
Remark 8. The local version of the upper bound in (12) was derived in [11] by exploiting a well selected quasi-interpolation operator I Q on hierarchical spline spaces that satisfies for w ∈ S D (Q * ). In order to do this, a class of L 2 -stable quasi-interpolation operators onto the space of splines on tensor-product meshes [9] was suitably combined with hierarchical quasi-interpolants expressed in terms of the truncated basis [48].

5.
The modules MARK and REFINE. The module MARK of the AIGM selects and marks a set of elements M ⊂ Q according to Dörfler's marking [16]. A set of marked elements M such that is suitably selected by considering a fixed parameter θ ∈ (0, 1]. Even if for the convergence of method, the cardinality of this set is not crucial, it plays a key role for deriving the optimality of the method in Section 6. In particular, its minimal cardinality is required in Lemma 14 below where a suitable bound of the number of elements marked at the iterative step k is provided. The REFINE module is based on a refinement strategy which exploits the properties of THB-splines for obtaining strictly admissible meshes between two successive steps of the adaptive loop. Note that the strict version of admissibility is considered in order to keep the refinement routine as simple as possible while still exploiting the effect of the truncation on the support of THB-splines. By recalling the definition of the support extension S(Q, k) of an element Q with respect to level k introduced in (8), the refinement should be recursively propagated to a certain neighborhood of every marked element. The automatic refinement of the AIGM is based on the definition of the RE-FINE and REFINE RECURSIVE modules presented in Figure 2, that were first introduced in [10]. The fundamental properties of these algorithms, summarized in Lemma 10 and Proposition 11 below, were analyzed and proved in [10].
A last remark concerning the hierarchical refinement should be devoted to the overlay mesh Q * := Q 1 ⊗Q 2 of two meshes Q 1 , Q 2 obtained as the coarsest common refinement of Q 1 and Q 2 . When strictly admissible meshes are considered, the overlay is still strictly admissible [12] and its cardinality satisfies [5,40] where Q 0 is the initial mesh configuration. Even if both the complexity estimate (13) and the overlay inequality (14) were obtained on the parametric domain Ω in [12], they also hold on physical meshes defined as images of parametric meshes. characterizes the quality of the best approximation in Q m M with respect to the best total error defined as A complete analysis of these approximation classes is still missing and goes beyond the scope of this paper. In particular, the connection and possible dependence of the class of approximation on the value of m have not been identified.
In the paper [11], the following two additional results associated to the considered AIGM were recently given.
Lemma 13. (Quasi-optimality of total error) Let Q ∈ Q m be a strictly admissible mesh. The total error associated to the Galerkin solution U ∈ S D (Q) of problem for any k ≥ 0.
The quasi-optimality result, summarized in the next theorem, was also proved in [11] by exploiting Lemmas 13 and 14, together with the complexity estimate (13), the (global) lower bound in (11), and the contraction property [10].
Theorem 15. Let the marking parameter θ satisfy θ ∈ (0, θ * ) with θ * small enough, and assume that the module MARK selects a set M k of marked elements with minimal cardinality. Let u be the solution of the model problem (1). If (u, f, A) ∈ A s , the AIGM generates a sequence {Q k , S D (Q k ), U k } k≥0 of strictly admissible meshes, hierarchical spline spaces, and discrete solutions so that for any k ≥ 1.

Numerical examples.
In this section we show the performance of the adaptive isogeometric method by applying it to Poisson problem where Γ D ∪Γ N = ∂Ω, Γ D ∩Γ N = ∅ and Γ D = ∅. We note that, to deal with Neumann boundary conditions, the error indicator must be modified in the following way where r is defined as in (10), and r N = ∂U/∂n − g N . Although Neumann and non-homogeneous Dirichlet conditions are not covered by our theory, we have included such conditions in some of the numerical tests, since their main purpose is to analyze the behavior of the adaptive method when considering different degrees and different values of the admissibility class m. All the results are obtained by using the isogeometric Matlab/Octave GeoPDEs library, see [21]. We consider three different numerical tests with different kind of solutions: the first one with smooth solution, and the other two with singular solution, for which local refinement is required to obtain optimal convergence rates. For all the numerical tests, we discretize with hierarchical B-splines of degrees p = (p, p) = (2, 2), (3,3), (4,4) and we consider the admissibility classes m = 2, 3, 4, ∞, where m = ∞ corresponds to strictly consider the set of marked elements without refining any elements in the neighborhood in the REFINE RECURSIVE module. The adaptivity algorithm is stopped when the hierarchical mesh reaches a certain number of levels (seven in the first example, nine in the second and eight in the third). At each step k ≥ 0 of the adaptive procedure, for the approximate solution U k in the hierarchical space S(Q k ), we compute the error estimator ε Q k (U k , Q k ) and, when an exact solution is available, the H 1 -seminorm of the error and the effectivity index I k e , defined as which is shown on the left of Figure 3. The solution is a smooth function with a peak in the center of the domain. For the adaptive algorithm the first iteration is computed in a Cartesian mesh of 4 × 4 elements, and Dörfler's marking strategy is performed setting the parameter θ = 0.25. On the left of Figure 4 we show the behavior of the error and of the estimator during adaptive refinement. With any of the considered values of degree p and admissibility class m we obtain the optimal convergence rate O(N −p/2 d ), where N d = #T (Q) is the number of degrees of freedom. Notice that since the solution is smooth, optimal convergence rate is also obtained with uniform refinement, but the adaptive method always requires less degrees of freedom.
On the right of Figure 4 we plot the effectivity index divided by 10. The effectivity index remains bounded from above and from below for all the values of the admissibility class, which indicates the good behavior of the estimator in all cases.
Although the convergence rates are optimal for all cases, independently of the admissibility class, there is a relative impact of the value of m in the number of refined elements, and consequently in the number of degrees of freedom. As expected, when the value of m is increased the refinement is more localized in the center of the domain, see Figure 5. This loss of locality for low values of m is also observed in the error plot of Figure 4, as the curve for m equal 2 remains slightly above those for m equal to 3 and 4. We also report, in Table 1, the total number of elements and degrees of freedom at the last iteration.  Table 1. Number of elements (#Q) and degrees of freedom (#T (Q)) for Example 1.
Example 2: Singular function. For the second numerical test, the domain is again the unit square Ω = (0, 1) 2 , and we impose homogeneous Dirichlet boundary conditions on Γ D = ∂Ω, and a source function f such that the exact solution is given by , which is shown on the right of Figure 3. The adaptive algorithm is run with the same parameters of the previous example. In this case the solution belongs to H s (Ω), for some 2 < s < 3, with two different singularities on the edges x = 0 and y = 0. Hence, uniform refinement will not attain optimal convergence rates for high degree p, and local refinement near the two edges is necessary. This expected behavior is indeed observed in the convergence errors of Figure 6 (left). Moreover, the optimal convergence rate for adaptive refinement is obtained for all the values of the admissibility class m, with a slight improvement in terms of the number of degrees of freedom for higher m, as in the previous example.
On the right of Figure 6 we show the effectivity index divided by 10. As in the previous numerical test, the effectivity index remains bounded from above and from below for any degree and any value of m, which indicates that the estimator works properly both for smooth and singular solutions.  We report in Table 2 the number of elements and degrees of freedom obtained at the last iteration of the adaptivity procedure. As expected, much coarser meshes can be used for high degree. This can be also verified from the plots of the hierarchical meshes in Figure 7. Moreover, increasing the value of m also reduces the number of elements, since refinement is less spread by the REFINE procedure.   Figure 6. Numerical error and estimator (left) and effectivity index (right) for the example with singular solution. On the left, the error (solid lines) and the estimator (dashed lines) are plotted for different degrees. The red, green, blue and black lines (star, diamond, cross and circle markers, respectively) represent admissibility classes m = 2, 3, 4, ∞, respectively, while the magenta line (square marker) corresponds to uniform refinement. The same coloring and marking is used for the effectivity index (divided by 10) on the right figure.
In both examples any value of m provides the optimal convergence rate. In order to maximize the efficiency of the adaptive method, and to keep the number of degrees of freedom as small as possible, the numerical tests suggest that higher values of m should be used. We remark that in these examples optimal convergence is also obtained for m = ∞, but it is important to note that this is not guaranteed by the theory.
Example 3: 3D example. For the last example we consider a three dimensional domain with non-trivial geometry. The domain is obtained by linear interpolation of two surfaces, where the first one is a quarter of a ring with inner and outer radius equal to one and two, respectively, while the second one is the same surface rotated 90 degrees around the z axis, and then translated by the vector (0.5, 0, 1) (see Figure 8(a)). We consider a null source function, and impose a Neumann condition on the upper boundary with g N = 1, and homogeneous Dirichlet conditions elsewhere. This generates singularities on the edges of the upper boundary. In this case the exact solution is not known, but an approximation is shown in Figure 8(b). For the adaptive algorithm we start with a mesh formed by one single element, and refine using Dörfler's marking strategy with a parameter θ = 0.75. Since the exact solution is not known, we only present in Figure 9 the convergence of the error indicator. The adaptive method attains optimal convergence rates in terms of the degrees of freedom, which is now O(N −p/3 d ). Moreover, as in the previous examples better convergence is obtained with higher values of the admissibility class m, and this good behavior is even more evident in the current three-dimensional example than in the previous ones. Finally, we show in Figure 10 the restriction to the boundary of the final hierarchical mesh, for degree p = (4, 4) and for all the chosen admissibility classes. The adaptive method localizes the refinement near the edges where the singularity is more pronounced. As in the two-dimensional examples, a lower value of m reduces the locality of the refinement, which propagates to the interior.
8. Closure. This paper is largely inspired by [11,10,12] and provides a comprehensive analysis of adaptive isogeometric methods based on truncated hierarchical B-splines with a residual-based error estimator. As it is natural, at each adaptive step, in order to restore the good properties of the mesh, i.e. to maintain its class of admissibility, refinement is done "around" marked elements. Although the complexity estimates derived in [12] prove that this procedure still enjoys optimal complexity, in practice we add a non-negligible number of degrees of freedom, and the question whether this is really needed remains open. Our numerical results give some insight about this problem: the larger m is chosen the less the proliferation of degrees of freedom is obtained while preserving optimal convergence rates. Additional numerical investigations on the best values of m and on the impact of such a proliferation on the computational cost for a given tolerance are object of ongoing research.    (c) Estimator for p = (4, 4) Figure 9. Convergence of the estimator for the 3D example. The red, green, blue and black lines (star, diamond, cross and circle markers, respectively) represent admissibility classes m = 2, 3, 4, ∞, respectively, while the magenta line (square marker) corresponds to uniform refinement.