Multi-level Bézier extraction of truncated hierarchical B-splines for isogeometric analysis

Multivariate basis splines (B-splines) and Non-uniform rational B-splines (NURBS) lack adaptivity due to their tensor product structure. Truncated hierarchical B-splines (THB-splines) provide a solution for this. THB-splines organize the parameter space into a hierarchical structure, which enables efficient approximation and representation of functions with different levels of detail. The truncation mechanism ensures the partition of unity property of B-splines and defines a more scattered set of basis functions without overlapping on the multi-level spline space. Transferring these multi-level splines into Bézier elements representation facilitates straightforward incorporation into existing finite element (FE) codes. By separating the multi-level extraction of the THB-splines from the standard Bézier extraction, a more general independent framework applicable to any sequence of nested spaces is created. The operators for the multi-level structure of THB-splines and the operators of Bézier extraction are constructed in a local approach. Adjusting the operators for the multi-level structure from an element point of view and multiplying with the Bézier extraction operators of those elements, a direct map between Bézier elements and a hierarchical structure is obtained. The presented implementation involves the use of an open-source Octave/MATLAB isogeometric analysis (IGA) code called GeoPDEs. A basic Poisson problem is presented to investigate the performance of multi-level Bézier extraction compared to a standard THB-spline approach.


Introduction
Isogeometric analysis (IGA) [2,3] is a numerical method for solving partial differential equations (PDEs) by using a unified representation for both the analysis and design of structures.IGA is promising to close the gap between analysis and CAD geometries of the structures.Traditionally, engineers have used separate methods for modeling the geometry of a structure and analyzing its behavior.However, IGA breaks through this limitation by employing non-uniform rational Bsplines (NURBS) functions, a standard technology employed in CAD systems.The isoparametric philosophy is invoked, and the solution space for dependent variables is represented in terms of the same functions that represent the geometry.
Refinement techniques for NURBS functions are inherently simple and do not require further communication with the CAD system.Order elevation and knot insertion are refinement techniques analogous to standard FE methods like h-and p-, and a new, higher-order methodology emerges, k-refinement, which combines order elevation followed by knot insertion.All subsequent meshes retain the exact geometry of the initial CAD structure.
Due to the tensor product structure of the B-spline and NURBS functions, local refinement can not be achieved.Various approaches have been developed to overcome this limitation of tensor product structures, such as T-splines [4], LR-splines [5], and HB-splines [6].HB-splines are based on a multi-level structure with different levels of detail and appear promising to implement local refinement.Those splines can also be applied to T-splines indicating hierarchical T-splines [7].The hierarchical B-splines recently extended with the truncation mechanism, namely, THBsplines [8].This work presents an adaptive isogeometric analysis framework based on THBsplines.This type of function preserves all the essential mathematical properties of standard B-splines and retains the partition of unity property in a hierarchical structure.Despite the good properties of THB-splines and their simplicity, implementing the hierarchical definition of shape functions into existing FE solvers can be rather challenging.The proposed approach for this implementation is based on a generalization of the classical Bézier extraction [9], extended into a hierarchical concept [10,11].In fact, the multi-level Bézier extraction maps the multilevel hierarchical functions, and the corresponding knot spans, into a sequence of elements being equipped with a standard single-level basis.

Preliminaries
In this section, some preliminary concepts need to be discussed in order to give a summarized overview of the work steps.In particular, B-splines and NURBS, the idea of Bézier extraction, and hierarchical splines, together with the truncation mechanism, are presented.

B-splines and NURBS
We start by introducing the concept of B-splines and NURBS functions, but further details can be found in [2], [3] and [12].B-splines basis functions are defined by degree p, and a knot vector: where n is the number of basis functions.The Cox-de Boor recursion formula defines the n functions, by starting from degree p = 0, up to the prescribed degree, The multiplicity k of a knot value can be found from the knot vector and defines the smoothness of the B-spline basis at this location, e.g., C p−k .In the case of open knot vectors, the first and the last knot values have multiplicity k = p + 1, and the basis is interpolatory at the ends of its definition domain.A B-spline curve can be constructed as the linear combination of the basis functions.

Standard Bézier extraction
This section provides a short introduction to Bézier extraction [9].Unlike the standard finite element approaches, B-spline basis functions are defined globally (see Fig. 1).Bézier extraction provides a mapping between the smooth spline space and a C 0 representation of Bézier elements.Those elements are defined by Bernstein polynomials, which are local, i.e., their support is nonzero only in a single element.Thus, conventional element-based FE routines can be enriched by splines in a straightforward manner.Bézier extraction is generally based on knot insertion algorithms.The Bézier decomposition of the B-spline curve into piece-wise Bézier curves is obtained by increasing the multiplicity of each internal knot in Ξ, equals to the degree p. Starting by inserting a single knot ξ into Ξ, so Ξ = Ξ ∪ ξ, where ξ ∈ [ξ s , ξ s+1 ].This procedure continues until the desired multiplicity of the internal knots is achieved, resulting a new set of basis functions and a new set of control points that retain the geometry of the curve.Further details for knot insertion techniques can be found in [2].The B-spline curve C is described by a set of Bernstein polynomials B i,p and the Bernstein control points Pi The control points Pi can be computed as a linear combination of the original control points P i , which are represented by a matrix operator E, the Bézier extraction operator.The same operator is used for the translation of the basis functions N i,p (ξ) to the Bernstein polynomials B i,p (ξ).The basis functions can be written in a vector format with N and B, respectively.So, globally defined B-splines transformed into a basis of Bézier elements reads

Truncated hierarchical B-splines
In this section, the hierarchical B-splines are presented, which are used to achieve local refinement and overcome the limitation of adaptivity due to the tensor product structure of splines.Then, the extension with the truncation mechanism is added to retain the partition of unity property [8,10].Hierarchical B-splines [6,13] introduce a multi-level structure with different levels of refinement.Consider a sequence of nested spline spaces of univariate splines S (0) (Ω) ⊂ S (1) (Ω)... ⊂ S (ℓmax) (Ω) with the same polynomial degree, where ℓ = 0 . . .ℓ max are the levels of refinement.In this contribution, the different levels are constructed with the dyadic refinement technique, where each new level has 2d • n el elements (d is the dimension and n el number of elements in ℓ − 1).In Fig. 2, a univariate quadratic B-spline with two additional levels of refinement is depicted.The geometry and the parametric representation do not change under the refinement process, but the relation between the levels is happening by an operator called, subdivision matrix S. The parent functions of the level ℓ are a linear combination of fine-scale functions (children functions) of level ℓ + 1 Each nested knot vector Ξ ℓ defines a new set of B-spline basis functions N ℓ j,p (ξ).So, the multi-level spline space is defined as in [14] M S(Ω) := span{N Therefore, the refinement process starts by deactivating parent functions of the corresponding element from the coarse level ℓ and activating their children on level ℓ + 1, until the desired adaptivity is achieved, see Fig. 2. The selection of the coarse and fine elements results from error estimator algorithms or manual selection of the preferred elements.This information is stored in the subdivision matrices [14,15].The truncated basis was first time introduced in [8].The truncation mechanism preserves all the nice properties of hierarchical B-splines, such as linear independence, local support, and non-negativity.In addition, truncated hierarchical B-splines (THB-splines) have smaller support and retain the partition of unity property of traditional B-splines.The truncated basis functions are defined as In general, it is straightforward to identify the truncated B-spline functions.Let's assume an evaluation of an element on the finest level ℓ + 1.The basis functions that describe the element are N ℓ+1 i (ξ), with i = s ℓ+1 − p . . .s ℓ+1 , where the s ℓ+1 is the knot span of the element in that level.This element holds truncated B-splines on the previous level, ℓ, only if at least one of its N ℓ+1 i (ξ) is inactive (see Fig. 2).

Multi-level Bézier extraction
Having a sequence of nested spline spaces into a hierarchical structure, every basis function of level ℓ can be expressed as a linear combination of basis functions of level ℓ + 1, with the subdivision matrix (7).The global multi-level extraction operator, M glob ℓ , can be defined by applying a sequence of nested spline spaces, where the goal is to extract the hierarchical functions supporting an element at a specific level.For levels ℓ = 0, .., ℓ max the M glob ℓmax can be obtained by joining the rows, corresponding to the active basis, of the refinement operators S (ℓ,ℓmax) .More details for constructing the global multi-level extraction operator can be found in [11].Next, localization of the global multi-level extraction to each element of the domain leads to a local element-by-element approach.By selecting the appropriate columns and rows associated with the basis functions of the element, the extraction operator effectively captures the hierarchical functions.The multi-level extraction operator is always of the same level as the element.All the different matrices, M loc e , correspond to each element of the hierarchical spline space, e = 0, ..., n el .They can be combined directly with the Bézier extraction operator, resulting the multi-level Bézier extraction technique (see Fig. 3).This creates a direct map C e from a standard set of reference basis functions equal for each element (Bernstein basis, B) to the multi-level local basis (hierarchical basis).As the hierarchical basis functions are mapped into a local approach with elements with a standard set of basis functions, the calculations into the Bézier elements are performed.This can happen using the existing FE algorithms, and there is no need for separate algorithms to handle the globally defined spline basis functions.

Numerical example
This contribution uses GeoPDEs [1,16], an open-source FE solver for IGA, based on Matlab/Octave.Let us first highlight an advantage of Bézier extraction over the direct use of B-splines basis functions.The problem refers to the evaluation of points at the boundaries of the elements.Due to the Cox-de Boor formula (2), boundary points are generally associated with the right element of the domain, which can lead to misinterpretations.Consider the numerical integration of an element using a Newton-Cotes quadrature rule.These rules contain the boundary points, and based on (2), the contribution of the last point would be associated with the next element instead of the current one.This can be easily overcome with the use of Bézier extraction.An example to illustrate the problem is in Fig. 4, where the Newton-Cotes quadrature rule is applied for the integration of the elements in a univariate problem, instead of the Gaussian quadrature method.Bézier extraction should be applied to evaluate accurately.In this case, improving the results is not the goal of Newton-Cotes, but to verify the problem of spline-based algorithms using Cox-de Boor formula.

Unit square boundary value problem
The numerical example for this work is a second-order elliptic boundary value problem in a unit square domain.The functions f , g and h are chosen from the exact solution (see Fig. 5) of the problem u ex = exp (−C • ((x − 0.5) 2 + (y − 0.5) 2 )), where C is constant.
For the calculation, two types of THB-splines are assumed, one quadratic and one cubic.In both cases, global and local refinement techniques are compared.For the local refinement, the refined mesh is shown in Fig. 6.Now, as for the multi-level Bézier extraction, the goal is to approach the error of the local refinement that GeoPDEs is resulting in.Calculating the multi-level Bézier extraction operator C e , as described in Section 3, and then calculating the stiffness matrix of a single Bézier element, the stiffness matrix of an element in the hierarchical structure can be derived.Gathering all the element-wise stiffness matrices, the stiffness matrix for the whole multi-level spline space is obtained.Fig. 7 shows the convergence of the multi-level Bézier extraction and the local refinement of GeoPDEs.When comparing the values of the graphs in both cases, the results of multi-level Bézier extraction fit the IGA-based algorithm.This verifies that by just evaluating a single Bézier element, the stiffness matrix can be directly mapped into a hierarchical structure and evaluate multi-level spline spaces.

Conclusion
This study demonstrates that Bézier extraction can be extended into a multi-level structure of truncated hierarchical B-splines (THB-splines), resulting the multi-level Bézier extraction.The calculation over Bézier elements is sufficient to extract results for hierarchical spline spaces which perform local refinement.The only calculation that needs to be derived is the multi-level Bézier extraction operators for each element of the domain and the standard Bézier extraction operator.Then those operators are applied to Bézier elements.The implementation of the isogeometric analysis concept into existing finite element solvers can be employed.Additionally, this general approach is not restricted to THB-splines but can be compatible with any nested space.

Figure 1 :
Figure 1: Basis functions of a univariate quadratic B-spline and the corresponding curve.C 0 continuity at the knot = 0.75.

Figure 2 :
Figure 2: Basis functions for the computational domain.The gray functions are the deactivated ones and the colored ones are the active ones.The N (0) 3 and N (1) 7 are truncated functions.

Figure 3 :
Figure 3: Translate from hierarchical spline space to local Bernstein polynomials.

Figure 4 :
Figure 4: Newton-Cotes integration points appear in the boundary of the elements and can not be evaluated correctly due to the Cox-de Boor formula.Bézier extraction, in the right figure, overcomes this limitation with a standard set of Bernstein polynomials.

Figure 5 :
Figure 5: Exact solution of the unit square boundary value problem.

Figure 6 :
Figure 6: Mesh through iterations of local refinement with THB-splines.

Figure 7 :
Figure 7: L2 norm error.Global refinement, local refinement, and multi-level Bézier extraction on the local refinement.On the left are quadratic splines, and on the right are cubic splines.The green line of multi-level Bézier extraction matches exactly the blue line of the local refinement from the patch-wise approach.