A simple strategy for defining polynomial spline spaces over hierarchical T-meses

A simple strategy for defining polynomial spline spaces over hierarchical T-meshes M. Brovkaa, J.I. L´opeza, J.M. Escobara, R. Montenegroa, J.M. Casc´onb aUniversity of Las Palmas de Gran Canaria, University Institute for Intelligent Systems and Numerical Applications in Engineering (SIANI), Spain. bUniversity of Salamanca, Department of Economics and Economic History, Faculty of Economics and Business, Spain. Abstract We present a new strategy for constructing spline spaces over hierarchical T-meshes with quadand octree subdivision scheme. The proposed technique includes some simple rules for inferring local knot vectors to define C2-continuous cubic tensor product spline blending functions. Our conjecture is that these rules allow to obtain, for a given T-mesh, a set of linearly independent spline functions with the property that spaces spanned by nested T-meshes are also nested, and therefore, the functions can reproduce cubic polynomials. In order to span spaces with these properties applying the proposed rules, the T-mesh should fulfill the only requirement of being a 0balanced mesh. The straightforward implementation of the proposed strategy can make it an attractive tool for its use in geometric design and isogeometric analysis. In this paper we give a detailed description of our technique and illustrate some examples of its application in isogeometric analysis performing adaptive refinement for 2D and 3D problems. Keywords: isogeometric analysis, multivariate splines, local refinement, T-mesh, nested spaces. 1. Introduction The main drawback of using B-splines and NURBS for geometric design is the impossibility to perform local refinement. T-splines were introduced by Email address: [email protected] (M. Brovka) Preprint submitted to Computer-Aided Design June 5, 2015 Sederberg et al. [1] as an alternative to NURBS. Based on the idea of admitting meshes with T-junctions and inferring local knot vectors by traversing T-mesh edges, T-splines have provided a promising tool for geometric modelling that allows to perform local refinement without introducing a large number of superfluous control points. Later, in [2], T-splines were incorporated to the framework of isogeometric analysis. Isogeometric analysis (IGA) was introduced in 2005 by Hughes et al. in [3, 4]. It has arisen as an attempt to unify the fields of CAD and classical finite element methods. The main idea of IGA consists in using for analysis the same functions that are used in CAD representation of the geometry. To use spline functions for numerical analysis and obtain a proper convergence behaviour, these functions must meet some requirements: linear independence, polynomial reproduction property, local supports and possibility to perform local adaptive refinement. This issue has been the object of numerous research works in recent years. Analysis-suitable T-splines, proposed by Scott et al. in [5], are a class of T-splines defined over T-meshes that should meet certain topological restrictions formulated in terms of T-junction extensions. Basis functions defined over an extended analysis-suitable T-mesh are linearly independent [6] and form a partition of unity. The refinement algorithm allows to accomplish highly localized refinements and constructs nested T-spline spaces, but it presents an elevated implementation complexity and, as far we know, the generalization of the strategy to 3D cases is still an open question. Another approach to the problem of local enrichment of the approximation space is the hierarchical refinement, originally introduced by Forsey and Bartels in [7] and later developed in [8]. Recently, hierarchical refinement technique in the context of isogeometric analysis was described in [9, 10, 11]. This approach is based on a simple and natural idea to construct multilevel spaces by replacing coarse level functions with finer basis functions. Starting from an initial uniform mesh, hierarchical refinement scheme leads to sequential construction of nested spline spaces with linearly independent basis functions. Simplicity of its implementation and straightforward generalization to 3D make it an attractive option for local refinement. However, a drawback of this strategy is the impossibility to define a spline space over a given arbitrary T-mesh as well as the presence of redundant basis functions and excessive support overlapping. An interesting theoretical approach to the latter problem was given in [12]. The truncation technique is applied to redefine the function supports and reduce their overlapping. 2 Other options for performing local refinement of spline spaces are C1continuous PHT-splines [13] or local refined splines (LR-splines) [14]. In the present paper we propose another possible alternative for the construction of spline functions that span spaces with nice properties. The technique we present here is designed for hierarchical T-meshes (multilevel meshes) with a quadand octree subdivision scheme. This type of meshes can be efficiently implemented with tree data structures [15] which are frequently used in engineering. Due to the elevated complexity of all current strategies, the main goal we pursue here is the simplicity and low computational cost of the implementation, both in 2D and 3D. For that, we have to assume a restriction on the T-mesh. Namely, the T-mesh should fulfill the requirement of being a 0-balanced mesh. A balanced mesh condition is usually imposed to have gradual transition from the coarse mesh to the finer zones. In addition, for our technique, this condition is an obligatory prerequisite that the T-mesh should fulfill. Assuming this reasonable restriction over the T-mesh, we can define easily cubic spline functions that span spaces with desirable properties: linear independence, C2-continuous, cubic polynomial reproduction property, nestedness of spanned spaces and a straightforward implementation. The key of the strategy lies in some simple rules used for inferring local knot vectors for each blending function. The paper is organized as follows. Some basic concepts about B-splines and T-meshes are given in Section 2. Section 3 includes the general scheme of our strategy and the description of its main stages. In Section 4 we explain in detail the key of our technique, that is, the rules used for inferring function supports in order to span spaces with desirable properties. In Section 5 the properties of the defined functions are given and the issue of support overlapping and sparsity of stiffness matrix is discussed. Computational examples of performing adaptive refinement for 2D and 3D Poisson problems are presented in Section 6. Conclusions are given in Section 7. 2. Basic concepts 2.1. B-spline basis functions We start with a brief summary of the main concepts about B-splines. A set of B-spline basis functions Bi,p (i= 1,2, . . . , n) of degree p, inferred from a non-decreasing sequence Ξ = {ξ1, ξ2, . . . , ξn+p+1}, called knot vector, is defined by the Cox-de Boor recursion formula 3 Bi,0(ξ) = {1 if ξi≤ξ < ξi+1, 0 othewise. Bi,p(ξ) = ξ−ξi ξi+p−ξi Bi,p−1(ξ) + ξi+p+1 −ξ ξi+p+1 −ξi+1 Bi+1,p−1(ξ). A knot vector Ξ is called open knot vector if the first and the last knots are repeated p+ 1 times. At each knot of multiplicity mthe basis functions are Cp−m. A B-spline curve is defined as a linear combination of B-spline basis functions S(ξ) = ∑ i∈I PiBi,p(ξ), where coefficients Pi∈Rsare called control points, typically s= 2 or 3. Multivariate B-splines are defined as a tensor product of univariate B-spline functions Bi,p(ξ) = d ∏ k=1 Bik,p(ξk), where ξ= (ξ1, . . . , ξd) and the multi-index i= (i1, . . . , id)∈I. The multi-index set is defined by I={1,2, . . . , n1} × · · · × {1,2, . . . , nd}. A B-spline surface (solid) is defined as a linear combination of bivariate (trivariate) B-spline functions S(ξ) = ∑ i∈I PiBi,p(ξ), where the control points Pi∈Rs(s= 2 or 3) form a control mesh. For more details about B-splines see [16]. 2.2. T-meshes and T-splines In order to overcome the drawback of tensor product structure, which does not allow to perform local refinement, it is necessary to admit T-junctions in the mesh. The concept T-junction is similar to hanging node in the classical finite element method. An axes-aligned grid that allows T-junctions is called T-mesh. As was mentioned in the previous section, there are several strategies to define tensor product spline functions over T-meshes and one of these 4 Figure 1: An example of T-mesh and inferring of local knot vectors for cubic T-spline functions by traversing T-mesh edges. strategies is T-splines. The underlying idea of T-splines consists in defining blending functions by means of a set of local knot vectors instead of a global knot vector, as in the case of B-splines or NURBS. A local knot vector for each bivariate function Bαis inferred by traversing the T-mesh edges in both parametric directions starting from a vertex vαof the mesh (the anchor), see Fig. 1. For a pair of local knot vectors Ξj α=(ξj 1, ξj 2, ξj 3, ξj 4, ξj 5), j = 1,2 the bicubic spline function Bαis defined as Bα(ξ1, ξ2) = B[Ξ1 α](ξ1)B[Ξ2 α](ξ2),where B[Ξj α](ξj) is an univariate B-spline corresponding to the knot vector Ξj α. In general, T-spline blending functions do not span a polynomial space. Some additional restrictions on the T-mesh configuration [5] should be satisfied in order to span a polynomial spline space. If these restrictions are not verified, the T-splines should be normalized in order to form a partition of unity. This leads to rational blending functions: Rα(ξ1, ξ2) = Bα(ξ1, ξ2) ∑β∈ABβ(ξ1, ξ2), where Ais the index set of the basis spanned by the T-mesh. These rational blending functions are capable of reproducing a constant function, but, in general, can not reproduce a polynomial of a higher order. A T-spline approximation is constructed as a linear combination of all blending functions: S(ξ1, ξ2) = ∑α∈APαRα(ξ1, ξ2). 5 (a) (b) (c) Figure 2: Motivation of the strategy. (a) Initial mesh and a basis function Nα. (b) Refined mesh, where the new T-spline space cannot reproduce the original function Nα. (c) We try to preserve the original basis function Nαin the new space if it cannot be recovered. 3. Strategy for the construction of polynomial spline spaces over hierarchical T-meshes. In this section we describe our strategy to define tensor product spline functions over hierarchical T-meshes. The strategy we propose has some similarity with T-splines inasmuch as we define the blending functions from local knot vectors that are inferred by traversing the T-mesh edges. Some additional rules and requirements are imposed for the local knot vectors in order to obtain spline spaces with nice properties. The strategy is mainly motivated by the idea of preserving an original basis function if it cannot be reproduced with the basis of the new space after a refinement. For example, in the initial mesh of the Fig. 2(a), we define a basis function Nαassociated to the vertex α. Then, we perform some cell refinements as shown in Fig. 2(b). It is easy to check that the new T-spline space obtained in Fig. 2(b) is not capable of reproducing the original basis function Nα. However, if we define the function ˆ Nαas shown in Fig. 2(c), keeping in this case the initial basis unaltered in the new refined space, we guarantee the nestedness of spline spaces. In our strategy, the process of spline space construction for a given Tmesh can be divided in the following three steps: 1. Mesh pretreatment (0-balancing) 2. Inferring local knot vectors 3. Modification of local knot vectors 6 (a) (b) Figure 3: Example of 0-balancing procedure. (a) A 0-unbalanced quadtree mesh. Cells to be refined are shaded in blue. (b) Resulting 0-balanced mesh. Next, we give a description of each step of the process. 3.1. Mesh pretreatment. 0-balanced quadtree and octree T-meshes. Due to its simplicity, quadtree and octree meshes are an attractive tool for performing adaptive refinement in IGA and geometric modelling. To guarantee a good quality of the approximation space constructed over a mesh, it is preferable to have a gradual transition from the coarse mesh to the finely refined zone. That is why it is common to work with balanced quadtree and octree meshes. The strategy we propose in this paper is designed exclusively for the 0-balanced T-meshes. A mesh with tree structure is said to be 0-balanced if for any k, no cell at level kshares a vertex (0-face) with a cell at level greater than k+ 1. In other words, a 0-balanced quadtree mesh implies that any cell has contact (through vertex, edge or face) only with cells that differ at most in one level of depth. An example of 0-balanced procedure quadtree is shown in Fig. 3(b). To obtain a 0-balanced quadtree, a standard balancing procedure is applied. Note that refinements performed during the 0-balancing procedure do not propagate, see [17]. It should be highlighted that 0-balancing the T-mesh is an essential prerequisite for the construction of spline spaces by means of our technique. In general, if the T-mesh is not 0-balanced, our rules for inferring local knot vectors do not lead to polynomial spaces. Also, it is important to emphasize that, for our 2D (3D) T-meshes, a subdivision of any cell is performed by subdividing the cell in 4 (8) equal subcells so that, all cells of the same level 7 (a) No T-junction is skipped. (b) Two T-junctions are skipped in each direction. (c) Two T-junctions are skipped and the knots are repeated when the boundary is reached. Figure 4: Inferring local knot vectors for a bivariate function by traversing the T-mesh edges. have the same size and the edge size of a k-level cell is twice larger than the edge size of a (k+ 1)-level cell. 3.2. Inferring local knot vectors Let us consider a T-mesh Tof the squared parametric domain Ω = [0,1]d, d= 2 or 3. We call regular node the node of the mesh that is not a Tjunction. We associate a blending function only to regular nodes of the mesh, as it is usual in classical finite element methods when working with hanging nodes. The skeleton of a d-dimensional mesh Tis the geometric set of points composed of the union of all (d−1)-faces of the mesh and it is denoted by skt(T). That is, for a 2D space, the mesh skeleton is the union of all the edges of the mesh, and the skeleton of a 3D mesh is the union of all its faces. To define our cubic tensor product spline blending functions over a given d-dimensional T-mesh, a local knot vector for dparametric directions should be assigned to each function Nα: Ξj α=(ξj 1, ξj 2, ξj 3, ξj 4, ξj 5),j= 1, ..., d. Similarly to [1], these knot vectors are inferred by traversing the T-mesh skeleton. For simplicity, let us describe this procedure for a two-dimensional T-mesh. Starting from the central knot (ξ1 3, ξ2 3), i.e., the anchor of the function, we walk across the T-mesh until intersecting perpendicularly a mesh edge. According to our strategy, we should skip over the T-junctions where the missing edge is perpendicular to the direction of our marching, see Fig. 8 4. When the boundary of the parametric domain is reached while walking across the mesh, we repeat knots creating an open knot vector structure along the boundary, see Fig. 4(c). Note that all interior knots have multiplicity 1. Thus, we obtain for mesh Ta set of blending functions {Nα}α∈AT, where ATis the index set. Figure 5 illustrates an example of T-mesh in the parameter space and the anchors of all blending functions defined over this mesh. Any interior regular node has exactly one function associated to it and the boundary nodes have more than one function associated due to the open knot vector structure. The process of inferring local knot vectors can be resumed as follows: •Blending functions are associated only to regular nodes of the mesh. •Local knot vectors are inferred by walking across the mesh until intersecting the mesh skeleton. This intersection should not coincide with a T-junction perpendicular to the marching direction. •Boundary knots are repeated to create an open knot vector structure along the boundary. Next, in order to span a spline space with good properties, some function supports should be modified. This issue is addressed in the next subsection. 3.3. Modification of local knot vectors The key of our strategy lies in some simple rules used for the modification of the function supports that lead to the construction of a polynomial spline space over a given 0-balanced T-mesh. In order to describe the idea, let us introduce some notation. For the local knot vectors Ξj α= (ξj 1, ξj 2, ξj 3, ξj 4, ξj 5), j = 1, ..., d let us denote the length of each knot interval as ∆j i=ξj i+1 −ξj i, j = 1, ..., d and i= 1, ..., 4. The support of a d-variate blending function Nαis a d-dimensional rectangular box: [ξ1 1, ξ1 5]× · · · × [ξd 1, ξd 5]. We are going to call frame of a function support the union of all (d−2)-faces of this box and we will denote it by frm(supp Nα). That is, for the rectangular support of a bivariate function, the frame is the union of the four vertices of this rectangle. For the cuboidal support of a trivariate function, its frame is composed of the union of the twelve edges of this cuboid. 9 Algorithm 2: Extension rule 2 in 2D. Input: A 0-balanced mesh Tand a pair of local knot vectors S={Ξ,H}. 1Function Modify2(T, S) 2S∗←S 3h= max(∆ξ 2, ∆ξ 3) 4for n∈ {1,5}do 5for m∈ {1,5}do 6if (ξn, ηm)/∈skt(T)then 7ξ∗ n←ξ3+ 3hsgn(ξn−ξ3) 8η∗ m←η3+ 3hsgn(ηm−η3) 9return S∗ Output: A modified support S∗={Ξ∗,H∗}that satisfies Condition 2. 4.2. Support extension for 3D meshes In this section we give a description and illustration of the proposed strategy for defining trivariate spline functions over 0-balanced octree T-meshes. The skeleton skt(T) of a three-dimensional mesh Tis the union of all faces of the mesh. For a trivariate function let us denote the vertices of its support as Vn,m,k = (ξn, ηm, ζk) where n, m, k ∈ {1,5}. And the edge of a support formed by the two vertices Vn,m,k and Vp,q,r is denoted as E(n,m,k),(p,q,r). Then, the frame frm(supp Nα) of a trivariate function support is the union of its twelve edges. Figure 9 illustrates the introduced notation for the support of a trivariate blending function. The formulation of Condition 1 for the local knot vectors of a trivariate function is analogue to the 2D case. Condition 2 adapted to 3D meshes is stated as follows: Edges of the cuboidal function support should be situated over the mesh faces. The implementation of the strategy for 3D is similar to the 2D case. To satisfy Condition 1, Extension rule 1 is applied to each of the three local knot vectors of a function analogously to the 2D case using Algorithm 1. Extension rule 2. In order to fulfill Condition 2, we check whether the edges of a function support are situated over the mesh faces. If not, the two knot vectors perpendicular to this edge should be modified by skipping over a knot for both parametric directions and placing this edge over the mesh faces. This checking is performed independently for each of the eight quadrants of the function support and, in each quadrant, three edges should 16 Figure 9: Support notation for 3D case. be checked. Figure 10 illustrates the support extension procedure for the quadrant of the vertex V5,1,1. Due to the octree structure only two cases can take place: (i) the quadrant contains one edge that is not situated over the mesh faces or (ii) the quadrant contains three edges and a vertex that are not situated over the mesh skeleton. Now we study each case. (i) If a quadrant contains one edge that does not fulfill Condition 2, then two knot vectors perpendicular to this edge are modified, see Fig. 10(a). For the function support shown in Fig. 10(a) left, the edge E(5,1,1),(5,1,5) is not situated over the mesh faces. Therefore, its two knot vectors Ξ and Hperpendicular to this edge are modified in order to place the edge over the mesh faces, namely, knots ξ5and η1are redefined as ξ∗ 5←ξ3+ 3hand η∗ 1←η3−3h, where h= max(∆ξ 2, ∆ξ 3) = max(∆η 2, ∆η 3) = max(∆ζ 2, ∆ζ 3). (ii) If a quadrant contains three edges that are not situated over the mesh faces, then the three knot vectors are modified by skipping over a knot for each of the three parametric directions, see Fig. 10(b). Vertex V5,1,1and the three edges connected to it are not situated over the mesh skeleton, so all the three knot vectors are modified to place the three edges over the mesh faces: ξ∗ 5←ξ3+ 3h,η∗ 1←η3−3h, ζ∗ 1←ζ3−3h. Algorithm 3 explains the Extension rule 2 used to modify a trivariate function support according to Condition 2. 17 (a) Only one edge violating Condition 2. Node V5,1,1is sited in the center of the face of size 2h. The support is extended in two directions. (b) Three edges violating Condition 2. Node V5,1,1is sited in the center of the cell of size 2h. The support is extended in three directions. Figure 10: Extension rule 2 for support modification of a trivariate function. 18 Algorithm 3: Extension rule 2 in 3D. Input: A 0-balanced T-mesh Tand three local knot vectors S={Ξ,H,Z}. 1Function Modify2 (T, S) 2S∗←S 3h= max(∆ξ 2, ∆ξ 3) 4mov(i) := i+ 4 sgn(3 −i) 5for n∈ {1,5}do 6for m∈ {1,5}do 7for k∈ {1,5}do 8if E(n,m,k),(mov(n),m,k)/∈skt(T)then 9η∗ m←η3+ 3hsgn(ηm−η3) 10 ζ∗ k←ζ3+ 3hsgn(ζk−ζ3) 11 if E(n,m,k),(n,mov(m),k)/∈skt(T)then 12 ξ∗ n←ξ3+ 3hsgn(ξn−ξ3) 13 ζ∗ k←ζ3+ 3hsgn(ζk−ζ3) 14 if E(n,m,k),(n,m,mov(k)) /∈skt(T)then 15 ξ∗ n←ξ3+ 3hsgn(ξn−ξ3) 16 η∗ m←η3+ 3hsgn(ηm−η3) 17 return S∗ Output: A modified support S∗={Ξ∗,H∗,Z∗}that satisfies Condition 2. 19 5. Properties of our spline functions, support overlapping, sparsity and condition number Here, we summarize the properties and some characteristics of the spline spaces constructed by means of our method. 5.1. Properties For any 0-balanced mesh T, the set of blending functions defined according to our strategy spans the space ST= span {Nα:α∈AT}with the following properties: 1. Functions {Nα}α∈ATare C2-continuous. 2. Functions {Nα}α∈ATare linearly independent. 3. Spaces spanned by nested T-meshes are also nested: T1⊂T2⇒ST1⊂ST2. 4. Polynomial reproduction property: P3(Ω) ∈ST. A rigorous proof of these properties is currently under preparation and we plan to give it in a future work. However, for a better understanding of the strategy, some clues should be mentioned. For the proof of linear independence of our blending functions it is essential to take into account that we do not repeat knots in the interior of the domain, and each regular node of the T-mesh has only one function assigned to it. Some reasoning from [18] can be used to proof the linear independence of such functions. A brief outline of the proof of the property (3) is given in the appendix of the present paper. It is worth mentioning that we have carried out numerous numerical experiments and the claimed properties were verified in all of them. 5.2. Support overlapping It should be pointed out that the proposed strategy can lead, in some cases, to increased overlapping of the function supports of different refinement levels, which can affect the conditioning and sparsity of the stiffness matrix. This effect can take place in some problems with a very sharp singularity when the area marked to refine by the error indicator at each iteration is smaller than the area occupied by the supports of the previous level. According to our strategy, a cell refinement always adds at least one new blending 20 (a) A function support (red) that intersects the refined cell (shaded cell). (b) New function support after applying additional refinements. (c) Three refinement steps around a singular point (black square) that lead to an excessive overlapping of the function supports. Note that the red support does not change. (d) The same three refinement steps applying additional refinements. Figure 11: Additional refinements in order to avoid excessive overlapping of the function supports. 21 function, while a function of the coarse level is replaced by the finer ones only if all the cells that share their anchor are refined. To avoid a possible accumulation of the functions of different refinement levels, the following strategy can be adopted: if a function support intersects a cell marked to be refined, then the coarsest cells sharing the anchor of the function are refined, see Fig. 11(a) and (b). Application of this strategy is illustrated in Fig. 11(c) and (d). In this example, let us suppose that the problem has a singularity at the center of the domain, and at each iteration, we refine only the four cells adjacent to the singularity point. Then, at each iteration there are functions whose support is not changed, see Fig. 11(c). However, if at each iteration we perform additional refinements, functions of the coarser level are replaced by the finer ones, see Fig. 11(d). This approach is similar to the one used in hierarchical refinement where all supports that intersect the marked cell are also refined in order to guarantee the replacing of the supports of the previous level by the finer ones at each refinement step. Application of this additional refinement, where needed, can reduce an excessive function overlapping and improve the stiffness matrix conditioning and sparsity. However, due to the unnecessary extension of the refined zone, the optimal rate of convergence can be lost. In practice, we did not observe a significant advantage of this approach, and taking into account its computational cost, we do not apply this strategy in the computational examples presented in this paper. Moreover, it was observed that for a large variety of problems, the excessive support overlapping is avoided naturally after various refinement steps due to the behaviour of the error indicator. On the other hand, another possible solution for this problem is a more accurate and selective definition of the function supports in order to obtain a space with better support locality. The strategy we propose in this paper can be improved by including more sophisticated rules for inferring function supports, but that would probably lead to a more expensive implementation. 6. Computational examples In this section we present computational examples of the application of our technique in geometric design and analysis for problems involving adaptive refinement. The proposed method is tested by performing surface interpolation and resolution of 2D and 3D Poisson problems using IGA. 22 (a) (b) (c) Figure 12: Igea’s face. Spline representation of the surface from its triangulation. (a) Input surface triangulation; (b) spline representation of the surface; (c) parametric T-mesh. 6.1. Spline representation of the surface from its triangulation Here we construct a spline representation of a surface given by a triangulation. First, a global parameterization of the surface triangulation is obtained by means of the method proposed by M. Floater in [19]. As a result, we have one-to-one mapping from planar triangulation of the parametric domain (unit square) to the surface triangulation. Then, a quadtree T-mesh adapted to the planar triangulation is constructed. For this purpose we have chosen the following simple criterion. We start from a coarse T-mesh and refine it until each cell of the mesh contains no more than a certain number of points of the input triangulation (3 in our case). The spline approximation of the surface is built as a linear combination of our blending functions S(ξ) = ∑ α∈AT CαNα(ξ). The control points Cαare found by imposing the interpolation conditions xβ=∑ α∈AT CαNα(ξβ),∀ξβ, β ∈AT, 23 where ξβare interpolation points in the parametric domain and xβare their images in the physical space determined by the triangular parameterization. As interpolation points we use the anchors of the functions, i.e., the regular nodes of the T-mesh, and some additional interpolation points associated to functions whose local knot vectors contain repeated knots on the boundary. These additional points are situated at the midpoint of the edges that have contact with the boundary and coincide approximately with the maximum point of the corresponding blending functions. See [20] for more details. The resulting spline surface, parametric T-mesh and input surface triangulation are shown in Fig. 12. 6.2. Adaptive refinement for surface interpolation In this example we interpolate the function u(r) = r1 2,(3) defined on the square domain [0,1]2, being r=√(x−0.5)2+ (y−0.5)2. The spline approximation of the surface (3) is built as a linear combination of our blending functions uh(ξ) = ∑ α∈AT CαNα(ξ). And the control points Cαare found by imposing the interpolation conditions u(ξβ) = ∑ α∈AT CαNα(ξβ),∀ξβ, β ∈AT, Adaptive refinement is performed according to the indicator based on exact L2interpolation error: η(Ωe)2=∥u−uh∥2 L2,Ωe=∫Ωe (u−uh)2dΩ. A cell Ωeis marked to be refined if η(Ωe)> γ maxi{η(Ωi)}, being γ∈ [0,1]. Some steps of adaptive refinement for the surface interpolation (3) are shown in Fig. 13, and Fig. 14 illustrates the convergence in L2-norm and H1-seminorm. 24 (a) (b) (c) (d) (e) Figure 13: Adaptive refinement for surface interpolation (3). (a) Function to interpolate. (b) Initial mesh. (c), (d) and (e) Some steps of adaptive refinement. 25 à H1error æ L2error æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à 3103 104 3104 101 100 10-1 10-2 10-3 10-4 10-5 10-6 number of degrees of freedom error 1 2 Figure 19: Convergence of L2-norm and H1-seminorm error for the Poisson problem of the section 6.4, for γ= 0.3. to Ω, i.e S:b Ω→Ω, using the algorithm described in our previous work [20]. This technique, based on a T-mesh untangling and optimization procedure, allows us to obtain a good quality parameterization suitable for application of IGA. The procedure is an extension of the ideas presented in our works [22, 23, 24, 25]. Another method for parameterization of computational domain, based on pillowing technique, was given in [26, 27]. The mean ratio Jacobian is used to evaluate the quality of the parameterization in the sense of its orthogonality and uniformity. Figures 17(a) and (b) show the colormap of the parameterization quality for our parametric and physical domains. It should be mentioned that for all numerical examples the parameterization of the computational domain is performed using the same blending functions that are used for the solution approximation, so isoparametric concepts holds during isogeometric analysis. That is, if VT(b Ω) is the finite dimensional space spanned by the spline functions associated to a parametric mesh T, then the parametric mapping Sis constructed as a linear combination of functions from VT(b Ω) and the discrete approximation space VT(Ω) in physical domain is defined as follows: VT(Ω) = {v∈H1(Ω) : v=bv◦S−1,for all bv∈VT(b Ω)}. 32 Now we consider a Poisson problem with Dirichlet boundary condition on Ω, which exact solution is a function with steep wave front given by u(r) = arctan(α(r−r0)), where r=√(x−xc)2+ (y−yc)2, the parameter αdetermines the steepness of the wave front and r0is its location. In this example α= 200 and r0= 0.6.The center of the wave front (xc, yc) = (0,0) is situated outside our computational domain, so the function is smooth in Ω. The numerical solution of the problem and the mesh corresponding to the final refinement iteration is shown in Fig. 18. As expected, the error estimator has marked for refinement the zone of the wave front. The evolution of the exact error in L2-norm and H1-seminorm are shown in Fig. 19. 6.5. Poisson problem on 3D domain The next computational example is the resolution of a 3D Poisson problem using IGA. The computational domain is a spline approximation of a sphere portion, see Fig. 20. The initial uniform mesh was composed by 4 ×4×4 cells. The approximation is constructed using our spline blending function. The Poisson problem with Dirichlet boundary condition is set up so that the analytical solution of the problem is u(r) = sin (1 α+r), where r=√x2+y2and parameter α= 1/10π. This is a smooth function with an oscillation near the origin. Results of the final refinement iteration are shown in Fig. 20. The convergence behaviour of the adaptive refinement is illustrated in Fig. 21. 7. Conclusions and future research In this paper we have proposed a strategy for defining C2-continuous cubic tensor product spline functions over quadtree and octree T-meshes. We only demand these T-meshes to be 0-balanced and this requirement can be easily satisfied by using a standard balancing procedure. The proposed strategy includes simple instructions used for inferring local knot vectors to define blending functions. We conjecture that the resulting spline spaces have nice properties: linear independence and the characteristic that spaces spanned 33 (a) (b) (c) (d) (e) (f) (g) (h) Figure 20: Mesh and numerical solution in the sixth adaptive refinement step for the 3D Poisson problem of the section 6.5. (a) Parametric mesh. (b) Physical mesh. (c) Numerical solution in the parametric domain. (d) Numerical solution in the physical domain. (e) A section of the parametric mesh. (f) A section of the physical mesh. (g) Numerical solution in a section of the parametric domain. (h) Numerical solution in a section of the physical domain. 34 æ æ æ æ æ æ æ æ æ æ æ æ æ æ à à à à à à à à à à à à à à 102 103 104 105 101 100 10-1 10-2 10-3 10-4 10-5 number of degrees of freedom error à H1error æ L2error 1.5 1 Figure 21: Convergence of L2-norm and H1-seminorm error for the Poisson problem of the section 6.5, for γ= 0.3. by nested T-meshes are also nested. The above mentioned properties were verified in numerous numerical experiments and we plan to provide a rigorous proof in our next work. The implementation of our technique is straightforward taking into account the balance mesh restriction. Examples of adaptive refinement using IGA for 2D and 3D Poisson problems have been presented. In all of them, optimal rates of convergence have been obtained. We believe that the simplicity of our technique can make it an attractive tool for its application in IGA and geometric design. In future works we also plan to improve the strategy obtaining a better locality of the function supports. Acknowledgements This work has been supported by the Spanish Government, “Secretar´ıa de Estado de Universidades e Investigaci´on,” “Ministerio de Econom´ıa y Competitividad,”, “Programa de FPU del Ministerio de Educaci´on, Cultura y Deporte”, “Programa de FPI propio de la ULPGC” and FEDER, grant contracts: CGL2011-29396-C03-01 and CGL2011-29396-C03-03; “Junta Castilla Le´on,” grant contract: SA266A12-2. It has also been supported 35 Figure 22: Elemental refinements. A possible sequence of 0-balanced meshes {Te i}, where Te iis obtained from Te i−1by refining only one cell. by CONACYT-SENER (“Fondo Sectorial CONACYT SENER HIDROCARBUROS,” grant contract: 163723). Appendix Here we include an outline of the proof (2D case) of the nestedness of the spline spaces constructed via our strategy. It is worth noting that the rules we propose and the 0-balanced restriction for the mesh lead to a very few number of all possible function supports that can be defined over a T-mesh. In order to proof the nestedness of the spaces for any two nested T-meshes T1⊂T2, we should verify that any blending function of ST1can be represented as linear combination of blending functions from ST2. To this end, we have to perform this verification for all types of function supports and for any possible mesh configuration. Applying some reasoning, it is easy to see that this verification can be reduced to a few number of basic cases. Firstly, to simplify the process, we assume that for any two 0-balanced meshes T1and T2such that T1⊂T2there exists a sequence of meshes {Te i}i=1,n such that (i) T1=Te 0⊂Te 1⊂Te 2⊂ · · · ⊂ Te n−1⊂Te n=T2, (ii) any mesh Te iis obtained from mesh Te i−1by refining only one cell, 36 (a) ... ... ... ... m m+1 m-1 m-2 (b) Figure 23: (a) All possible configurations for a knot vector (without repeated knots). (b) On the left, a function support of level mdefined over a 0-balanced T-mesh. The level of this function is determined by the cell marked in red. On the right, the quadtree representation of the T-mesh. The tree node represented in red corresponds to the marked cell of the mesh. (iii) all the meshes {Te i}i=1,n are 0-balanced. A pair of 0-balanced meshes Te 2and Te 1is called elemental refinement if Te 2is obtained from Te 1by refining only one cell, see Fig. 22. Since it is sufficient to verify the nestedness of spaces for any elemental refinement, from now on we assume that T1⊂T2is an elemental refinement. We have to verify that, after refining a cell, the blending functions of the new space ST2are able to represent any blending function of the previous space ST1. The proof can be divided in two steps. First, we analyze some basic cases corresponding to the simplest mesh configuration, then, more general mesh configuration should be considered applying some recursivity idea to proof the nestedness of spaces for arbitrary elemental refinement. To make the necessary studies more graphic, first, we classify all function supports that can take place according to our strategy. Then, we study for each type all possible situations when a cell refinement affects this type function support. Classification of function supports Let us consider any knot vector Ξ = (ξ1, ξ2, ξ3, ξ4, ξ5) and let ∆ξ i=ξi+1−ξi, i= 1, ..., 4., so that, a set ∆ξ= (∆ξ 1, ∆ξ 2, ∆ξ 3, ∆ξ 4) corresponds to each knot vector Ξ. Due to the extension rule 1 and the balanced mesh condition we have the following possible configurations for the set ∆ξ: (h, h, h, h), (h, h, h, 2h), (2h, h, h, h) and (2h, h, h, 2h), see Fig. 23(a). All possible bivariate function supports are obtained by combining all configurations for Ξ 37 (a) (b) Figure 24: Two possible scenarios for a function support according to the underlying local mesh. (a) The anchor of the m-level function coincides with the center of a (m−1)-level cell of the quadtree. (b) The rest of the functions. and H, i.e., combining two sets ∆ξand ∆η. Note that since we ignore the hanging nodes while inferring two knot vectors, the value of his the same for ∆ξand ∆η, i.e., h=∆ξ 2=∆ξ 3=∆η 2=∆η 3and the value of hdetermines the level of a function. A blending function is said to be of level mif the largest cell contacting its anchor is of level m, see Fig. 23(b). Also, we can distinguish different function supports according to the underlying local mesh that leads to this support.Basically, it depends on the position of the anchor of the function, see Fig. 24. Scenario (a) the anchor of the m-level function coincides with the center of a (m−1)-level cell of the quadtree, scenario (b) the rest of the functions, i.e., the anchor of the m-level function does not correspond to any center of a (m−1)-level cell of the quadtree. Finally, taking into account all the above considerations about values of ∆ξand ∆ηand the underlying local mesh, we can classify all function supports in nine types that are illustrated in Fig. 25. First six types correspond to scenario (a) and the other three types correspond to scenario (b). Note that type 6 and type 9 function supports can be considered as the same type, although they come from different scenarios. Function supports with repeated knots across the domain boundary can be included in the nine types shown in Fig. 25. Next, we have to see how an elemental refinement can affect a function support of each type. As it was said before, the proof of the nestedness can be reduced to the proof of nestedness for some basic cases. 38 (a) Type 1 (b) Type 2 (c) Type 3 (d) Type 4 (e) Type 5 (f) Type 6 (g) Type 7 (h) Type 8 (i) Type 9 Figure 25: Function support classification according to the underlying local mesh and configurations for ∆ξand ∆η. Function recuperation for the basic cases For each support type, we have to analyze all possible elemental refinements that affect this function and verify that the blending functions of the finer mesh are able to reproduce the function under consideration. Let see in detail an example of this study. Here we illustrate type 2 function support. This type of function support has the largest number of cases to study. Namely, we should study four different situations that take place after a cell refinement. In these situations the original type 2 function N0associated to the anchor v0is not a blending function of the new space ST2. The new function ˆ N0associated to v0can be 39 (a) N0=ˆ N0+1 4ˆ N1(b) N0=ˆ N0+1 4ˆ N1 (c) N0=ˆ N0+1 4ˆ N1+1 4ˆ N2+1 16 ˆ N3(d) N0=ˆ N0−1 4ˆ N1 Figure 26: Basic cases to study for type 2 function. Support of the function N0defined over T1is represented by pink shadowed area in the leftmost and the rightmost images. Blue discontinuous lines in the leftmost figure mark the refinement of T1that leads to a new mesh T2. Supports and the anchors of the blending from ST2necessary to recuperate the function N0are marked by different colors in the rightmost image. (a) Type 2 support turns into a type 3 support after a cell refinement. (b) Type 2 support turns into a type 4 support. (c) Type 2 support turns into a type 5 support after a cell refinement. (d) Type 2 support turns into a type 5 support. type 3, 4, 5 or 1, as shown below. Type 2 to type 3. In this case, type 2 function N0of level mcan be recuperated by two new functions of the same level: type 3 function ˆ N0and type 7 function ˆ N1. Only one of local knot vectors undergoes a knot insertion. Without loss of generality, suppose that a new knot ˆ ξis inserted between the fourth and the fifth knot of Ξ, see Fig. 26(a). Then, using the knot insertion formula, it is easy to see that N0(ξ, η) = B[Ξ]B[H] = B[Ξ1]B[H] + 1 4B[Ξ2]B[H] = ˆ N0+1 4ˆ N1, where Ξ1= (ξ1, ξ2, ξ3, ξ4,ˆ ξ) and Ξ2= (ξ2, ξ3, ξ4,ˆ ξ, ξ5). 40 6 ≡ 9 8 7 5 4 3 2 1 Figure 27: Graph of relations between different function types. An arrow coming out of a circle with type number and pointing to another type means that this type support can turn into the second one after a cell refinement. Type 2 to type 4. Function N0is split into two new functions of the same level: type 4 function ˆ N0and type 8 function ˆ N1. Local knot vector H undergoes a knot insertion. A new knot ˆηis inserted in its local knot vector H= (η1, η2, η3, η4, η5) between the fourth and the fifth knot, see Fig. 26(b). Then N0(ξ, η) = B[Ξ]B[H] = B[Ξ]B[H1] + 1 4B[Ξ]B[H2] = ˆ N0+1 4ˆ N1, where H1= (η1, η2, η3, η4,ˆη) and H2= (η2, η3, η4,ˆη, η5). Type 2 to type 5. For this case both local knot vectors undergo a knot insertion. New knot ˆ ξis inserted between the fourth and the fifth knot of Ξ and knot ˆηis inserted between the fourth and the fifth knot of H, see Fig. 26(c). Then function N0is split into four functions of the same level: type 5 function ˆ N0, type 8 function ˆ N1, type 9 function ˆ N2and type 9 (or 6) function ˆ N3. N0(ξ, η) = B[Ξ]B[H] = (B[Ξ1] + 1 4B[Ξ2])(B[H1] + 1 4B[H2]) = =B[Ξ1]B[H1] + 1 4B[Ξ2]B[H1] + 1 4B[Ξ1]B[H2] + 1 16B[Ξ2]B[H2] = =ˆ N0+1 4ˆ N1+1 4ˆ N2+1 16 ˆ N3. where H1= (η1, η2, η3, η4,ˆη), H2= (η2, η3, η4,ˆη, η5), Ξ1= (ξ1, ξ2, ξ3, ξ4,ˆ ξ) and Ξ2= (ξ2, ξ3, ξ4,ˆ ξ, ξ5). 41