Truncated hierarchical Catmull–Clark subdivision with local refinement

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20 www.elsevier.com/locate/cma

Truncated hierarchical Catmull–Clark subdivision with local refinement Xiaodong Wei a , Yongjie Zhang a,∗ , Thomas J.R. Hughes b , Michael A. Scott c a Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, USA b Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712, USA c Department of Civil and Environmental Engineering, Brigham Young University, Provo, UT 84602, USA

Received 4 September 2014; received in revised form 18 March 2015; accepted 24 March 2015 Available online 1 April 2015

Abstract In this paper we present a new method termed Truncated Hierarchical Catmull–Clark Subdivision (THCCS), which generalizes truncated hierarchical B-splines to control grids of arbitrary topology. THCCS basis functions satisfy partition of unity, are linearly independent, and are locally refinable. THCCS also preserves geometry during adaptive h-refinement and thus inherits the surface continuity of Catmull–Clark subdivision, namely C 2 -continuous everywhere except at the local region surrounding extraordinary nodes, where the surface continuity is C 1 . Adaptive isogeometric analysis is performed with THCCS basis functions on a benchmark problem with extraordinary nodes. Local refinement on complex surfaces is also studied to show potential wide application of the proposed method. c 2015 Elsevier B.V. All rights reserved. ⃝

Keywords: Catmull–Clark subdivision; Local refinement; Truncation mechanism; Hierarchical B-splines; Extraordinary nodes; Isogeometric analysis

1. Introduction Local refinement has been an important aspect of isogeometric analysis research [1,2] since the introduction of isogeometric analysis with NURBS (Non-Uniform Rational B-Spline) basis functions. Several different schemes have been developed that break the tensor product structure of NURBS in order to support local refinement, namely T-splines [3,4], PHT-splines [5], hierarchical B-splines [6,7], and locally refined splines [8,9]. T-splines have been widely used in geometric modeling [10–12] and isogeometric analysis [13,14]. However, the initial definition of T-splines proved too general in that linear independence [15,16] of the basis functions was not guaranteed. Linear independence is a requisite property in isogeometric analysis. Analysis suitable Tsplines [17] were subsequently developed to ensure this property. However, local refinement in analysis-suitable Tsplines may propagate beyond the domain of interest. Hierarchical B-splines were originally introduced by Forsey ∗ Corresponding author. Tel.: +1 412 268 5332; fax: +1 412 268 3348.

E-mail address: [email protected] (Y. Zhang). http://dx.doi.org/10.1016/j.cma.2015.03.019 c 2015 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

and Bartels [18], supporting local refinement of coarse B-spline patches through the overlapping of fine ones. Later Kraft [19] developed a global selection mechanism to select linearly independent basis functions from different B-spline hierarchical levels. Recently Giannelli et al. [20] developed a truncation mechanism for hierarchical Bsplines to satisfy partition of unity and to decrease the overlapping of basis functions from different levels for better numerical conditioning. However, a major disadvantage of hierarchical B-splines is its restriction to a global rectangular parametric domain. Complex topologies necessarily involve extraordinary nodes (an extraordinary node has other than four faces adjacent to it) and cannot be represented in terms of global rectangular parametric domains. As such, hierarchical B-splines cannot support arbitrary topology. Subdivision surfaces have been used broadly in free-form surface design and are advantageous in handling extraordinary nodes [21–23]. A thorough survey of different subdivision schemes was given by Zorin et al. [24] and Sabin [25]. Among them, Catmull–Clark subdivision [26] was developed based on bicubic B-splines. With the advent of explicit Catmull–Clark basis functions derived by Stam [27], hierarchical Catmull–Clark subdivision has been implemented in commercial products such as Pixar’s Renderman [28], where coarse patches are simply overlaid with fine ones and linear independence is not guaranteed. Before isogeometric analysis emerged, Cirak et al. [29,30] developed a unified framework to integrate geometric design and simulation with subdivision surfaces for thin shell structures. Catmull–Clark subdivision solids were recently introduced into isogeometric analysis [31]. However, adaptive simulation was not discussed in these studies. Grinspun et al. [32] developed adaptive simulation with basis refinement rather than traditional element refinement for subdivision surfaces. However, basis refinement may produce ill-shaped control mesh that is not suitable for geometric design [33]. In this paper we develop a new method, termed Truncated Hierarchical Catmull–Clark Subdivision (THCCS), to incorporate extraordinary nodes and support local refinement via a similar hierarchical structure as in truncated hierarchical B-splines. We introduce the local selection into hierarchical Catmull–Clark subdivision to ensure linear independence of basis functions. Instead of rationalizing the basis functions, we also introduce the truncation idea into subdivision basis functions to make them satisfy partition of unity. Partition of unity is essential to the convex hull property in geometric design. THCCS preserves the exact geometry when adaptive h-refinement is performed and it inherits the surface continuity of Catmull–Clark subdivision, namely C 2 -continuity everywhere except C 1 -continuity at extraordinary nodes. THCCS basis functions have several nice properties suitable for both geometric design and isogeometric analysis, including partition of unity, linear independence, convex hull at each hierarchical level, supporting local refinement and extraordinary nodes, as well as smooth surface continuity. We have applied our THCCS basis functions in isogeometric analysis on a 2D benchmark problem with extraordinary nodes and several complex geometries. The simulation results show potential wide application of the proposed method in integrating design and analysis. The remainder of this paper is organized as follows. Section 2 begins with the basic idea of hierarchical B-splines and the truncation mechanism. Section 3 introduces the explicit Catmull–Clark basis functions and investigates their refinability. In Section 4, the main idea of THCCS is discussed. Subsequently, in Section 5 we prove geometry preservation during local refinement. Section 6 describes a benchmark problem and several complex geometries with isogeometric analysis results, and Section 7 concludes the paper and suggests future work. 2. Review of hierarchical B-splines To facilitate the development, in the following we first review several essential concepts of hierarchical B-splines and the truncation mechanism. For related details we refer to related literature [18,19,6,7,20,24]. 2.1. Hierarchical B-splines Univariate B-spline basis functions Bi, p (ξ ) (i = 0, . . . , n − 1) of a given degree p are defined on a knot vector Ξ = {ξ0 , ξ1 , . . . , ξn+ p }, where ξi ∈ R is the ith knot and n is the number of basis functions. Bi, p (ξ ) has local support on [ξi , ξi+ p+1 ]. B-spline basis functions are refinable, which enables the construction of hierarchical Bsplines. Suppose Ξ 0 is an initial uniform knot vector, Ξ ℓ (ℓ = 0, 1, . . .) can be obtained by ℓ subdivisions of Ξ 0 . Refinability indicates that a basis function Bi,ℓ p defined on Ξ ℓ can be represented as a linear combination of p + 2 basis functions defined on Ξ ℓ+1 , Bi,ℓ p (ξ ) =

p+1  k=0

ℓ+1 ck B2i+k, p (ξ ) p

p

with ck =

  1 p+1 , 2p k

i = 0, 1, . . . , n ℓ − 1,

(1)

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

ℓ . (a) B0,3

ℓ+1 (b) Bk,3 (k = 0, . . . , 4).

3

ℓ+1 (c) ck3 Bk,3 (k = 0, . . . , 4).

ℓ is defined on the knot vector Ξ ℓ = {0, 1, 2, 3, 4}; (b) B ℓ+1 (k = 0, . . . , 4) Fig. 1. Refinability of a B-spline basis function of degree p = 3. (a) B0,3 k,3   ℓ+1 1 3 5 7 ℓ+1 are defined on Ξ = {0, 2 , 1, 2 , 2, 2 , 3, 2 , 4}; and (c) Bk,3 is weighted with ck3 = 13 k4 and k = 0, . . . , 4. Each color represents one basis 2 3 ℓ+1 ℓ = 4 function, and we have B0,3 k=0 ck Bk,3 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

where ck are the refinement coefficients and n ℓ is the number of basis functions defined on Ξ ℓ . The p + 2 basis ℓ+1 ℓ functions B2i+k, p are called the children of Bi, p , denoted as p

ℓ+1 chd Bi,ℓ p = {B2i+k, p | k = 0, 1, . . . , p + 1}.

(2)

ℓ and its children B ℓ+1 , where k = 0, . . . , 4. Fig. 1 shows a uniform univariate cubic B-spline basis function B0,3 k,3 ℓ According to Eq. (1), B0,3 can be represented by a weighted summation of its five children. Given two knot vectors Ξ = {ξ0 , ξ1 , . . . , ξn+ p } and H = {η0 , η1 , . . . , ηm+q }, and polynomial degrees p and q respectively, a bivariate basis function Bi,p (ξ, η) is defined by a tensor product of two univariate basis functions,

Bi,p (ξ, η) = Bi, p (ξ )B j,q (η),

(3)

where i = (i, j) and p = ( p, q) are two dimensional indices and degrees. Bi, p (ξ ) and B j,q (η) are defined on Ξ and H, respectively. Similar to Eq. (2), the children of Bi,p (ξ, η) are denoted as chdBi,p (ξ, η) with ( p + 2) × (q + 2) basis functions. Bi,p (ξ, η) has local support, supp Bi,p (ξ, η) = [ξi , ξi+ p+1 ] × [η j , η j+q+1 ].

(4)

The bivariate B-spline basis is B = {Bi,p | i = 0, . . . , n − 1; j = 0, . . . , m − 1}.

(5)

With the above notations, we then construct bivariate hierarchical B-splines. Starting from the initial parametric domain Ω 0 with equally spaced knots, we define a set of B-spline basis functions B 0 . The union of support of all the basis functions in B 0 is Ω 0 at Level 0, Ω 0 = supp B 0 . Assume that the degree is fixed during the construction. According to Kraft [19], the function space spanned by B 0 can be enlarged by replacing the selected B-spline basis functions with their children. Such replacement indicates a local refinement of basis functions. The refinement will be recursively performed until the maximum level, ℓmax . In the following, we only study two consecutive levels and construct Level-(ℓ + 1) basis functions from Level ℓ, where ℓ ≥ 0. Fig. 2 illustrates the construction process of univariate cubic hierarchical B-spline basis functions in three steps: • Step 1. Identify a set of basis functions Brℓ ⊆ B ℓ to be refined at Level ℓ (blue dashed curve) whose supports define the subdomain Ω ℓ+1 = supp Brℓ at Level ℓ + 1, and designate Brℓ as passive and the remaining basis functions in B ℓ (blue solid curves) as active, Baℓ = B ℓ \Brℓ . • Step 2. Obtain all the children of Brℓ at Level ℓ + 1 (orange solid curves) and define them as active; we have B ℓ+1 = Baℓ+1 = chd Brℓ . • Step 3. Collect all the active basis functions at Levels ℓ and ℓ + 1 (blue and orange solid curves) to obtain the hierarchical B-spline basis functions, ℓ+1 Bhb = Baℓ ∪ Baℓ+1 .

(6)

4

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

(a) Step 1 (no truncation).

(d) Step 1 (with truncation).

(b) Step 2 (no truncation).

(e) Step 2 (with truncation).

(c) Step 3 (no truncation).

(f) Step 3 (with truncation).

Fig. 2. Three steps to construct univariate cubic hierarchical B-spline basis functions without truncation (a)–(c) and with truncation (d)–(f). (a), (d): In Step 1 (Level ℓ), the blue dashed curves are identified as basis functions Brℓ to be refined and they are designated as passive, all the other basis functions (blue and red solid curves) are active (Baℓ ). Red curves are truncated basis functions; (b), (e): In Step 2 (Level ℓ + 1), the orange solid curves represent the five children of Brℓ and they are designated as active (Baℓ+1 ), while the orange dashed curves are passive; and (c), (f): In Step 3, all the active basis functions from Levels ℓ and ℓ+1 (blue, orange and red solid curves) are collected to form the hierarchical B-spline ℓ+1 basis functions Bhb . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Eq. (6) refers to the so-called global selection of active basis functions, and the above recursive construction terminates at Level ℓmax − 1. The active basis functions are updated in each recursive step. Finally all the active basis functions from different levels form the basis of hierarchical B-splines. As proved by Kraft [19], the hierarchical B-spline basis functions are globally linearly independent. Generally, the hierarchical B-spline basis functions in non-rational form do not satisfy partition of unity. Moreover, local refinement of hierarchical B-splines focuses on basis refinement, which may produce ill-shaped control meshes at the refined level. To overcome this deficiency, a truncated mechanism for hierarchical B-splines [20] was developed. 2.2. Truncated basis functions The truncated mechanism was first developed by Giannelli et al. [20] for the hierarchical B-spline basis functions to form a partition of unity and to decrease the overlapping of basis functions for better numerical conditioning. Very similar to hierarchical B-splines, the construction of truncated hierarchical B-splines also consists of three steps as shown in Fig. 2(d)–(f). The only difference is that in Step 1 the remaining active basis functions (Baℓ ) surrounding Brℓ need to be modified or truncated; this is illustrated in red curves in (d). All other steps remain the same. In the following, we will focus on the construction of a truncated basis function. Definition. Given a set of basis functions Brℓ to be refined, we define the subdomain Ω ℓ+1 = suppBrℓ . Provided that Biℓ ̸∈ Brℓ is refinable, according to the refinability in Eq. (1) we have  Biℓ = ci,j Bjℓ+1 , (7) suppBjℓ+1 ⊆suppBiℓ

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

5

where ci,j ∈ R are the refinement coefficients from mid-knot insertions [7], and Bjℓ+1 ∈ chdBiℓ . The truncated basis function of Biℓ with respect to Brℓ is defined as  ci,j Bjℓ+1 . (8) trunBiℓ = suppBjℓ+1 ̸⊆Ω ℓ+1

Eq. (8) indicates that among all the children of Biℓ , we discard those with support fully contained in Ω ℓ+1 when we construct the truncated basis function trunBiℓ . As shown in Fig. 2(d), the blue dashed curve is still the basis function Brℓ to be refined, and we have Ω ℓ+1 = [3, 7]. For univariate cubic hierarchical B-splines, each basis function at Level ℓ has five children at Level ℓ + 1, and only the four basis functions surrounding Brℓ need to be truncated because they have children with supports fully contained in Ω ℓ+1 . For the two basis functions adjacent to Brℓ , we discard three children. And for the other two basis functions, we only need to discard one child. Far away basis functions have no children within Ω ℓ+1 , so we do not need to truncate them. These truncated basis functions (red curves) are still designated as active, and they are collected to form the truncated hierarchical B-spline basis functions (Step 3) together with the active basis functions from Step 2 (Baℓ+1 = chd Brℓ ); see Fig. 2(e)–(f). The hierarchical B-spline basis with truncation has been proven to form a partition of unity and therefore achieves strong stability [34]. It gives a sparser connectivity among basis functions at different levels, and it can preserve geometry when local refinement is performed. 3. Catmull–Clark basis functions and their refinability Our THCCS scheme uses Catmull–Clark basis functions to handle extraordinary nodes, so we study them here and investigate their refinability. 3.1. Explicit Catmull–Clark basis functions By virtue of the explicit basis functions [27], a Catmull–Clark surface can be constructed from a valid control mesh. A valid control mesh requires that it only consists of quadrilaterals and each quadrilateral contains at most one extraordinary point. The appearance of extraordinary points makes it impossible to employ a global rectangular parameter domain. Instead, we need to use a local parametric domain to each quadrilateral. Such local parametric domain is called an element. There are two kinds of elements in the Catmull–Clark subdivision: regular elements and irregular elements. An element is regular if it does not contain an extraordinary node,1 otherwise it is irregular. Each node is associated with a basis function, which is derived from uniform bi-cubic B-spline basis functions. Thus a Catmull–Clark basis function has local support on its two-ring neighborhood elements. A patch is then defined for an element by a set of basis functions that have support on it, together with the corresponding control points. A patch corresponding to a regular element is called a regular patch, otherwise it is irregular. As shown in Fig. 3(a), a regular patch contains 16 nodes whose basis functions have support over the element marked in blue. These functions are uniform B-spline basis functions, and the element is a unit square parametric domain. Fig. 3(b) shows an irregular patch corresponding to Ω0ℓ , where the indexed 2N + 8 basis functions have support (N is the valence number of the extraordinary node). To evaluate an irregular patch, we subdivide it into an arbitrarily small irregular patch and an infinite number of regular patches. In Fig. 3(c), Ω0ℓ is subdivided into one irregular patch Ω0ℓ+1 and three regular patches Ωkℓ+1 (k = 1, 2, 3), with 2N + 17 basis functions having support over them. The four ¯ (2N +17)×(2N +8) . In particular, the patches Ωkℓ+1 (k = 0, 1, 2, 3) can be obtained from Ω0ℓ via a subdivision matrix A ℓ+1 ¯ irregular patch Ω0 can be obtained from the first 2N + 8 rows of A, denoted as the square matrix A(2N +8)×(2N +8) . For computational efficiency, the eigenstructure (3, V) of the subdivision matrix A (AV = V3) [27] is employed and only diagonal matrix multiplication is needed. In this manner, the basis functions at Level ℓ over an irregular element can be derived explicitly by ¯ T b(ξ, η), Bℓ (ξ, η) = (V−1 )T 3n−ℓ−1 (Pk AV) 1 A node is the correspondence of a control point in the parametric domain.

(9)

6

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

Fig. 3. (a) A regular patch; (b) an irregular patch (corresponding to Ω0ℓ ) with 2N + 8 nodes where Node 1 is a valence-3 extraordinary node (N = 3); and (c) one irregular patch (Ω0ℓ+1 ) and three regular patches (Ωkℓ+1 , k = 1, 2, 3) after subdividing the irregular patch in (b), where

Ω0ℓ = ∪3k=0 Ωkℓ+1 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

where rows of Bℓ are the resulting basis functions, and rows of b(ξ, η) are the uniform cubic B-spline basis functions. Pk (k = 1, 2, 3) are selection matrices to locate which regular patch (Ωkn , where k = 1, 2, 3) certain parametric ¯ is the values (ξ, η) belong to after n times of subdivisions. The configuration of the subdivision matrix A and A same as Stam’s [27]. As (ξ, η) approaches arbitrarily close to (0, 0), n will become infinite and only the first diagonal element is nonzero in 3n−ℓ−1 . This is because all the eigenvalues are non-negative and smaller than 1 except the first eigenvalue which equals 1. Therefore in the limit, as n goes to infinity, we obtain Bℓ with only its first 2N + 1 basis functions nonzero due to the configuration of V−1 ; see Eq. (12) of Stam’s development [27]. These nonzero basis functions are the one-ring neighboring basis functions of the extraordinary node. In isogeometric analysis, we only need to evaluate basis functions at specific quadrature points which typically are a bit away from extraordinary nodes, and in practice only a few levels of subdivision are needed. 3.2. Refinability of Catmull–Clark basis functions Refinability of basis functions is fundamental to local refinement of hierarchical B-splines and the development of truncated basis functions, and THCCS development as well. We now study how coarse-level basis functions can be represented by fine-level basis functions in each patch for Catmull–Clark basis functions. It is trivial for regular patches where only B-spline basis functions are involved, so here we focus on the basis functions on irregular patches, that is, Ω0ℓ and Ω0ℓ+1 in Fig. 3(b), (c). Denoting the basis functions over Ω0ℓ and Ω0ℓ+1 as Bℓ and Bℓ+1 , respectively, and we have ¯ Tb Bℓ+1 = (V−1 )T 3n−(ℓ+1)−1 (Pk AV) ¯ Tb = (VT )−1 3−1 3n−ℓ−1 (Pk AV) ¯ T b. = (3VT )−1 3n−ℓ−1 (Pk AV)

(10)

Recall that AV = V3 and 3 is a diagonal matrix, so 3 = 3T and 3VT = VT AT .

(11)

By plugging Eq. (11) into Eq. (10), we obtain ¯ Tb Bℓ+1 = (VT AT )−1 3n−ℓ−1 (Pk AV) ¯ Tb = (AT )−1 (VT )−1 3n−ℓ−1 (Pk AV) = (AT )−1 Bℓ .

(12)

Then Bℓ = AT Bℓ+1 .

(13)

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

7

Eq. (13) indicates a linear relationship between basis functions over Ω0ℓ and Ω0ℓ+1 . We can observe that refinability holds for the Catmull–Clark basis functions via the subdivision matrix A, which is the basis of our THCCS construction. 4. Construction of THCCS With truncated hierarchical B-splines and the refinability relationship of Catmull–Clark basis functions, we now discuss how to develop THCCS. Instead of regular control meshes in hierarchical B-splines, irregular control meshes are the main focus in the construction of THCCS. Here we adopt the same notations as in Sections 2 and 3. Starting with an initial valid Catmull–Clark control mesh at Level 0, we define all its basis functions in a set, B 0 . The Level-0 domain Ω 0 is defined as the union of all the elements. It is also the union of the support of all basis functions in B 0 , 0

Ω =

e −1 n

Ωi0 = supp B 0 ,

(14)

i=0

where n e is the number of elements and Ωi0 denotes the ith element at Level 0. Recall that the support of a Catmull– Clark basis function is its two-ring neighborhood elements. In addition, we denote all Level-0 elements as the set E 0 = {Ωi0 | i = 0, 1, . . . , n e − 1}.

(15)

We then recursively construct THCCS in a similar manner as for truncated hierarchical B-splines. The main challenge is how to handle irregular control mesh in each step. 4.1. Three steps in recursive construction In the following, we discuss the construction process in three main steps: identification and truncation (Step 1), refinement (Step 2) and collection (Step 3). Step 1 aims at identifying to-be-refined basis functions and elements, defining active basis functions and constructing truncated basis function; Step 2 aims at refining the identified elements and defining active basis functions at the finer level; and finally, Step 3 collects all the active basis functions and elements. Step 1 (Identification and Truncation). The recursive construction allows us to consider two consecutive levels at one time: Level ℓ (ℓ ≥ 0) with basis functions B ℓ and Level ℓ + 1. First of all, we identify a set of basis functions to be refined at Level ℓ, denoted as Brℓ , by comparing the local geometry or simulation error with a given threshold. Then, we set all the to-be-refined basis functions (Brℓ ) to be passive, and the remaining basis functions to be active (Baℓ = B ℓ \Brℓ ). The fine-level domain Ω ℓ+1 is defined by the support of all the to-be-refined basis functions, Ω ℓ+1 = supp Brℓ . We can also define a set of to-be-refined elements Erℓ using all the Level-ℓ elements inside Ω ℓ+1 , Erℓ = {Ωiℓ | Ωiℓ ∈ E ℓ , Ωiℓ ⊆ Ω ℓ+1 },

(16)

where E ℓ contains all the elements at Level ℓ and Ωiℓ denotes a level-ℓ element. Due to the centrosymmetric property with respect to the extraordinary node, there only exist three types of special Catmull–Clark basis functions which can be identified as to-be-refined basis functions. The three types are marked in red, blue and green squares in Fig. 4(a), and the other unmarked basis functions are standard uniform cubic B-spline basis functions. The supports of these special basis functions are marked in pink in Fig. 4(b)–(d), denoting their tworing neighborhood elements. If the marked basis function in Fig. 4(b)–(d) is identified as to-be-refined (Brℓ ), all the elements inside the pink region will be identified as to-be-refined elements (Erℓ ). In this situation the pink region will also serve as the fine-level domain Ω ℓ+1 . With Brℓ , we also need to identify to-be-truncated basis functions (Btℓ ) from the remaining active basis functions ℓ Ba . If a basis function Biℓ ∈ Baℓ has any children with support fully contained in Ω ℓ+1 (or Biℓ and Brℓ have any shared children), then it is to be truncated. We have Btℓ = {Biℓ | Biℓ ∈ Baℓ , chd Biℓ ∩ chd Brℓ ̸= ∅}.

(17)

For Catmull–Clark basis functions, all the two-ring neighboring basis functions around Brℓ should be truncated, and they are marked in green circles as shown in Fig. 4(b)–(d). We truncate each basis function Biℓ ∈ Btℓ by discarding the

8

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

Fig. 4. (a) Three types of special Catmull–Clark basis functions around a valence-3 extraordinary node are marked in red (type-1), green (type-2) and blue (type-3), respectively; and (b)–(d) the support of type-1, type-2 and type-3 basis functions (pink region). Green circles mark the basis functions to be truncated. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

children with support fully contained in Ω ℓ+1 , and obtain  trunBiℓ = ci j B ℓ+1 j ,

(18)

ℓ+1 supp B ℓ+1 j ̸⊆Ω

¯ according to Eq. (13). Eq. (18) is different from Eq. (8) in the basis functions involved where ci j comes from A (or A) and the corresponding coefficients. A special Catmull–Clark basis function has more children when the associated extraordinary node has a larger valence number. Moreover, the coefficients ci j here have values which are not all the same as those from mid-knot insertion; see the example in Section 4.2 for details. Step 2 (Refinement). This step aims at refining all the elements in Erℓ and obtaining all the children of Brℓ . The tobe-refined elements Erℓ may consist of regular and irregular elements. They are refined in terms of the corresponding patches. The refinement of a regular patch is straightforward where only B-spline basis functions are involved. Without loss of generality, here we consider local refinement of an irregular patch corresponding to an irregular element Ω0ℓ , marked in blue in Fig. 3(b). The 2N +8 control points of the patch are locally labeled and they are written in the matrix ℓ+1 ℓ , . . . , Pℓ ℓ T form, Pℓe = [Pe,1 (k = 0, 1, 2, 3) with 2N + 17 e,2N +8 ] . Ω0 is subdivided into four smaller elements Ωk ℓ+1 ℓ+1 T ℓ+1 ℓ+1 can be obtained by control points (Pe ) as shown in Fig. 3(c). We have Pe = [Pe,1 , . . . , Pe,2N +17 ] and Pℓ+1 e ¯ using the subdivision matrix A, ¯ ℓe . Pℓ+1 = AP e

(19)

Recall that the evaluation of these four patches corresponding to elements Ωkℓ+1 (k = 0, 1, 2, 3) gives the same surface as the evaluation of the original patch corresponding to Ω0ℓ . After the refinement, Ω0ℓ is defined as passive while the fine-level elements Ωkℓ+1 (k = 0, 1, 2, 3) are defined as active in Eaℓ+1 . Then the Level-(ℓ + 1) basis functions in these four patches can be locally derived as in Fig. 3(b), (c). Among them, we define all the children basis functions of Brℓ as active, Baℓ+1 = chd Brℓ . The active element set at Level ℓ is Eaℓ = E ℓ \Erℓ , and the active element set at Level ℓ + 1 is Eaℓ+1 = {Ωkℓ | k = 0, 1, 2, 3}.

Step 3 (Collection). In this step, we collect all the active basis functions from Levels ℓ and ℓ + 1 to form the THCCS basis functions at Level ℓ, ℓ ℓ+1 BTℓ+1 H CC = Ba ∪ Ba ,

(20)

and also collect all the active elements to obtain ℓ ℓ+1 ETℓ+1 H CC = Ea ∪ Ea .

(21)

The above three construction steps are recursively performed until ℓmax − 1. In each recursive step, the active basis ℓ+1 functions BTℓ+1 H CC and the active elements ET H CC are updated as discussed above. Finally we obtain all the active basis functions and active elements from different hierarchical levels, which form the THCCS basis functions and the corresponding THCCS control meshes. During the THCCS construction, we build the control mesh locally at each level. Local selection is employed to select active basis functions and active elements at each level. It also ensures

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

9

the linear independence property of the obtained THCCS basis functions. The proof of linear independence can be obtained by a direct generalization from THB-spline basis functions; see Proposition 2 in the Appendix. In addition, the property of partition of unity is verified in Proposition 3 in the Appendix. The construction of truncated basis functions in Step 1 requires the children basis functions obtained in Step 2. A truncated basis function is actually the same as the original basis function over the active elements at Level ℓ, where the discarded children basis functions have no support. On the other hand, a Level-(ℓ + 1) element within the support of truncated basis functions is defined as a truncated element. A truncated element can be either regular or irregular. Only over truncated elements exist basis functions from different levels. Eq. (18) indicates a tree structure to evaluate a truncated basis function. For each truncated basis function, we first store its children together with the corresponding coefficients ci j , and then set certain ci j to be zero according to Eq. (18). Remark 4.1. Hierarchical B-splines are developed for arbitrary degrees and the degree can vary for different hierarchical levels. However, our development is restricted to cubic polynomials at all levels due to the limitation of a cubic Catmull–Clark basis. In this respect, THCCS is not as flexible as hierarchical B-splines. This is a trade-off for efficiently handling extraordinary nodes for complex geometries. A development based on other degrees requires other types of subdivision basis functions. Remark 4.2. Specific templates were also developed to handle extraordinary nodes in T-splines using the idea of polynomial capping [11,35], where a conforming one-ring neighborhood is required around an extraordinary node. Within the one-ring neighborhood, no T-junctions are allowed in order to obtain a gap-free T-spline geometry. Local refinement of THCCS has no such requirement around extraordinary nodes, and THCCS achieves C 1 -continuity there. 4.2. Example of truncated Catmull–Clark basis functions The truncated Catmull–Clark basis functions are of primal importance to the construction of THCCS. In the following, we take a valence-3 extraordinary node (N = 3), Node 1 in Fig. 5(a), as an example and derive the truncated basis function associated with it. In this study, we only need to consider the two-ring neighborhood elements around the extraordinary node because its associated Catmull–Clark basis function has support contained within them. In Step 1 of our algorithm, we start by defining a set of to-be-refined basis functions Brℓ . To simplify the exposition, we consider only one basis function to be refined at one time. Particularly, we study three cases as shown in ℓ ℓ ℓ Fig. 5(b)–(d): Brℓ = {B4ℓ } (Case 1), Brℓ = {B2N +6 } (Case 2) and Br = {B5 } (Case 3). Corresponding to each case, the ℓ+1 ℓ fine-level domain Ω is the two-ring neighborhood of Br as marked in pink. All the elements in Ω ℓ+1 are identified ℓ as to-be-refined (Er ). According to Eq. (17), the two-ring neighborhood basis functions around Brℓ are identified as to-be-truncated because they have children shared with Brℓ . Without loss of generality, here we study how to construct the truncated basis function for B1ℓ ∈ Btℓ in all three cases; see Fig. 5(a), which is a special basis function associated with the extraordinary node. The other truncated basis functions can be constructed in a similar manner. Fig. 5(e) refers to the refinement of all the elements in Fig. 5(a), where the labeled 6N + 1 basis functions are the children of B1ℓ . According to refinability, B1ℓ can be represented by a linear combination of its children, B1ℓ =

6N +1 

c1 j B ℓ+1 j ,

(22)

j=1

¯ as shown in Fig. 6. where the coefficients c1 j can be obtained from the subdivision matrix A, For each to-be-refined basis function, the refinement is locally performed for its corresponding to-be-refined elements (Erℓ ), as shown in Fig. 5(f)–(h), where dashed lines indicate refinement outside Ω ℓ+1 . Let us consider one to-be-refined basis function at a time, for example B4ℓ (Case 1) in Fig. 5(b), (f). We check all the children of B1ℓ to find out those with support fully contained in Ω ℓ+1 (the pink region), and mark them in green dots. The remaining children are marked in black. All the children basis functions associated with these green dots should be discarded in constructing trunB1ℓ . We then define an index set I ℓ+1 to include all the green dots; we have Case 1: I ℓ+1 = {1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16}; Case 2: I ℓ+1 = {10, 11, 12, 13, 14}; and Case 3: I ℓ+1 = {1, 4, 5, 6, 12, 13, 14, 15, 16}.

10

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

Fig. 5. Construction of the truncated basis function trunB1ℓ . (a) The two-ring neighborhood around a valence-3 extraordinary node at Level ℓ; ℓ ℓ (b) B4ℓ is to be refined (Case 1); (c) B2N +6 is to be refined (Case 2); (d) B5 is to be refined (Case 3); and (e)–(h) are refinement of (a)–(d), respectively. Green dots are discarded in constructing trunB1ℓ , while black dots are kept. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. Coefficients corresponding to the children of B1ℓ .

The truncated basis function trunB1ℓ is then derived in each case by setting c1 j = 0 in Eq. (22) if j ∈ I ℓ+1 . In the same way, we can obtain the truncated basis function for a valence-5 extraordinary node by changing N = 3 to N = 5. The corresponding I ℓ+1 are Case 1: I ℓ+1 = {1, 2, 3, 4, 5, 6, 12, 13, 14, 15, 16, 17, 18, 19, 20}; Case 2: I ℓ+1 = {14, 15, 16, 17, 18}; and Case 3: I ℓ+1 = {1, 4, 5, 6, 15, 16, 17, 18, 19}. The obtained truncated basis functions for a valence-3 and valence-5 extraordinary node are plotted and compared with non-truncated ones in Fig. 7. From Fig. 5(f)–(h), we can observe that we need to discard fifteen, five and nine basis functions (green dots) for Cases 1, 2 and 3, respectively (when N = 3). As shown in Fig. 7(b)–(d) and (f)–(h),

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

11

Fig. 7. Basis functions associated with extraordinary nodes of valence-3 (a)–(d) and valence-5 (e)–(h). (a), (e) are original non-truncated basis functions; (b)–(d) are three truncation cases for (a); and (f)–(h) are three truncation cases for (e). In (b)–(d), (f)–(h), the red dots represent to-berefined basis functions and non-truncated basis functions are visualized in transparency for comparison.

the more basis functions discarded in building the truncated basis function, the more truncation. We do not need to truncate basis functions beyond a two-ring neighborhood (pink domain) of Brℓ . Similarly, we can derive ci j and I ℓ+1 for the other types of basis functions in Fig. 4. The number of children of B1ℓ and the corresponding coefficients c1 j vary with the valence number of the extraordinary node, in contrast with the truncation of standard B-spline basis functions that have a fixed number of children and coefficients. Truncation of a type-2 or type-3 basis function involves ¯ a fixed number of children but various coefficients obtained from A. 5. Geometry preservation In this section, we theoretically study the geometry preservation during h-refinement. THCCS construction starts with an initial Catmull–Clark control mesh, denoted as M0 . We define a sequence of Catmull–Clark control meshes, M0 , . . . , Mℓmax . Mℓ (0 ≤ ℓ ≤ ℓmax ) is generated by ℓ subdivisions of M0 , and we denote the surface obtained from 0 . Recall that evaluation of each Mℓ gives the same Catmull–Clark surface, Sℓ = S0 (ℓ = 0, . . . , ℓ M0 as SCC max ). CC CC After recursive construction, THCCS contains hierarchical levels up to ℓmax . Denote the surface evaluated with the ℓmax 0 THCCS basis functions as SℓTmax H CC . Geometry preservation means that ST H CC = SCC . Note that each level of THCCS ℓ only occupies a portion of the entire surface SℓTmax H CC , and the Level-ℓ THCCS control mesh is only a subset of M . ℓ Correspondingly, the Level-ℓ THCCS basis functions are also a subset of those associated with M . Therefore if we ℓmax 0 ℓmax , then we can prove SℓTmax H CC is exactly the same as the Catmull–Clark surface SCC (= SCC ) evaluated with M can guarantee geometry preservation for THCCS. Proposition 1. Given an initial valid Catmull–Clark control mesh M0 and corresponding Catmull–Clark surface ℓmax 0 , the THCCS surface obtained at levels up to ℓ SCC max (ST H CC ) is exactly the same as the Catmull–Clark surface ℓmax ℓmax ℓmax 0 . evaluated after ℓmax subdivisions (SCC ). We have ST H CC = SCC = SCC Proof. We prove the proposition using induction. It is trivial to verify the initial step since no local refinement is performed. The THCCS basis functions and control mesh are exactly the same as the Catmull–Clark basis functions and 0 . Now we assume that the proposition holds for Level ℓ, 0 ≤ ℓ ≤ ℓ control mesh M0 , we have S0T H CC = SCC max − 1. This assumption means that the THCCS surface evaluated with levels up to ℓ is the same as the Catmull–Clark surface evaluated after ℓ times of subdivision with the control mesh Mℓ , we have  0 ℓ , (23) SℓT H CC = SCC = Piℓ Biℓ = SCC i∈ I¯ℓ

where I¯ℓ is an index set denoting all the basis functions in Mℓ . Among them, the THCCS basis functions B ℓ at Level ℓ are denoted as an index set I ℓ ⊆ I¯ℓ . The set I¯ℓ \I ℓ denotes the basis functions in Mℓ that are not the Level-ℓ THCCS basis functions. Eq. (23) can be rewritten as   SℓT H CC = Piℓ Biℓ + Piℓ Biℓ . (24) i∈I ℓ

i∈ I¯ℓ \I ℓ

12

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

During the construction of THCCS from Level ℓ to Level ℓ + 1, we perform refinement on the to-be-refined basis functions (Brℓ with the index set R ℓ ) and truncation on the to-be-truncated basis functions (Btℓ with the index set T ℓ ). We have R ℓ ⊂ I ℓ and T ℓ ⊂ I ℓ . Therefore in Eq. (24) only the first term on the right needs to be updated. Then, the THCCS surface is evaluated as     Piℓ Biℓ + Piℓ Biℓ , P jℓ+1 B ℓ+1 + Piℓ trunBiℓ + Sℓ+1 j T H CC = (25) ℓ ℓ ℓ ℓ ℓ+1 i∈I \(T ∪R )

i∈T

j∈I

i∈ I¯ℓ \I ℓ

where I ℓ+1 is the index set of all the THCCS basis functions at Level ℓ + 1 (denoted as Baℓ+1 because they are all active). We know that (I ℓ \(T ℓ ∪ R ℓ )) ∪ ( I¯ℓ \I ℓ ) = I¯ℓ \(T ℓ ∪ R ℓ ), so the last two terms in Eq. (25) can be merged together to obtain    Sℓ+1 P jℓ+1 B ℓ+1 + Piℓ trunBiℓ + Piℓ Biℓ . (26) T H CC = j i∈T ℓ

j∈I ℓ+1

i∈ I¯ℓ \(T ℓ ∪R ℓ )

In the following, we check each term in Eq. (26). In the first term on the right, P jℓ+1 at Level ℓ + 1 can be derived ¯ in Eq. (19) and we have from Piℓ via the subdivision matrix A  ci j Piℓ for j ∈ I ℓ+1 . (27) P jℓ+1 = i∈ I¯ℓ

Note that in Eq. (27), ci j = 0 if i ∈ I¯ℓ \(T ℓ ∪ R ℓ ). We now check the third term of Eq. (26) on the right. According to Eq. (13), we can represent Biℓ using their children,  Biℓ = ci j B ℓ+1 for i ∈ I¯ℓ \(T ℓ ∪ R ℓ ). (28) j j∈ I¯ℓ+1 \I ℓ+1

This implies that for the basis function Biℓ with i ∈ I¯ℓ \(T ℓ ∪ R ℓ ), its children are not active THCCS basis functions at Level ℓ + 1. The other Biℓ with i ∈ T ℓ ∪ R ℓ can be either to-be-refined or truncated basis functions. The truncated basis functions in the second term of Eq. (26) can be expressed as  trunBiℓ = ci j B ℓ+1 for i ∈ T ℓ . (29) j j∈ I¯ℓ+1 \I ℓ+1

Note that Eq. (29) also holds when i ∈ R ℓ because all the children B ℓ+1 = 0 when j ∈ I¯ℓ+1 \I ℓ+1 . Therefore in j Eq. (29), we replace i ∈ T ℓ with i ∈ T ℓ ∪ R ℓ . Then, we plug Eqs. (27)–(29) into Eq. (26) and obtain       Sℓ+1 ci j Piℓ B ℓ+1 + Piℓ ci j B ℓ+1 + Piℓ ci j B ℓ+1 T H CC = j j j i∈T ℓ ∪R ℓ j∈ I¯ℓ+1 \I ℓ+1

j∈I ℓ+1 i∈ I¯ℓ

=

  j∈I ℓ+1

=

ci j Piℓ B ℓ+1 + j

i∈ I¯ℓ

 



i∈ I¯ℓ

j∈ I¯ℓ+1 \I ℓ+1

ci j Piℓ B ℓ+1 j

i∈ I¯ℓ \(T ℓ ∪R ℓ ) j∈ I¯ℓ+1 \I ℓ+1

Piℓ ci j B ℓ+1 j



(30)

j∈ I¯ℓ+1 i∈ I¯ℓ

=



ℓ+1 0 = SCC = SCC . P jℓ+1 B ℓ+1 j

j∈ I¯ℓ+1

Eq. (30) indicates that the proposition also holds for Level ℓ + 1. Therefore, it holds for all the levels up to ℓmax , and ℓmax 0 we have SℓTmax  H CC = SCC = SCC . Based on this proposition, we can conclude that THCCS preserves geometry and a THCCS surface inherits the continuity of the Catmull–Clark subdivision surface, namely C 2 -continuity everywhere except C 1 around extraordinary nodes. We adopt C continuity analogous to G continuity, emphasizing the continuity of the parametric space.

13

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20 Table 1 Statistics of all the tested models. Models

Nodes

Elements

Extraordinary nodes

Invalid elements

Refinement steps

Total levels

L-domain1

153 221 3068 3993 3023 1559 2909

128 192 3072 3991 3021 1552 2907

0 4 4 40 36 695 1219

0 0 0 4 4 899 1572

20 40 20 20 20 20 20

5 3 2 2 3 3 3

L-domain2 Genus-3 Hand Bunny Venus Head

Note: See Section 6.2 for elaboration of invalid elements. There are extraordinary nodes in L-domain2 , but none in L-domain1 .

a

ΓD

y

b

c

e

f

r θ

2 Δu =

x

∂ 2u ∂ 2u =0 + ∂x2 ∂x2

2

d

Fig. 8. Laplace equation on the L-shaped domain solved with THCCS. (a) Geometry and problem settings of the L-shaped domain; (b), (c) input control meshes of the case without and with extraordinary nodes, respectively; (d, e) L 2 -error distributions; and (f) convergence curves in both uniform and adaptive refinement.

6. Examples and discussion In this section, we perform adaptive isogeometric analysis using THCCS basis functions. We first start with a widely used benchmark problem for adaptive analysis: solving the Laplace equation over an L-shaped domain. We study cases with and without extraordinary nodes for this problem. Subsequently, we also solve the Laplace equation on complex geometries. Table 1 summarizes the statistics of each model. 6.1. L-shaped domain As shown in Fig. 8(a), a widely used benchmark problem is tested in our adaptive isogeometric analysis: solving the Laplace equation ∆u = 0 over the L-shaped domain [−1, 1]2 \[0, 1]2 with Dirichlet boundary conditions (Γ D ). The analytical solution in polar coordinates (r, θ ) is u(r, θ ) = r 2/3 sin(2θ/3 − π/3),

where r > 0 and π/2 ≤ θ ≤ 2π.

(31)

We duplicate each boundary edge and connect the corresponding nodes with zero-length edges in the parametric space. Then, the obtained elements along the boundary have zero parametric area. They are not evaluated and always set to

14

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

be passive in the analysis. The basis functions over boundary elements are derived with the aid of local knot vectors, which can be obtained in the T-spline manner [3]. Generally speaking, boundary elements have non-uniform local knot vectors. The truncated basis functions can also be applied over boundary elements with the coefficients obtained from mid-knot insertion. Extraordinary nodes are currently not allowed on the boundary in THCCS. They actually need to be separated at least three elements away from the boundary. The reason is that the derivation of explicit Catmull–Clark basis functions is based on uniform cubic B-spline basis functions where knot intervals are required to be the same. Two initial control meshes are utilized: one without extraordinary nodes and the other with four extraordinary nodes; see Fig. 8(b), (c). For each initial control mesh, we conduct both uniform and adaptive refinement using our THCCS scheme. We then obtain a set of active basis functions and active elements. Only active elements need to be evaluated in the analysis. An active element can be either truncated or non-truncated. The evaluation of a nontruncated element, whether regular or irregular, has no difference from the evaluation of an element in Catmull–Clark subdivision [27]. It only requires the information of one hierarchical level and no communication between different hierarchical levels is involved. On the other hand, the evaluation of a truncated element is much more complicated since basis functions from different hierarchical levels may have support over it. The total number of active basis functions varies and truncated basis functions are involved. In our algorithm, we develop a fast bottom-up searching method to collect all the active basis functions for each truncated element. The L 2 -norm error is evaluated using the analytical solution for each element. Usually elements with error larger than a given threshold are identified to be refined. Instead of using the element-wise error Err(Ω ℓj ), in our adaptive approach we convert it to a basis-wise error Err(Biℓ ),  Err(Biℓ ) = Err(Ω ℓj ). (32) Ω ℓj ⊆supp Biℓ

Adaptive refinement is performed and at each refinement step the basis function with the largest error is identified to be refined. Its two-ring neighboring elements are then identified and refined. We perform 20 refinement steps for the case without extraordinary nodes and 40 steps for the case with extraordinary nodes. Fig. 8(d), (e) show the element error distributions for these two cases, respectively. We also plot the convergence curve with respect to the number of degrees of freedom in Fig. 8(f), where we can observe that the case with extraordinary nodes possesses a much larger error in both uniform and adaptive refinement. This is to be expected as it was observed previously [36] that the presence of extraordinary nodes introduces quadrature problem and thus produces larger error, which we will discuss more in Section 6.3. As for the uniform refinement of both cases, we obtain two somewhat different convergence rates: 0.63 and 0.75. Compared to the corresponding uniform refinement, adaptive refinement in both cases achieves the same accuracy with many fewer degrees of freedom. 6.2. Adaptive refinement on complex geometries Here we solve the Laplace equation over complex geometries using THCCS basis functions. Irregular quadrilateral meshes are taken as input control meshes. Recall that a valid control mesh requires that each element contains at most one extraordinary node. However, an irregular quadrilateral mesh has a considerable number of elements that do not satisfy this topology requirement and such elements are invalid. Therefore we need a preprocessing procedure to generate all-valid elements for analysis. In THCCS, we identify a set of to-be-refined basis functions with support covering all the invalid elements, and then we designate all their two-ring neighborhood elements (including valid or invalid ones) as to-be-refined elements. After one step of refinement, extraordinary nodes in invalid elements are separated and all the fine-level elements are valid. In the analysis, Dirichlet boundary conditions are strongly imposed. The L 2 -error is obtained with the aid of socalled bubble functions [6]. Five models are studied: genus-3 model (Fig. 9), hand model (Fig. 10), bunny model (Fig. 11), Venus model (Fig. 12) and head model (Fig. 13). We perform 20 refinement steps for each model. The areas where the Dirichlet boundary conditions are applied are indicated by the arrows in Figs. 9(a)–13(a). Figs. 9(b)–13(b) display the simulation results on the initial valid control meshes whereas Figs. 9(c)–13(c) show the results after 20 refinement steps. Since the genus-3 model has no invalid element, the initial control mesh can be directly used in analysis. However, invalid elements exist in the other models and the input meshes need to be preprocessed. In particular, the Venus and head models contain a large number of such elements, which result in an almost global refinement; see

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

15

Fig. 9. Solving Laplace equation over a genus-3 model. (a) Input quadrilateral mesh with boundary conditions; and (b), (c) simulation results on the initial valid control mesh and the control mesh after 20 refinement steps.

Fig. 10. Solving Laplace equation over a hand model. (a) Input quadrilateral mesh with boundary conditions; and (b), (c) simulation results on the initial valid control mesh and the control mesh after 20 refinement steps.

Fig. 11. Solving Laplace equation over a bunny model. (a) Input quadrilateral mesh with boundary conditions; and (b), (c) simulation results on the initial valid control mesh and the control mesh after 20 refinement steps.

16

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

Fig. 12. Solving Laplace equation over a Venus model. (a) Input quadrilateral mesh with boundary conditions; and (b), (c) simulation results on the initial valid control mesh and the control mesh after 20 refinement steps.

Fig. 13. Solving Laplace equation over a head model. (a) Input quadrilateral mesh with boundary conditions; and (b), (c) simulation results on the initial valid control mesh and the control mesh after 20 refinement steps.

Figs. 12(b) and 13(b). From Figs. 9–13, we can observe that the proposed THCCS method works robustly for complex geometries, which provides a potential wide application in isogeometric analysis with extraordinary nodes involved. 6.3. Limitations Mathematically, the Catmull–Clark basis function and its first order derivatives can be calculated at the extraordinary node. Moreover, the curvature of Catmull–Clark subdivision is square integrable around the extraordinary node [37]. The approximation properties of subdivision surfaces were also studied using interpolation estimates [38], showing that the convergence of approximating a smooth function with the subdivision functions is suboptimal in the H 1 norm. In the above analysis, we always use 4 × 4 Gauss points for the integration over each element, no matter it is regular or irregular. Such numerical integration fails the patch test when extraordinary nodes are involved. This problem was also addressed in Scott’s Ph.D. dissertation [36], and recently it was studied in solving Poisson’s equation on the disk [39]. One solution is to keep subdividing the irregular element into smaller ones until the irregular element is sufficiently small. It was observed that the patch test result can reach the machine precision after about 25 subdivision

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

17

levels [36], which requires excessive Gauss points for the integration. To resolve this problem, a better integration scheme is needed. Although the patch test fails when 4-point Gauss integration rule is applied, in the analysis over the L-shaped domain we observe a similar convergence rate for cases with and without extraordinary nodes. Further study is also needed to better understand the influence of quadrature rule on the convergence behavior when extraordinary nodes are present. 7. Conclusion and future work In this paper we propose a new method termed THCCS to handle extraordinary nodes for complex geometries with support of local refinement. THCCS basis functions satisfy partition of unity and they are linear independent, which are essential properties for geometric design and isogeometric analysis. In addition, THCCS preserves the geometry and inherits the surface continuity property of Catmull–Clark subdivision, namely C 2 -continuity everywhere except C 1 -continuity at extraordinary nodes. Complex surfaces with extraordinary nodes are constructed with the proposed method and adaptive isogeometric analysis is also performed with the THCCS basis functions. THCCS supports a flexible local refinement and there is no propagation of the refined region. From an isogeometric analysis perspective, however, obtaining optimal convergence rates when extraordinary nodes are present remains an open problem. In the future, we will apply our THCCS algorithm to CAD models with sharp features and also extend to trivariate parameterization. Acknowledgments We would like to thank T. Liao for the help in generating the bunny and the hand models, and J.W. Chen for the help in proofreading. X. Wei and Y. Zhang were supported in part by the PECASE Award N00014-14-1-0234 and NSF CAREER Award OCI-1149591. T.J.R. Hughes was supported by grants from the Office of Naval Research (N00014-08-1-0992), the National Science Foundation (CMMI-01101007) and SINTEF (UTA10-000374), with the University of Texas at Austin. Appendix A Throughout the paper, we use basis functions other than blending functions in THCCS, because linear independence is ensured by the local selection mechanism during the THCCS construction. It is worth mentioning that a parallel work on truncated hierarchical subdivision splines was recently developed by Zore et al. [40]. For linear independence, here we refer to global linear independence over the entire parametric domain Ω 0 . As proved by Peters [41] in Lemma 4.1, the Catmull–Clark subdivision functions are globally linearly independent over Ω 0 except for a special case that all the control points are valence-3. Although local linear independence is not valid for Catmull–Clark subdivision functions over any subdomain of Ω 0 , linear independence over an element (regular or irregular) is satisfied for the Catmull–Clark functions that have support on it, unless all the nodes of the element are valence-3. This was proved in Theorem 5.3 by Peters [41]. In practice, an element with all the nodes valence-3 is rare in our control meshes, and we make an assumption to exclude such cases in the construction of THCCS. The following proof of linear independence (Proposition 2) is a direct generalization of THB-spline basis functions [20]. Proposition 2. Given an initial valid Catmull–Clark control mesh M0 in which none of its elements has all nodes with valence-3, the constructed THCCS surface with levels up to ℓmax possesses globally linearly independent basis functions. Proof. We still use induction to prove this proposition. The given valid Catmull–Clark control mesh does not have any element with all its nodes having valence-3. Therefore according to Peters [41], the basis functions associated with M0 are linearly independent on Ω 0 . Then we assume that Proposition 2 holds for levels up to ℓ (0 ≤ ℓ ≤ ℓmax − 1), which means it also holds for Level ℓ. Basis functions at Level ℓ form the domain Ω ℓ and we have  αiℓ Biℓ = 0 if and only if αiℓ |i∈I ℓ = 0, (33) i∈I ℓ

where I ℓ denotes the index set of all the Level-ℓ basis functions. Now we construct Level ℓ + 1 from Level ℓ. Only the Level-ℓ basis functions are influenced, so the remaining basis functions are linearly independent according to the

18

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

assumption. In the following, we only need to prove that basis functions over Ω ℓ are linearly independent, namely    αiℓ trunBiℓ + αiℓ Biℓ + αiℓ+1 Biℓ+1 = 0 (34) i∈T ℓ

i∈I ℓ \(T ℓ ∪R ℓ )

i∈I ℓ+1

if and only if αiℓ |i∈I ℓ \R ℓ = 0

and

αiℓ+1 |i∈I ℓ+1 = 0,

(35)

where αiℓ ∈ R and αiℓ+1 ∈ R are coefficients, T ℓ and R ℓ (T ℓ ⊂ I ℓ and R ℓ ⊆ I ℓ ) are the index sets of truncated and refined basis functions at Level ℓ respectively, I ℓ+1 is the index set of active basis functions at Level ℓ + 1 (they form the domain Ω ℓ+1 ). Let us check the index sets of the first two terms in Eq. (34) and we have T ℓ ∪ (I ℓ \(T ℓ ∪ R ℓ )) = I ℓ \R ℓ ⊆ I ℓ . According to Eq. (18), trunBiℓ is different from Biℓ only in the refined domain Ω ℓ+1 and trunBiℓ stays the same in Ω ℓ \Ω ℓ+1 , that is, trunBiℓ |Ω ℓ \Ω ℓ+1 = Biℓ . Note that Ω ℓ \Ω ℓ+1 is the union of active Level-ℓ elements, and none of such elements has all four nodes with valence-3. According to Theorem 5.3 in Peters [41], the Level-ℓ basis functions that have support on an active Level-ℓ element are linearly independent over this element. Therefore we can conclude that the corresponding coefficients are zero. Looping through all the active Level-ℓ elements, we cover the entire domain Ω ℓ \Ω ℓ+1 , over which we obtain an index set of linearly independent basis functions, I ℓ \R ℓ . Therefore, we can obtain αiℓ |i∈I ℓ \R ℓ = 0. For the third term in Eq. (34), Biℓ+1 are standard Catmull–Clark basis functions, so it is obvious that they are linearly independent. Therefore we have αiℓ+1 |i∈I ℓ+1 = 0. To this end, we prove the THCCS construction from Level ℓ to Level ℓ + 1 results in linearly independent basis functions and Proposition 2 holds for all the levels up to ℓmax . Our THCCS basis functions also satisfy partition of unity. This property can be proved by (1) verifying partition of unity for Catmull–Clark basis functions; and (2) verifying partition of unity for THCCS basis functions. The latter step is exactly the same as the proof of THB-spline basis functions [20] as long as the first step is verified to be true. In the following we focus on the first step. Proposition 3. Explicit Catmull–Clark basis functions satisfy partition of unity. Proof. It is trivial to verify Proposition 3 on regular elements, where Catmull–Clark basis functions are B-spline basis functions and they satisfy partition of unity. On irregular elements, basis functions are expressed in Eq. (9), which can be reorganized as ¯ T (PT b). B = [(V−1 )T 3n−1 VT ]A k

(36)

Note that without loss of generality we set level ℓ = 0. Recall that (3, V) is the eigenstructure of the subdivision ¯ T (PT b). PT b is a vector with 2N +17 elements, matrix A. Therefore (V−1 )T 3n−1 VT = (An−1 )T , and B = (An−1 )T A k k T T Pk b = [α1 , . . . , α2N +17 ] . Among them 16 elements are non-zero and they are values of uniform bi-cubic B-spline 2N +17 basis functions. Therefore, elements of PkT b satisfy partition of unity, and we have i=1 αi = 1. ¯ Denote the ith row and jth column element of A(2N +17)×(2N +8) as a¯ i j . The summation of the elements in the      ¯ T (PT b) is 2N +8 2N +17 a¯ ji αi = 2N +17 αi ( 2N +8 a¯ ji ). Note that 2N +8 a¯ ji is the row summation vector A j=1 i=1 i=1 j=1 j=1 k ¯ A ¯ is the Catmull–Clark subdivision matrix and in each row of A, ¯ only non-zero elements are weight masks. of A.  ¯ equals 1 and we have 2N +17 The summation of each type of masks equals 1. Therefore, the row summation of A i=1 2N +8 2N +17 ¯ T (PT b). This finding indiαi ( j=1 a¯ ji ) = i=1 αi = 1, that is, partition of unity is satisfied for the vector A k cates that if a vector v satisfies partition of unity and the row summation of a matrix M equals 1, the vector MT v ¯ T (PT b), if we can prove the row summation of An−1 is also satisfies partition of unity. Therefore in B = (An−1 )T A k ¯ the row summation 1, then partition of unity for Eq. (36) is satisfied. Since A consists of the first 2N + 8 rows of A, of A equals 1. Let ai j be the ith row and jth column element of A(2N +8)×(2N +8) . The ith row summation of A2 is 2N +8 2N +8 2N +8 2N +8 2N +8 j=1 k=1 aik ak j = k=1 aik ( j=1 ak j ) = k=1 aik = 1. Following the same manner, we can prove that the row summation of An−1 is always 1 no matter what n is. Partition of unity at the extraordinary node (0, 0) is a limit case. Given a parametric value (ξ, η) arbitrarily close to the extraordinary node, we subdivide the irregular element n times until (ξ, η) is located within a regular element. n can be arbitrarily large and the row summation of

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

19

An−1 still remains 1. That is, in the limit case, partition of unity is also satisfied. To this end, we prove partition of unity is satisfied for explicit Catmull–Clark basis functions. References [1] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [2] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogemetric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons, 2009. [3] T.W. Sederberg, J. Zheng, A. Bakenov, A. Nasri, T-splines and T-NURCCs, ACM Trans. Graphics 22 (2003) 477–484. [4] T.W. Sederberg, D.L. Cardon, G.T. Finnigan, N.S. North, J. Zheng, T. Lyche, T-spline simplification and local refinement, ACM Transactions on Graphics 23 (2004) 276–283. [5] J. Deng, F. Chen, X. Li, C. Hu, W. Tong, Z. Yang, Y. Feng, Polynomial splines over hierarchical T-meshes, Graphical Models 70 (2008) 76–86. [6] A.-V. Vuong, C. Giannelli, B. J¨uttler, B. Simeon, A hierarchical approach to adaptive local refinement in isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 200 (2011) 3554–3567. [7] P.B. Bornermann, F. Cirak, A subdivision-based implementation of the hierarchical B-spline finite element method, Comput. Methods Appl. Mech. Engrg. 253 (2013) 584–598. [8] T. Dokken, T. Lyche, K.F. Pettersen, Polynomial splines over locally refined box-partitions, Comput. Aided Geom. Design 30 (2013) 331–356. [9] K.A. Johannessen, T. Kvamsdal, T. Dokken, Isogeometric analysis using LR B-splines, Comput. Methods Appl. Mech. Engrg. 269 (0) (2014) 471–514. [10] Y. Zhang, Y. Bazilevs, S. Goswami, C. Bajaj, T.J.R. Hughes, Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2943–2959. [11] W. Wang, Y. Zhang, G. Xu, T.J.R. Hughes, Converting an unstructured quadrilateral/hexahedral mesh to a rational T-spline, Comput. Mech. 50 (2012) 65–84. [12] W. Wang, Y. Zhang, L. Liu, T.J.R. Hughes, Trivariate solid T-spline construction from boundary triangulations with arbitrary genus topology, in: A Special Issue of Solid and Physical Modeling 2012, in: Comput. Aided Design, 45, 2013, pp. 351–360. [13] Y. Bazilevs, V.M. Calo, J.A. Cottrell, J. Evans, T.J.R. Hughes, S. Lipton, M.A. Scott, T.W. Sederberg, Isogeometric analysis using T-splines, Comput. Methods Appl. Mech. Engrg. 199 (2010) 229–263. [14] M.A. Scott, M.J. Borden, C.V. Verhoosel, T.W. Sederberg, T.J.R. Hughes, Isogeometric finite element data structures based on B´ezier extraction of T-splines, Internat. J. Numer. Methods Engrg. 88 (2011) 126–156. [15] A. Buffa, D. Cho, G. Sangalli, Linear independence of the T-spline blending functions associated with some particular T-meshes, Comput. Methods Appl. Mech. Engrg. 199 (2010) 1437–1445. [16] X. Li, J. Zheng, T.W. Sederberg, T.J.R. Hughes, M.A. Scott, On the linear independence of T-spline blending functions, Comput. Aided Geom. Design 29 (2012) 63–76. [17] M.A. Scott, X. Li, T.W. Sederberg, T.J.R. Hughes, Local refinement of analysis-suitable T-splines, Comput. Methods Appl. Mech. Engrg. 213–216 (2012) 206–222. [18] D.R. Forsey, R.H. Bartels, Hierarchical B-spline refinement, Comput. Graphics 22 (1988) 205–212. [19] R. Kraft, Adaptive and linearly independent multilevel B-splines, in: A.L. M´ehaut´e, C. Rabut, L.L. Schumaker (Eds.), Surface Fitting and Multiresolution Methods, Vanderbilt University Press, 1997, pp. 209–218. [20] C. Giannelli, B. J¨uttler, H. Speleers, THB-splines: The truncated basis for hierarchical splines, Comp. Aided Geom. Design 29 (2012) 485–498. [21] T.J. Cashman, U.H. Augsd¨orfer, N.A. Dodgson, M.A. Sabin, NURBS with extraordinary points: high-degree, non-uniform, rational subdivision schemes, ACM Transactions on Graphics 28 (46) (2009) 1–46. 9, July. [22] U. Reif, A unified approach to subdivision algorithms near extraordinary vertices, Comp. Aided Geom. Design 12 (1995) 153–174. [23] Q. Pan, G. Xu, Y. Zhang, Unified method for hybrid subdivision surface design using geometric partial differential equations, in: A Special Issue of Solid and Physical Modeling 2013, in: Comput. Aided Design, 46, 2014, pp. 110–119. [24] D. Zorin, P. Schr¨oder, Subdivision for modeling and animation, in: ACM Siggraph Course Notes, 2000. [25] M. Sabin, Recent progress in subdivision: a survey, in: Advances in Multiresolution for Geometric Modelling, 2005, pp. 203–230. [26] E. Catmull, J. Clark, Recursively generated B-spline surfaces on arbitrary topological meshes, Comput. Aided Design 10 (1978) 350–355. [27] J. Stam, Exact evaluation of Catmull-Clark subdivision surfaces at arbitrary parameter values, in: Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques, 1998, pp. 395–404. [28] Hierarchical subdivision surfaces, 2004, http://renderman.pixar.com/resources/current/rps/hierarchsubdiv.html. [29] F. Cirak, M. Ortiz, P. Schr¨oder, Subdivision surfaces: a new paradigm for thin shell analysis, Internat. J. Numer. Methods Engrg. 47 (2000) 2039–2072. [30] F. Cirak, M.J. Scott, E.K. Antonsson, M. Ortiz, P. Schr¨oder, Integrated modeling, finite-element analysis, and engineering design for thin shell structures using subdivision, Comput. Aided Design 34 (2002) 137–148. [31] D. Burkhart, B. Hamann, G. Umlauf, Iso-geometric finite element analysis based on Catmull-Clark subdivision solids, Comput. Graph. Forum 29 (2010) 1575–1584. [32] E. Grinspun, P. Krysl, P. Schr¨oder, CHARMS: A simple framework for adaptive simulation, ACM Transactions on Graphics 21 (2002) 281–290. [33] M.A. Scott, D.C. Thomas, E.J. Evans, Isogeometric spline forests, Comput. Methods Appl. Mech. Engrg. 269 (0) (2014) 222–264. [34] C. Giannelli, B. J¨uttler, H. Speleers, Strongly stable bases for adaptively refined multilevel spline spaces, Adv. Comput. Math. (2013) 1–32.

20

X. Wei et al. / Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20

[35] L. Liu, Y. Zhang, T.J.R. Hughes, M.A. Scott, T.W. Sederberg, Volumetric T-spline construction using Boolean operations, Eng. Comput. (2014) http://dx.doi.org/10.1007/s00366-013-0346-6. [36] M.A. Scott, T-splines as a Design-Through-Analysis technology (Ph.D. thesis), The University of Texas at Austin, 2011. [37] U. Reif, P. Schr¨oder, Curvature integrability of subdivision surfaces, Adv. Comput. Math. 14 (2) (2001) 157–174. [38] G. Arden, Approximation properties of subdivision surfaces (Ph.D. thesis), University of Washington, 2001. [39] T. Nguyen, K. Karˇciauskas, J. Peters, A comparative study of several classical, discrete differential and isogeometric methods for solving Poisson’s equation on the disk, Axioms 3 (2014) 280–300. [40] U. Zore, B. J¨uttler, J. Kosinka, On the linear independence of (truncated) hierarchical subdivision splines, Geometry + Simulation Report (17) (2014). [41] J. Peters, X. Wu, On the local linear independence of generalized subdivision functions, SIAM J. Numer. Anal. 44 (6) (2006) 2389–2407.