C2 interpolation T-splines

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 362 (2020) 112835 www.elsevier.com/locate/cma

C 2 interpolation T-splines Yuanpeng Zhua ,∗, Xuli Hanb a b

School of Mathematics, South China University of Technology, Guangzhou, 510641, PR China School of Mathematics and Statistics, Central South University, Changsha, 410083, PR China Received 7 October 2018; received in revised form 6 January 2020; accepted 6 January 2020 Available online xxxx

Abstract In this work, we construct a class of C 2 interpolation T-spline (called IT-spline) basis functions over regular T-meshes. A new rule for determining the local knot vectors of the IT-spline basis functions is designed. The resulting IT-spline basis functions are always linearly independent for regular T-meshes. In order to guarantee the partition of unity and compact support properties of the IT-spline basis functions, several modified schemes are developed. An iteration adjustment algorithm for adjusting the points on IT-spline surface to the targeted points is proposed. Some applications of the new IT-spline basis functions in fitting bivariate functions and open 3D meshes, solving numerical PDEs and IGA are given. c 2020 Elsevier B.V. All rights reserved. ⃝ MSC: 65D07; 65D17 Keywords: T-mesh; T-splines; Interpolation; Isogeometric analysis

1. Introduction In Computer Aided Geometric Design (CAGD) and Isogeometric Analysis (IGA), locally refinable splines have gained widespread attention. The motivation to apply such splines in design and analysis is to overcome the bottleneck of global tensor product structure of NURBS. By allowing T-junctions in the control meshes, T-splines based on local knot vectors were proposed by Sederberg and his colleagues [1,2], which have shown great potential applications in data fitting, geometric modeling and adaptive IGA [3–5]. However, as pointed out by Buffa et al. [6], the blending functions of T-splines associated with several particular T-meshes may be linearly dependent. This implies that the whole class of T-splines are not suitable for IGA. In order to fix the linear dependence of T-spline blending functions, by constraining the associated T-meshes having no intersections of T-junction extensions, analysis suitable T-splines (AST-splines for short) were proposed in [7–9]. AST-splines are a subset of T-splines and possess some desired properties for design and analysis, such as partition of unity, compact support and linear independence, see [10]. However, the topology of the T-meshes of AST-splines is relatively restricted and local refinement in AST-splines may propagate beyond the domain of interest. Recently, ∗ Corresponding author.

E-mail addresses: [email protected] (Y. Zhu), [email protected] (X. Han). https://doi.org/10.1016/j.cma.2020.112835 c 2020 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

analysis-suitable++ (AS++) T-splines were developed in [11,12], which include AST-splines as a special case and maintain all the good mathematical properties of AST-splines. In [13,14], modified T-splines and truncated T-splines were proposed. To address the problem of local refinement of splines on rectangular domain, the concept of hierarchical B-splines (HB-splines for short) was firstly introduced by Forsey and Bartels as an accumulation of B-splines with nested knot vectors [15]. HB-splines can be locally refined using overlays. In [16], Giannelli et al. normalized HB-splines by using a truncation mechanism to reduce the support of basis functions. They call such HB-splines as truncated hierarchical B-splines (THB-splines for short). The basis functions of THB splines are non-negative, locally supported, linearly independent and form a partition of unity, and allow an effective local control of refinement. LR-splines [17] also possess local refinement property by inserting line segment in a mesh. Recently, Johannessen et al. discussed the similarities and differences between HB-splines, THB-splines and LR-splines in [18]. In 2006, Deng, Chen and Feng introduced the concept of polynomial spline spaces over T-meshes and established a dimension formula for the spline space with some constrains on the degrees of the polynomials and the smoothness [19]. Later, in [20], a set of bi-cubic C 1 polynomial splines defined on hierarchical T-meshes (PHT-splines for short) was constructed by using truncation mechanism. PHT-splines possess many nice properties, such as linear independence, non-negativity, partition of unity, compact support, and existing an efficient and simple local refinement algorithm. These properties make PHT-splines attractive for geometric processing, adaptive finite elements and IGA [21–25]. Despite the excellent results reported, the limited smoothness of PHT-splines is a major shortcoming with regard to geometric modeling. NURBS, T-splines, AST-splines, AS++ T-splines, HB-splines, THB-splines, LR-splines and PHT-splines mentioned above are approximation splines. For some applications or some users, approximation splines are preferable, whereas for others, interpolation splines are imperative. In geometric design, approximation splines possessing partition of unity and non-negativity properties are used usually, by which the resulting parametric curves and surfaces have convex hull property. In finite element analysis, basis functions are typically interpolatory and may take on positive and negative values, see [26]. In IGA, the non-interpolatory property of NURBS necessitates special treatment the essential boundary conditions, see [27,28]. Recently, Escobar and his colleagues proposed a new method to construct a trivariate T-spline representation of complex genus-zero solids for the application of IGA [29]. In the method, the control points of the trivariate T-spline are calculated by imposing the interpolation conditions on points sited both on the inner and on the surface of the solid, which involves solving linear systems. Since the blending functions of T-splines may be linearly dependent for some special T-meshes, the method has some potential limits for practical applications. For scattered data interpolation, a sufficient condition for multilevel bi-cubic B-spline approximation to generate an interpolation function through scattered data points was developed in [30]. In the method, interpolation is achieved by solving linear systems when the control point spacing in the finest lattice becomes sufficiently small relative to the data distribution. In [31], Zhang and Ma made use of sinc function and a Gaussian multiplier to construct a set of interpolation basis functions. The resulting interpolatory curves and surfaces are C ∞ continuous. Nevertheless, the interpolation basis functions have no exact partition of unity as well as no compact support. Later, Yuan and Ma further used sinc function, a Gaussian multiplier and a Cosine multiplier to construct a kind of C 2 truncated interpolation basis functions and they showed that the new bases have some merits in representing special geometry and IGA [32]. The C 2 truncated interpolation basis functions have compact support but they still have no exact partition of unity. The novel contributions of this work are the constructions of C 2 interpolation T-spline (IT-spline for short) basis functions over regular T-meshes with exact partition of unity and linear independence. The rest of this paper is organized as follows. Section 2 gives the construction and properties of a kind of univariate C 2 quartic polynomial interpolation basis functions. In Section 3, a kind of bi-quartic C 2 IT-spline basis functions over regular T-meshes is constructed in detail. Several new modified schemes for guaranteeing the partition of unity and compact support properties of the C 2 IT-spline basis functions are developed. And a new iteration adjustment algorithm for adjusting the points on IT-spline surface to the targeted points is given in detail. In Section 4, applications of the new IT-spline basis functions in fitting bivariate functions and open 3D meshes, solving numerical PDEs and IGA are shown. Finally, we provide conclusions in the Section 5.

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

3

2. C 2 quartic polynomial interpolation splines 2.1. Quartic polynomial interpolation basis Given a knot vector U = [u −1 , u 0 , . . . , u n+1 , u n+2 ], for u ∈ [u i , u i+1 ], let h i = u i+1 −u i , then for i = 1, 2, . . . , n, a class of quartic polynomial interpolation basis functions is given as follows Bi (u) : = B [u i−2 , u i−1 , u i , u i+1 , u i+2 ] (u) ⎧ u−u Bi,0 ( h i−2 ), u ∈ [u i−2 , u i−1 ), ⎪ ⎪ i−2 ⎪ ⎪ ⎪ u−u i−1 ⎪ ⎪ ⎪ Bi,1 ( h i−1 ), u ∈ [u i−1 , u i ), ⎨ i = Bi,2 ( u−u ), u ∈ [u i , u i+1 ), hi ⎪ ⎪ ⎪ u−u ⎪ ⎪ Bi,3 ( h i+1 ), u ∈ [u i+1 , u i+2 ), ⎪ i+1 ⎪ ⎪ ⎩ 0, u∈ / [u i−2 , u i+2 ),

(1)

where 2 −h i−2 (1 − x) x 3 , (h i−2 + h i−1 )h i−1 [ 1 Bi,1 (x) = h i−2 h i (1 − x)3 x + 3 (h i−2 + h i−1 ) h i (1 − x)2 x 2 (h i−2 + h i−1 )h i ] + (h i−2 + h i−1 ) (h i−1 + 3h i ) (1 − x) x 3 + (h i−2 + h i−1 ) h i x 4 , [ 1 h i−1 (h i + h i+1 ) (1 − x)4 + (3h i−1 + h i ) (h i + h i+1 ) (1 − x)3 x Bi,2 (x) = h i−1 (h i + h i+1 ) ] + 3h i−1 (h i + h i+1 ) (1 − x)2 x 2 + h i−1 h i+1 (1 − x) x 3 ,

Bi,0 (x) =

Bi,3 (x) =

2 −h i+1 (1 − x)3 x . h i (h i + h i+1 )

Remark 1. In [33], a class of general quartic polynomial splines with three local shape parameters was constructed. For special choices of the three local shape parameters, the general quartic polynomial splines given in [33] include the quartic polynomial interpolation basis functions given in (1) as a special case. For multiple knots, it is worth noticing that we shrink the corresponding intervals to zero and drop the corresponding pieces of the interpolation basis function. For example, if u i−2 = u i−1 is a double knot, that is h i−2 = 0, we define ⎧ u−u i−1 ⎪ ⎪ Bi,1 ( h i−1 ), u ∈ [u i−1 , u i ), ⎪ ⎪ ⎨ B ( u−u i ), u ∈ [u i , u i+1 ), i,2 h i B [u i−1 , u i−1 , u i , u i+1 , u i+2 ] (u) = u−u i+1 ⎪ Bi,3 ( h ), u ∈ [u i+1 , u i+2 ), ⎪ ⎪ i+1 ⎪ ⎩ 0, u∈ / [u i−2 , u i+2 ). It is easy to check that B [u i−1 , u i−1 , u i , u i+1 , u i+2 ] (u) has C 1 continuity at u = u i−1 . For 1 ≤ i ≤ n and distinct knots, by directly computing, we have ( ± ) ( ± ) = 0, Bi u i−2 = 0, Bi′ u i−2 ( ± ) ( ) h i−2 ± Bi u i−1 = 0, Bi′ u i−1 = , (h i−2 + h i−1 ) h i−1 ( ) ( ) h i − h i−1 Bi u i± = 1, Bi′ u i± = , h i−1 h i ( ± ) ( ± ) h i+1 Bi u i+1 = 0, Bi′ u i+1 =− , h i (h i + h i+1 )

( ± ) Bi′′ u i−2 = 0, ( ) 6 ± Bi′′ u i−1 , = (h i−2 + h i−1 ) h i−1 ( ) 6 Bi′′ u i± = − , h i−1 h i ( ± ) 6 Bi′′ u i+1 = , h i (h i + h i+1 )

4

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

( ± ) Bi u i+2 = 0,

( ± ) = 0, Bi′ u i+2

( ± ) = 0, Bi′′ u i+2

which implies that Bi (u) has C 2 continuity at each of the knots and has the interpolation property Bi (u j ) = δi j . Moreover, for u ∈ [u i , u i+1 ], t = (u − u i )/ h i , i = 2, 3, . . . , n − 2, we have n ∑ Bi (u) = (1 − t)4 + 4(1 − t)3 t + 6(1 − t)2 t 2 + 4 (1 − t) t 3 + t 4 = [(1 − t) + t]4 ≡ 1, (2) i=1

which implies that Bi (u), i = 1, 2, . . . , n form a partition of unity on [u 2 , u n−1 ]. 2.2. Quartic polynomial interpolation curve Given a knot vector U = [u −1 , u 0 , . . . , u n+1 , u n+2 ] and control points Pi (i = 1, . . . , n) in Rd , d ∈ {1, 2, 3}, then, for u ∈ [u 1 , u n ], by using the given C 2 interpolation basis functions Bi (u), i = 1, 2, . . . , n, we define a kind of C 2 quartic polynomial interpolation curves as follows n ∑ P(u) = Pi Bi (u). (3) i=1

From (3), it is obvious that P(u i ) = Pi , for all 1 ≤ i ≤ n. We further give a geometric construction of the interpolation curve P(u) by a geometric description of each curve segment in terms of B´ezier control points. For u ∈ [u i , u i+1 ], t = (u − u i )/h i , 2 ≤ i ≤ n − 2, after some manipulations, we can rewrite the interpolation curve P(u) as the following quartic B´ezier curve segment form ) ) ( ( h i ∆i h i di 4(1 − t)3 t + Pi + 6(1 − t)2 t 2 P(u) = Pi (1 − t)4 + Pi + 4 2 ) ( h i di+1 4 (1 − t) t 3 + Pi +1 t 4 , (4) + Pi+1 − 4 where 1 di = (h i−1 ∆i + h i ∆i−1 ) , h i−1 + h i and ∆i := (Pi+1 − Pi )/h i . We give some shape analysis on the interpolation curve segment P(u)(u ∈ [u i , u i+1 ]) for the case of very small h i . For very small h i , it is obvious that both di and di+1 tend to ∆i . From (4), it follows that for u ∈ [u i , u i+1 ], t = (u − u i )/h i , the interpolation curve segment P(u) tends to a line segment (1 − t)Pi + t Pi+1 . Therefore, when the distance between Pi and Pi+1 is small, it is valid to set small h i correspondingly. From (4), it is easy to check that the given interpolation basis functions reproduce local linear polynomials, that is if the data set is given locally by a linear polynomial function, then the interpolant will locally reproduce the linear polynomial function. In fact, for any constants A, C in R, if P j = Au j + C, j = i − 1, i, i + 1, i + 2, then ∆i−1 = ∆i = ∆i+1 = A, thus for u ∈ [u i , u i+1 ], t = (u − u i )/h i , 2 ≤ i ≤ n − 2, from (4), we have ) ( ) [ ( hi hi 4(1 − t)3 t + u i + 6(1 − t)2 t 2 P(u) = A u i (1 − t)4 + u i + 4 2 ( ) ] 3h i + ui + 4 (1 − t) t 3 + (u i + h i ) t 4 + C 4 = A (u i + h i t) + C = Au + C,

(5)

which implies that the interpolant can locally reproduce linear polynomial. From this result, by using the standard estimates for polynomial approximation, see [34], it can be easily proved that, for any f ∈ C 2 ([u 1 , u n ]) and interpolation data (u j , P j ) with P j = f (u j ), we have f (u) − P (u) = O(h 2 ), where u ∈ [u 2 , u n−1 ], h = max2≤i≤n−2 h i . Fig. 1 shows the C 2 quartic polynomial interpolation basis functions (on the left) and the corresponding parametric interpolation curves (on the right) for the Akima’s data points, which can be found in [35]. The non-uniform knot vector is [−6, −5, 1, 2, 3, 4.8, 5.3, 5.5, 5.7, 7.1, 9, 10, 11, 15, 16].

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

5

Fig. 1. C 2 quartic polynomial interpolation curves.

3. C 2 interpolation T-splines 3.1. T-meshes and T-splines In this section, we will briefly recall some pertinent concepts concerning T-meshes and T-splines. More details can be found in [1–3]. The control grid of a T-spline surface is called a T-mesh, which is basically a rectangular grid allowing Tjunctions. The pre-image of each edge in a T-mesh is a line segment of constant u (which we will call an u-edge) or of constant v (which we will call a v-edge). Fig. 2(a) illustrates the pre-image of a T-mesh in the (u, v) parameter space. In the pre-image, the knot intervals di and ei are positive numbers indicating the differences between two knots, and these are assigned to edges of the T-mesh. A valid T-mesh requires the sum of all knot intervals along one side of any face to equal the sum of the knot intervals on the opposite side. For instance, in Fig. 2(a), e3 = e6 + e7 on face f 1 , and d2 = d6 + d7 on face f 2 . In the pre-image of a T-mesh, if a T-junction on edge of a face can be connected to a T-junction on an opposing edge of the face thereby splitting the face into two rectangles, that edge must be included in the T-mesh. From the knot intervals of a T-mesh, we can designate the origin (0, 0) of the (u, v) parameter space so as to establish a knot coordinate system. For example, in Fig. 2(a), if we choose (u 1 , v1 ) = (0, 0), then u 2 = d1 , u 3 = d1 + d2 , v2 = e1 , v3 = e1 + e2 , and so on. Therefore, the knot coordinate for P1 is (0, 0), for P2 is (u 2 , v2 + e5 ), and for P3 is (u 4 , v3 + e6 ). A T-spline is a point-based spline. Each of its control points Pi with knot coordinate (u i2 , vi2 ) corresponds to a C 2 bi-cubic T-spline blending function as follows Fi (u, v) = N [u i0 , u i1 , u i2 , u i3 , u i4 ] (u) N [vi0 , vi1 , vi2 , vi3 , vi4 ] (v) .

(6)

6

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

Fig. 2. An example of Pre-image of a T-mesh and inferring of local knot vectors for the T-spline blending function F3 (u, v) corresponding to P3 by Rule 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In (6), N [u i0 , u i1 , u i2 , u i3 , u i4 ] (u) is a cubic B-spline basis function defined on the u-directional knot vector, ui = [u i0 , u i1 , u i2 , u i3 , u i4 ] ,

(7)

and N [vi0 , vi1 , vi2 , vi3 , vi4 ] (v) is another cubic B-spline basis function defined on the v-directional knot vector, vi = [vi0 , vi1 , vi2 , vi3 , vi4 ] .

(8)

For example, ⎧ (u−u i0 )3 ⎪ , u ∈ [u i0 , u i1 ) , ⎪ ⎪ −u (u )(u i1 i0 i2 −u i0 )(u i3 −u i0 ) ⎪ ⎪ 2 ⎪ (u−u i0 ) (u i2 −u) ⎪ i0 )(u−u i1 )(u i3 −u) ⎪ + (u (u−u ⎪ (u i2 −u i0 )(u i3 −u i0 )(u i2 −u i1 ) ⎪ i3 −u i0 )(u i2 −u i1 )(u i3 −u i1 ) ⎪ ⎪ ⎪ (u−u i1 )2 (u i4 −u) ⎪ ⎪ ⎪+ (u i2 −u i1 )(u i3 −u i1 )(u i4 −u i1 ) , u ∈ [u i1 , u i2 ) , ⎨ (u−u i0 )(u i3 −u)2 i3 −u)(u i4 −u)(u−u i1 ) N [u i0 , u i1 , u i2 , u i3 , u i4 ] (u) = + (u (u−u (u i3 −u i0 )(u i3 −u i1 )(u i3 −u i2 ) ⎪ i3 i2 )(u i4 −u i1 )(u i3 −u i1 ) ⎪ ⎪ ⎪ (u i4 −u)2 (u−u i2 ) ⎪ + (u −u )(u −u )(u −u ) , u ∈ [u i2 , u i3 ) , ⎪ ⎪ i3 i2 i4 i2 i4 i1 ⎪ ⎪ ⎪ (u i4 −u)3 ⎪ ⎪ , u ∈ [u i3 , u i4 ) , ⎪ (u i4 −u i1 )(u i4 −u i2 )(u i4 −u i3 ) ⎪ ⎪ ⎪ ⎩0, u∈ / [u , u ) . i0

i4

The knot vectors ui and vi given in (7) and (8) are inferred from the T-mesh neighborhood of Pi by the following rule, see [1–3] (Fig. 2). Rule 1. The knot vectors ui and vi for the T-spline blending function Fi (u, v) corresponding to Pi are determined as follows (see Fig. 2). Suppose the knot coordinate of Pi is (u i2 , vi2 ). Consider a ray in the parameter space, R(α) = (u i2 + α, vi2 ). Then, u i3 and u i4 are the u coordinates of the first two u-edges intersected by the ray (not including the initial point (u i2 , vi2 )). The other knots in ui and vi are found in a similar manner. By u-edge, we mean a vertical line segment of constant u. For instance, in Fig. 2(a), by Rule 1, the knot vectors for P3 are u3 = [u 2 + d6 , u 3 , u 4 , u 5 , u 6 ], v3 = [v2 +e5 , v3 , v3 +e6 , v4 , v5 ]. P1 is a boundary control point. In this case, u 10 , u 11 , v10 and v11 do not matter, so we can take u1 = [u 1 −d, u 1 −d, u 1 , u 2 , u 3 ] and v1 = [v1 −e, v1 −e, v1 , v2 , v3 ]. For P2 , we take u2 = [u 1 −d, u 1 , u 2 , u 3 , u 4 ] and v2 = [v1 , v2 , v2 +e5 , v3 , v4 ]. For P4 , u4 = [u 2 +d6 , u 5 , u 5 +d8 , u 6 , u 6 +d] and v4 = [v2 +e5 , v4 , v4 +e8 , v5 , v5 +e]. And for P5 , u5 = [u 4 , u 5 , u 6 , u 6 + d, u 6 + d] and v5 = [v1 − e, v1 , v2 + e5 , v4 , v4 + e8 ]. In general, T-spline blending functions do not span a polynomial space. From Fig. 2(b), it is observed that in some cells of the mesh F3 (u, v) is not a polynomial, but a piecewise polynomial. Some additional restrictions on the T-mesh configuration should be satisfied in order to span a polynomial spline space, see [7]. If these restrictions are not verified, T-spline blending functions should be normalized in order to form a partition of unity. This leads

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

7

Fig. 3. An example of regular T-mesh and non-regular T-mesh.

to rational T-spline blending functions ρi Fi (u, v) Ti (u, v) = ∑n , i=1 ρi Fi (u, v) where ρi , i = 1, 2, . . . , n are positive weights. After getting all the rational T-spline blending functions Ti (u, v) corresponding to the control points Pi , i = 1, 2, . . . , n, a T-spline surface T (u, v) can be generated by blending the control points Pi and the functions Ti (u, v), i = 1, 2, . . . , n, i.e., n ∑ T (u, v) = Pi Ti (u, v) . i=1

Recently, Buffa et al. [6] have analyzed the linear dependence of the bi-cubic T-spline blending functions corresponding to some particular T-meshes, which means that not all T-splines are suitable as a basis for IGA. In order to fix the linear dependence of T-splines, a kind of analysis suitable T-splines was proposed in [8,9], which is a mildly topological restricted subset of T-splines. T-splines are approximation splines over T-meshes, which implies that T (u i2 , vi2 ) ̸= Pi in general. For some applications or some users, approximation splines are preferable, whereas for others, interpolation splines are imperative, see [29]. In the following subsection, we shall construct a class of C 2 interpolation T-splines over regular T-meshes by using the C 2 interpolation basis functions given in (1). 3.2. Construction of C 2 IT-spline basis functions In the following discussion, we add two constrains on the T-mesh. We assume that each vertex in the T-mesh must be on two other grid lines, and in the pre-image of the T-mesh, each cell or facet in the grid of must be a rectangle. These two constrains on the T-mesh have been used in [13,19]. We call such T-mesh as regular T-mesh. Fig. 3 gives an example of regular T-mesh and non-regular T-mesh. There are two reasons why the T-mesh given in Fig. 3(b) is non-regular. The first reason is that it includes three vertexes P2 , P6 and P7 , which are not on any two other grid lines. And the second reason is that it includes a cell or facet marked with gray color, which is a L-shape domain other than a rectangle. For constructing C 2 IT-spline basis functions, similar to that done in constructing T-spline blending functions, for each of control points Pi with knot coordinate (u i2 , vi2 ), i = 1, 2, . . . , n, by using the C 2 interpolation basis functions given (1), we designate a corresponding C 2 IT-spline basis function as follows Bi (u, v) = B [u i0 , u i1 , u i2 , u i3 , u i4 ] (u) B [vi0 , vi1 , vi2 , vi3 , vi4 ] (v) . (9) ( ) In order to ensure the interpolation conditions Bi u j2 , v j2 = δi j , the knot vectors ui = [u i0 , u i1 , u i2 , u i3 , u i4 ] and vi = [vi0 , vi1 , vi2 , vi3 , vi4 ] in (9) are inferred from the T-mesh neighborhood of Pi by the following new Rule 2, whose main idea is to ensure every vertex lying in the support domain of an IT-spline basis function must be on the knot lines. Rule 2 is slightly different from Rule 1.

8

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

Fig. 4. An example of Pre-image of a T-mesh and inferring of local knot vectors for the IT-spline basis function B3 (u, v) corresponding to P3 by Rule 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Rule 2. The knot vectors ui and vi for the IT-spline basis function Bi (u, v) corresponding to Pi are determined as follows (see Fig. 4). Suppose the knot coordinate of Pi is (u i2 , vi2 ). Consider a ray in the parameter space, R(α) = (u i2 + α, vi2 ). Then, u i3 and u i4 are the u coordinates of the first two extended u-edges intersected by the ray (not including the initial point (u i2 , vi2 )). The other knots in ui and vi are found in a similar manner. By extended u-edge we mean that extend an u-edge in two directions till it intersects the first v-edge in each of directions (refer to Fig. 4). Similarly, an extended v-edge means that extend an v-edge in two directions till it intersects the first u-edge in each of directions (refer to Fig. 4). For instance, by Rule 2, in Fig. 4(a), the knot vectors for P3 are u3 = [u 2 + d6 , u 3 , u 4 , u 5 , u 5 + d8 ], v3 = [v2 + e5 , v3 , v3 + e6 , v4 , v4 + e8 ]. P1 is a boundary control point. In this case, u 10 , u 11 , v10 and v11 do not matter, we take u1 = [u 1 −d, u 1 −d, u 1 , u 2 , u 3 ] and v1 = [v1 −e, v1 −e, v1 , v2 , v2 +e5 ]. For P2 , u2 = [u 1 −d, u 1 , u 2 , u 2 +d6 , u 3 ] and v2 = [v1 , v2 , v2 +e5 , v3 , v3 +e6 ], for P4 , u4 = [u 4 , u 5 , u 5 +d8 , u 6 , u 6 +d] and v4 = [v3 +e6 , v4 , v4 +e8 , v5 , v5 +e], and for P5 , u5 = [u 5 , u 5 + d8 , u 6 , u 6 + d, u 6 + d] and v5 = [v1 , v2 , v2 + e5 , v3 + e6 , v4 ]. We discuss some properties satisfied by the IT-spline basis functions Bi (u, v), i = 1, 2, . . . , n. Theorem 1. For a regular T-mesh, with the knot vectors ui = [u i0 , u i1 , u i2 , u i3 , u i4 ], vi = [vi0 , vi1 , vi2 , vi3 , vi4 ] determined by Rule 2, the resulting IT-spline basis functions Bi (u, v), i = 1, 2, . . . , n have the following properties (A) C 2 continuity. (B) Compact support, i.e. suppBi (u, v) = {(u, v) |u ∈ [u i0 , u i4 ] , v ∈ [vi0 , vi4 ] }. ( ) (C) Interpolation property, i.e. Bi u j2 , v j2 = δi j , where (u j2 , v j2 ) is the knot coordinate of P j . (D) Linear independence. (E) For 3 × 3 tensor-product submesh, the IT-spline basis functions Bi (u, v), i = 1, 2, . . . , n possess local partition of unity and local linear polynomials reproduction property. Proof. The properties (A), (B) are obvious. We shall prove (C), (D) and (E). (C) For a regular T-mesh, suppose that by Rule 2, the resulting knot vectors ui and vi of Bi (u, v) are ui = [u i0 , u i1 , u i2 , u i3 , u i4 ] and vi = [vi0 , vi1 , vi2 , vi3 , vi4 ], respectively. Then for any knot coordinate (u j2 , v j2 ) ∈ / suppBi (u, v) or (u j2 , v j2 ) ∈ ∂suppBi (u, v), it is obvious that Bi (u j2 , v j2 ) = 0. Moreover, from the interpolation property of the univariate interpolation basis functions given in (1), for any knot coordinate (u j2 , v j2 ) = (u ik , vil ), 1 ≤ k ≤ 3 and 1 ≤ l ≤ 3, we have Bi (u j2 , v j2 ) = δi j . We can conclude the result by further proving that there are at most nine possible knot coordinates (u ik , vil ), k = 1, 2, 3 and l = 1, 2, 3 within the interior of suppBi (u, v). In fact, if there is a possible knot coordinate (u ∗j2 , v ∗j2 ) differing from any one of (u ik , vil ), k = 1, 2, 3 and l = 1, 2, 3 within the interior of suppBi (u, v), then for a regular T-mesh, by Rule 2, we can see that there must at least be an intersection between the extended u-edge passing (u ∗j2 , v ∗j2 ) and one of the two rays R(β) = (u i2 , vi2 + β) or between the extended v-edge passing (u ∗j2 , v ∗j2 ) and one of the two rays

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

9

R(α) = (u i2 + α, vi2 ), which implies that at least one of the two u ∗j2 and v ∗j2 should be one of the components of the resulting knot vectors ui and vi . While this contradicts the hypothesis since we have supposed that by Rule 2, the resulting knot vectors ui and vi are ui = [u i0 , u i1 , u i2 , u i3 , u i4 ] and vi = [vi0 , vi1 , vi2 , vi3 , vi4 ], respectively. (D) For λi ∈ R, i = 1, 2, . . . , n, we consider a linear combination n ∑ λi Bi (u, v) = 0, i=1

then for (u, v) = (u j2 , v j2 ), we have n ∑

( ) λi Bi u j2 , v j2 = λ j = 0,

i=1

where j = 1, 2, . . . , n. (E) Among all the control points Pi , i = 1, 2, . . . , n, if there are exactly 4 × 4 control points Pkl (k = s −1, s, s +1, s +2, l = m −1, m, m +1, m +2) laying on a rectangular mesh, that is Pkl = Au k + Bvl +C, with the corresponding knot coordinates (u k , vl ) on a 3 × 3 tensor-product submesh and the associated IT-spline basis functions denoted as Bkl (u, v), then from the partition of unity and local linear polynomials reproduction property of the univariate interpolation basis functions, see (2) and (5), for any (u, v) ∈ [u s , u s+1 ]×[vm , vm+1 ], we have n s+2 m+2 ∑ ∑ ∑ Bi (u, v) = Bkl (u, v) ≡ 1, i=1

k=s−1 l=m−1

and n ∑ i=1

Pi Bi (u, v) =

s+2 m+2 ∑ ∑

Pkl Bkl (u, v) = Au + Bv + C,

k=s−1 l=m−1

which implies that for 3 × 3 tensor-product submesh, the resulting IT-spline basis functions Bi (u, v) possess local partition of unity and local linear polynomials reproduction property. This result also indicates that for 3 × 3 tensor-product submesh, the IT-spline basis functions Bi (u, v), i = 1, 2, . . . , n give a second-order approximation to C 2 continuous functions. These imply the theorem. □ In general, however, the resulting IT-spline basis functions Bi (u, v), i = 1, 2, . . . , n may not form a partition of unity. In order to obtain IT-spline basis functions with exact partition of unity, we propose the following scheme to modify each of the IT-spline basis functions Bi (u, v), i = 1, 2, . . . , n as [ ] n ∑ Mi (u, v) = Bi (u, v) + wi (u, v) 1 − Bi (u, v) , (10) i=1

∑n 2 where the weight functions ∑n wi (u, v), i = 1, 2, . . . , n are C continuous and satisfy i=1 wi (u, v) = 1. It is of interest to note that if i=1 Bi (u, v) = 1 for (u, v) ∈ D, then from (10), we have Mi (u, v) = Bi (u, v) for (u, v) ∈ D. Remark 2. To guarantee the exact partition of unity, we may use the following common normalizing method ρi Bi (u, v) Ri (u, v) = ∑n . (11) i=1 ρi Bi (u, v) Then the resulting normalized functions Ri (u, v), i = 1, 2, . . . , n have compact supports and form a partition of unity. However, since ∑ the IT-spline basis functions Bi (u, v), i = 1, 2, . . . , n are not always positive within their own n supports, the sum i=1 ρi Bi (u, v) may take the value of zero in somewhere. These imply that the normalizing method (11) may be invalid as it involves rational form.

10

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

It is clear that the resulting modified IT-spline basis functions Mi (u, v), i = 1, 2, . . . , n, have the following properties (A) C 2 continuity. ( ) ( ) ( ) ∑n (B) Interpolation property. Since Bi u j2 , v j2 = δi j and i=1 Bi u j2 , v j2 = 1, we have Mi u j2 , v j2 = δi j , where (u j2 , v j2 ) is the knot coordinate of P j . (C) Linear independence. ∑n ∑n (D) Partition of unity. i=1 Mi (u, v) = i=1 wi (u, v) = 1. We now give some discussions on the choosing of the C 2 weight functions wi (u, v), i = 1, 2, . . . , n. For simplicity, we can choose the weight functions wi (u, v) as wi (u, v) =

1 , i = 1, 2, . . . , n, n

(12)

or { 1, i = k, wi (u, v) = 0, i = 1, 2, . . . , n, and i ̸= k.

(13)

Then the resulting modified IT-spline basis functions Mi (u, v), i = 1, 2, . . . , n are piecewise bi-quartic polynomial functions. By (12), for a large n, we can see that suppMi (u, v) ≈ suppBi (u, v) for all 1 ≤ i ≤ n, that is all the resulting modified IT-spline basis functions Mi (u, v) nearly have compact supports for a large n. By (13), one of the resulting modified IT-spline basis functions Mk (u, v) has global support and all the rest of the modified IT-spline basis functions have the same compact supports as that of the corresponding IT-spline basis functions. In order to insure that all the modified IT-spline basis functions Mi (u, v), i = 1, 2, . . . , n have the same compact supports as that of the corresponding IT-spline basis functions Bi (u, v), that is suppMi (u, v) = suppBi (u, v) for all 1 ≤ i ≤ n, from (10), it is sufficient to request that suppwi (u, v) = suppBi (u, v) for all 1 ≤ i ≤ n. To this end, we can choose the weight functions wi (u, v), i = 1, 2, . . . , n as the following classical C 2 rational T-spline blending functions ρi Fi (u, v) wi (u, v) = ∑n , (14) i=1 ρi Fi (u, v) where Fi (u, v) is the bi-cubic T-spline blending functions given in (6) with the knot vectors ui = [u i0 , u i1 , u i2 , u i3 , u i4 ] and vi = [vi0 , vi1 , vi2 , vi3 , vi4 ] inferred by Rule 2, the same as that of the IT-spline basis function Bi (u, v). With the weight functions wi (u, v) given by (14), it is obvious that suppwi (u, v) = suppBi (u, v) for all 1 ≤ i ≤ n and thus suppMi (u, v) = suppBi (u, v) for all 1 ≤ i ≤ n. After getting all the modified IT-spline basis functions Mi (u, v) corresponding to the control points Pi , i = 1, 2, . . . , n, an IT-spline surface S(u, v) can be generated by blending the control points Pi and the functions Mi (u, v), i = 1, 2, . . . , n, i.e., S (u, v) =

n ∑

Pi Mi (u, v) .

(15)

i=1

( ) Since Mi u j2 , v j2 = δi j , it is obvious that S (u i2 , vi2 ) =Pi for i = 1, 2, . . . , n. We now discuss the approximation property of the modified IT-spline basis functions over a hierarchical regular n T-mesh T . Let ℑ = span {Mi (u, v)}i=1 be the modified IT-spline space defined over a hierarchical regular T-mesh T with the weight functions wi (u, v), i = 1, 2, . . . , n given by (14). For any two neighboring cells in T , we assume that the level difference is at most one. Ω is the rectangular parameter domain occupied by T . Ωk is the parametric domain that consists of all the cells in level k. h k is the maximum length of the cells in level k. For a given continuous function f (u, v) ∈ C (Ω ), we would like to know the approximation error of f (u, v) from ℑ. With the weight functions wi (u, v), i = 1, 2, . . . , n given by (14), since the resulting modified IT-spline basis functions Mi (u, v), i = 1, 2, . . . , n have the properties of partition of unity and compact support, we have the following theorem.

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

11

Theorem 2. For any continuous function f (u, v) ∈ C (Ω ), there exists an IT-spline function g(u, v) ∈ ℑ such that ) ( √ | f (u, v) − g(u, v)| ≤ 2 3 2 + 1 ω ( f, h l ) , (u, v) ∈ Ωl , where w ( f, h) is the modulus of continuity of f , i.e. ω ( f, h) = classical ℓ2 -norm.

max | f (x) − f (y)|, with ∥·∥2 denoting the

∥x−y∥2 ≤h

Proof. For any function f (u, v) ∈ C (Ω ), we define g(u, v) =

n ∑

f (u i , vi ) Mi (u, v),

i=1

where (u i , vi ) are the knot coordinates in T corresponding to the modified IT-spline basis functions Mi (u, v), i = 1, 2, . . . , n. We choose the weight functions wi (u, v) given in (14). Let θ ∈ Ωl be a cell in level l. For any (u, v) ∈ θ, we have f (u, v) − g(u, v) =

n ∑

[ f (u, v) − f (u i , vi )] Mi (u, v)

i=1

=

∑

[ f (u, v) − f (u i , vi )] Mi (u, v)

i∈K

≤ max | f (u, v) − f (u i , vi )| , i∈K

where K = { i| Mi (u, v) ̸= 0, i = 1, 2, . . . , n, (u, v) ∈ θ }. Now suppose θ is contained in a cell θ ′ at level l − 1. Denote the corresponding tensor-product mesh of T at level l −1 as Tl−1 . Then there is a 3 × 3 tensor-product submesh H in Tl−1 whose central cell is θ ′ and the modified IT-spline basis functions at the vertices of H do not vanish in θ . By the local support properties of the modified IT-spline basis functions, we have {(u i , vi ) |i ∈ K } ⊂ H . For a hierarchical regular T-mesh T , generally we have h l ≥ h l−1 /2. These together, we get ) ( √ max | f (u, v) − f (u i , vi )| ≤ w f, 3 2h l−1 i∈K ) ( √ ≤ 3 2 + 1 ω ( f, h l−1 ) ) ( √ ≤ 2 3 2 + 1 ω ( f, h l ) . These imply the theorem. □ Fig. 5 shows the IT-spline basis function B3 (u, v) and the modified IT-spline basis function M3 (u, v) with the weight functions wi (u, v) given by (12), (13) and (14) respectively, all of which are associated to the control point P3 as shown in Fig. 4. In Fig. 5, for the weight functions (13), we have set k = 3 and for the weight functions (14), we have set all ρi = 1. Figs. 6 and 7 show the modified IT-spline basis functions Mi (u, v) with different weight functions and the corresponding C 2 IT-spline surfaces for the control points Pi marked with solid blue dots. The surfaces S1 (u, v), S2 (u, v), S7 (u, v) and S8 (u, v) are generated by using the weight functions (12). The surfaces S3 (u, v), S4 (u, v), S9 (u, v) and S10 (u, v) are generated by using the weight functions (13) with k = 1. And the surfaces S5 (u, v), S6 (u, v), S11 (u, v) and S12 (u, v) are generated by using the weight functions (14) with all ρi = 1. 3.3. Refinement algorithm and local adjustment algorithm We now∑ give some discussions concerning the refinement of the C 2 IT-spline surface. Given an IT-spline surface n S (u, v) = i=1 Pi Mi (u, v) over a regular T-mesh T , and Ω is the rectangular parameter domain occupied by T . After some refinements, we assume that the resulting regular T-mesh is T ∗ and the associated IT-spline surface is ∑ S ∗ (u, v) = mj=1 P j∗ M ∗j (u, v). The rectangular parameter domain occupied by T ∗ is the same as that of T , that is Ω .

12

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

Fig. 5. The IT-spline basis function B3 (u, v) and the modified IT-spline basis function M3 (u, v) with the weight functions wi (u, v) given by (12), (13) and (14), respectively.

For exact refinement, we require that S ∗ (u, v) = S (u, v) for all (u, v) ∈ Ω . For general refinement case, however, this exact refinement requirement for S ∗ (u, v) is not always possible. By approximation of exact refinement, we may require that S ∗ (u k , vk ) = S (u k , vk ) for all 1 ≤ k ≤ N , where (u k , vk ), k = (1, 2, . . . ), N are the sampling ∗ points of set the sampling points (u k , vk ) = u ∗k2 , vk2 for all 1 ≤ k ≤ m, ( ∗Ω and ) N ≥ m. For N = m, we may ∗ ∗ ∗ where u , v is the knot coordinate of P in the regular T-mesh T , then from the requirement conditions k k2 ) k2 ( ( ) ( ) S ∗ u ∗k2 , u ∗k2 = S u ∗k2 , u ∗k2 for all 1 ≤ k ≤ m, we can get Pk∗ = S u ∗k2 , u ∗k2 , k = 1, 2, . . . , m. For N > m, the requirement conditions S ∗ (u k , vk ) = S (u k , vk ) for all 1 ≤ k ≤ N form an over-determined linear system, then by using the classical least square method, we can obtain the optimal least square solution as follows ( )−1 P∗ = MT M MT S, [ ]T where P∗ = P1∗ , P2∗ , . . . , Pm∗ , S = [S (u 1 , v1 ) , S (u 2 , v2 ) , . . . , S (u N , v N )]T and ⎡ ⎤ M1∗ (u 1 , v1 ) M2∗ (u 1 , v1 ) · · · Mm∗ (u 1 , v1 ) ⎢ M ∗ (u 2 , v2 ) M2∗ (u 2 , v2 ) · · · Mm∗ (u 2 , v2 ) ⎥ 1 ⎢ ⎥ M=⎢ ⎥. .. .. .. .. ⎣ ⎦ . . . . ∗ ∗ ∗ M1 (u N , v N ) M2 (u N , v N ) · · · Mm (u N , v N )

13

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

Fig. 6. Modified IT-spline basis functions with different weight functions and the corresponding IT-spline surfaces for the control points marked with solid blue dots.

∑n We now give some discussions on the local adjustment of IT-surface S (u, v) = i=1 Pi Mi (u, v). By using the weight functions (14), with the changing of a control point Pi , the shape of the surface S (u, v) will be changed locally in the domain (u, v) ∈ suppMi (u, v). Without changing the control points, we can also adjust the shape of surface S (u, v) locally by using the positive weights ρi in the weight functions (14). † † † We further consider to adjust the surface points S(u l , vl ) to the targeted points Pl , l = 1, 2, . . . , L. For this purpose, we propose the following iteration adjustment algorithm to adjust the surface S (u, v) { † S0 (u, v) = S (u, v) , [ ( )] (16) † † † † † † † Sl (u, v) = Sl−1 (u, v) + Pl − Sl−1 u l , vl Bl (u, v) , l = 1, 2, . . . , L , †

†

†

†

†

†

†

where Bl (u, v) is an IT-spline basis function as defined in (9), with the knot vectors ul = [u l0 , u l1 , u l , u l3 , u l4 ] † † † † † † † and vl = [vl0 , vl1 , vl , vl3 , vl4 ] chosen in the way that Bl (u, v) vanishes at the knot coordinates of the points † † l−1 n {Pi }i=1 ∪ {Pk }l−1 k=1 (notice that for l = 1, we set the point set {Pk }k=1 as an empty set ∅). For 1 ≤ l ≤ L, after l † † † n iterations, from (16), it is obvious that Sl (u, v) interpolates the point set {Pi }i=1 ∪ {Pl }lk=1 and Sl (u, v) = S (u, v) † † l for any (u, v) ∈ Ω \ ∪k=1 {suppBk (u, v)}. And thus after L iterations, the resulting surface SL (u, v) interpolates † L † † n L the point set ( {Pi }i=1 ) ∪ {Pl }l=1 and SL (u, v) = S (u, v) for any (u, v) ∈ Ω \ ∪l=1 {suppBl (u, v)}. Apparently, if † † † † † Pl = Sl−1 u l , vl for all 1 ≤ l ≤ L, then SL (u, v) = S (u, v) for any (u, v) ∈ Ω .

14

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

Fig. 7. Modified IT-spline basis functions with different weight functions and the corresponding IT-spline surfaces for the control points marked with solid blue dots.

Fig. 8 gives an example of adjusting the IT-spline surface points to the targeted points, where the targeted points have been marked with red star points. 4. Applications 4.1. Fitting bivariate functions Given a bivariate function f (u, v) defined on [0, 1] × [0, 1], we first sample data points (u i , vi , f (u i , vi )), 2 i = 1, 2, . . . , n from it and then apply the ∑nC modified IT-spline basis functions Mi (u, v), i = 1, 2, . . . , n to obtain an interpolation function g(u, v) = i=1 f (xi , yi )Mi (u, v). We choose the weight functions wi (u, v) given in (14) with all ρi = 1. The difference between g(u, v) and f (u, v) is then measured by computing the following root mean square error between the function values on a dense grid √∑ ∑ [ ( ) ( )]2 M N i=0 j=0 f u i , v j − g u i , v j , RMS = (M + 1) (N + 1) where u i = i/M,v j = j/N ,M = N = 100. To normalize the error values, we divided the RMS error by the difference of the maximum and minimum values of f (u, v) over the domain. For convenience, we use DOF to abbreviate for degree of freedom (the number of the modified interpolation basis functions).

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

15

Fig. 8. Adjust the points on the surface S13 (u, v) to the targeted points (marked with red star points). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1 RMS and normalized RMS (NRMS) errors between test functions and the approximations. i

f1

f2

f3

f4

RMS NRMS

0.006727 0.005523

0.001621 0.007296

0.001875 0.005115

0.001259 0.003594

Example 1. In this experiment, we consider the following four test functions f i (u, v), i = 1, 2, 3, 4 used in [30] for bivariate scattered data interpolation [ ] [ ] (9u − 2)2 + (9v − 2)2 (9u + 1)2 (9v + 1)2 f 1 (u, v) = 0.75 exp − + 0.75 exp − − 4 49 10 [ ] 2 2 ] [ (9u − 7) + (9v − 3) + 0.5 exp − − 0.2 exp −(9u − 4)2 − (9v − 7)2 , 4 1 [tanh (9 − 9u − 9v) + 1] , 9 1.25 + cos(5.4v) f 3 (u, v) = , 6 + 6(3u − 1)2 √ [ ] 1 f 4 (u, v) = 64 − 81 (u − 0.5)2 + (v − 0.5)2 − 0.5, 9 f 2 (u, v) =

where (u, v) ∈ [0, 1] × [0, 1]. The approximate ranges of the four functions are f 1 (u, v) ∈ [0, 1.3] , f 2 (u, v) ∈ [0, 0.22] , f 3 (u, v) ∈ [0, 0.38] , f 4 (u, v) ∈ [0, 0.39], respectively, see [30]. For each test function f i (u, v), we derive an associated interpolation function gi (u, v). Fig. 9 shows the four interpolation functions gi (u, v), i = 1, 2, 3, 4. Table 1 demonstrates that the developed interpolation scheme reconstructs a test function quite accurately regardless

16

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

Fig. 9. Four resulting interpolation functions gi , i = 1, 2, 3, 4.

of the type of the function. It is clear from Table 1 that the normalized RMS error is less than 0.73 percent of the range of the function values. 4.2. Fitting open 3D meshes Given an open surface triangulation with vertices V j , j = 1, 2, . . . , m in 3D space, the corresponding parameter values (u j , v j ) of V j , j = 1, 2, . . . , m are obtained by the method of shape-preserving parametrization developed in [36], and the parameter domain is assumed to be [0, 1] × [0, 1]. In order to construct an IT-spline surface to fit the given surface triangulation, we need only to compute all the control points Pi , i = 1, 2, . . . , n associated with the modified IT-spline basis functions Mi (u, v), i = 1, 2, . . . , n. ∑ We choose the weight functions wi (u, v) given n in (14) with all ρi = 1. The fitting IT-spline surface is S(u, v) = i=1 Pi Mi (u, v). We compute the control points Pi , i = 1, 2, . . . , n in two steps. Firstly, we construct regular T-meshes level by level over the parameter domain [0, 1]×[0, 1] and determine the associated knot coordinates (u i2 , vi2 ), i = 1, 2, . . . , n. Secondly, for each 1 ≤ i ≤ n, we find the corresponding parameterized triangular ∆i′ where the knot coordinates (u i2 , vi2 ) locate in and determine the corresponding barycentric coordinates (xi , yi , z i ) of (u i2 , vi2 ) in ∆i′ . In the open surface triangulation, we have a triangular ∆i corresponding to the parameterized triangular ∆i′ , whose three vertices are denoted by Vi1 , Vi2 and Vi3 . Then the control points Pi are computed by Pi = xi Vi1 + yi Vi2 + z i Vi3 , i = 1, 2, . . . , n. Example 2. In this example, we fit cat head model by using the modified IT-splines basis functions Mi (u, v) with the weight functions wi (u, v) given in (14). We set all ρi = 1. The model has 135 points and 257 faces. Fig. 10 shows the numerical results. 4.3. Solving numerical PDEs n given in (9) to In this subsection, we provide three examples of using the IT-spline basis functions {Bi (u, v)}i=1 solve the Poisson’s equation with Dirichlet boundary conditions. We consider the standard Poisson’s equation with Dirichlet boundary conditions as follows { −∆µ = f, on Ω , (17) µ = 0, on ∂Ω . { } The weak form is to seek µ ∈ V := η ∈ H 1 (Ω ) , η = 0 on ∂Ω such that

a(µ, η) = F(η), ∀η ∈ V,

(18)

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

17

Fig. 10. Fitting cat head mesh model with modified IT-spline basis functions.

where a is a bilinear form and F is a linear functional defined by ∫ ∫ a(µ, η) = ∇µ · ∇ηdΩ , F(η) = f ηdΩ . Ω

Ω

n Among all the IT-spline basis functions {Bi (u, v)}i=1 , suppose Bi (u, v), i = 1, 2, . . ., m (m < n) satisfy Bi (u, v) |∂ Ω = 0. Then we define a new finite dimensional subspace Vh as follows

Vh = span {Bi (u, v) : Bi (u, v) |∂Ω = 0 , i = 1, 2, . . . , m} . The Galerkin projection replaces the infinite-dimensional space V by the finite dimensional subspace Vh ⊂ V . Then the discretization reads: Find µh ∈ Vh , such that a(µh , ηh ) = F(ηh ),

∀ηh ∈ Vh .

∑m The numerical approximation µh is constructed as a linear combination µh = i=1 ci Bi (u, v) with unknown real coefficients ci ∈ R, i = 1, 2, . . . , m. Upon inserting µh into the weak form and testing with ηh = Bi (u, v) for i = 1, 2, . . . , m, we can obtain the following linear system of equations Ac = b, with stiffness matrix A = (a(Bi , B j ))i, j=1,2,...,m , unknowns c = (c1 , . . . , cm )T , and the right-hand side vector b = (F (Bi ))i=1,2,...,m . In the following three examples, we suppose the physical domain Ω = [0, 1]2 . Example 3. In this experiment, we choose the following exact solution µ(x, y) = e10(x

2 −x)(y 2 −y)

− 1,

and f (x, y) is determined by (17). We compute the solution by using the IT-spline basis functions Bi (u, v). As a comparison, we also compute the solution by using the classical bi-cubic B-splines. Uniform refinement is performed to show the convergence rate.

18

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

Fig. 11. Numerical results with uniform refinement for Example 3. Table 2 Numerical results with uniform refinement for Example 3. DOF 9 81 361 1521 6241

IT-spline basis functions

Bi-cubic B-splines

∥µ − µh ∥ L 2

CR

COND

∥µ − µh ∥ L 2

CR

COND

1.596728e − 4 2.681097e − 6 1.531331e − 7 9.059169e − 9 4.607873e − 10

– 3.720057 3.831128 3.931916 4.219655

4.125992e+0 2.489946e+1 9.991575e+1 4.001353e+2 1.601051e+3

6.039291e − 3 9.290129e − 4 2.235017e − 4 5.481567e − 5 1.356340e − 5

– 1.703894 1.906686 1.954386 1.978507

7.702860e+0 1.501194e+1 2.157510e+1 8.716943e+1 3.509343e+2

Table 2 and Fig. 11 show the numerical results for Example 3. For convenience, we use COND to abbreviate for condition number of the stiffness matrix A. The convergence rate CR with respect to the L 2 -norm is computed by ( ) 2 log eh,k /eh,k−1 CR = , log (m k−1 /m k ) where eh,k and m k denote the L 2 -norm error ∥µ − µh ∥ L 2 and DOF at level k, respectively. From the numerical results as shown in Table 2, we can see that the IT-spline basis functions have better performances in the L 2 -norm error and CR, whereas the classical bi-cubic B-splines have better performances in the COND. Example 4. In this experiment, we choose the following exact solution x (1 − x) y (1 − y) µ(x, y) = , (x − 0.5)2 + (y − 0.5)2 + 0.01 and f (x, y) is determined by (17). We compute the solution by using the IT-spline basis functions Bi (u, v). Fig. 12 shows the numerical results with adaptive refinement. Example 5. In this experiment, we choose the following exact solution x (1 − x) y (1 − y) x (1 − x) y (1 − y) µ(x, y) = − , 2 2 (x − 0.25) + (y − 0.25) + 0.01 (x − 0.75)2 + (y − 0.75)2 + 0.01 and f (x, y) is determined by (17). We compute the solution by using the IT-spline basis functions Bi (u, v). Fig. 13 shows the numerical results.

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

19

Fig. 12. Numerical results with adaptive refinement for Example 4.

Fig. 13. Numerical results with adaptive refinement for Example 5.

4.4. Isogeometric analysis We consider again to solve the standard Poisson’s equation (17). Suppose now that the physical domain Ω is parameterized by a global geometry function G : (u, v) ∈ Ω0 = [0, 1]2 → (x, y) ∈ Ω defined as follows (x, y) = G (u, v) :=

n ∑

Pi Mi (u, v) ,

i=1

where Pi := (xi , yi ) ∈ R2 and Mi (u, v), i = 1, 2, . . . , n are the modified IT-spline basis functions with the weight functions wi (u, v) given by (12). And Ω0 is the parameter domain. n Among ⏐ all the IT-spline basis functions {Bi (u, v)}i=1 , suppose thath Bi (u, v), i = 1, 2, . . ., m (m < n) satisfy ⏐ Bi (u, v) G −1 (∂Ω) = 0. Then we define a finite dimensional subspace V as follows { } V h = span G i (x, y)| G i (x, y) = Bi ◦ G −1 , G i (x, y)|∂Ω = 0, i = 1, 2, . . . , m .

20

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

The isogeometric approximation of the weak form in (18) is given as: Find µh ∈ V h , such that a(µh , ηh ) = F(ηh ),

∀ηh ∈ V h .

The approximate solution µh can be written as µh (x, y) =

m ∑

ci G i (x, y) =

i=1

m ∑

ci Bi ◦ G −1 (x, y),

i=1

with unknown coefficients ci ∈ R, i = 1, 2, . . . , m. Thus the weak form solution of problem (18) is converted into solving the following linear system Lc = r, with stiffness matrix L = (a(G i , G j ))i, j=1,2,...,m , unknowns c = (c1 , . . . , cm )T , and the right-hand side vector r = (F (G i ))i=1,2,...,m , where ∫ ( ) ( ) ( ) a Gi , G j = ∇ Bi ◦ G −1 · ∇ B j ◦ G −1 d xdy ∫Ω [ ] [ ] DG(u, v)−T ∇ Bi · DG(u, v)−T ∇ B j |det DG(u, v)| dudv, = ∫Ω0 F(G i ) = f (x, y)Bi ◦ G −1 (x, y)d xdy Ω ∫ f (G(u, v)) Bi (u, v) |det DG(u, v)| dudv. = Ω0

and DG(u, v) is a 2 Jacobian matrix ] [ DG(u, v) = ⎡

∂x ∂u ∂y ∂u n ∑

⎢ ⎢ i=1 =⎢ n ⎢ ∑ ⎣ i=1

∂x ∂v ∂y ∂v

xi

∂ Mi (u, v) ∂u

yi

∂ Mi (u, v) ∂u

n ∑

⎤ ∂ Mi (u, v) ⎥ ∂v ⎥ i=1 ⎥. n ∑ ∂ Mi (u, v) ⎥ ⎦ yi ∂v i=1 xi

Example 6. In this example, we solve the Poissons equation on an L-shape physical domain Ω = [−1, 1] × [−1, 0] ∪ [−1, 0] × [0, 1] with zero Dirichlet boundary conditions. The equation has an exact solution µ(x, y) = sin(π x)sin(π y). The L-shape physical domain is parameterized by the modified IT-spline basis functions Mi (u, v) with the weight functions wi (u, v) given by (12). Fig. 14 shows the numerical results. Fig. 15 shows the convergence plots for Examples 4–6 concerning the L 2 -norm errors with respect to DOF. From the numerical results, the advantage of the adaptive refinement is distinct. 5. Conclusions We have proposed a new class of C 2 IT-spline basis functions. The resulting IT-spline basis functions are always linearly independent. Several new modified schemes for guaranteeing the partition of unity and compact support properties of the C 2 IT-spline basis functions are developed. By the new proposed local iteration adjustment algorithm, we can adjust the points on the IT-spline surface to the targeted points conveniently. We demonstrate several numerical examples in fitting bivariate functions and open 3D meshes, solving Poisson’s equation with Dirichlet boundary conditions and IGA. Experimental results show that the new IT-spline basis functions are suited for geometric design and IGA. In the future, we will concentrate on developing exact local refinement algorithm for the proposed C 2 IT-spline basis functions and explore more applications.

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

21

Fig. 14. Numerical results with adaptive refinement for Example 6.

Fig. 15. Convergence plots for Examples 4–6.

Acknowledgments We thank the anonymous reviewers for their valuable remarks for improvements. The first author thank Prof. Falai Chen, University of Science and Technology of China, for several fruitful and insightful discussions on the technical aspects of this work. The research is supported by the National Natural Science Foundation of China (Nos. 61802129, 11771453), the Natural Science Foundation Guangdong Province, China (No. 2018A030310381), the Postdoctoral Research Foundation of China (No. 2015M571931) and the Fundamental Research Funds for the Central Universities (No. 2017MS121). References [1] T.W. Sederberg, J. Zheng, A. Bakenov, A. Nasri, T-splines and TNURCCs, ACM Trans. Graph. 22 (3) (2003) 161–172. [2] T.W. Sederberg, D.L. Cardon, G.T. Finnigan, N.S. North, J. Zheng, T. Lyche, T-spline simplification and local refinement, ACM Trans. Graph. 23 (2004) 276–283. [3] H.W. Lin, Z.Y. Zhang, An efficient method for fitting large data sets using T-splines, SIAM J. Sci. Comput. 35 (6) (2013) 3052–3068. [4] M. Dörfel, B. Jüttler, B. Simeon, Adaptive isogeometric analysis by local h-refinement with T-splines, Comput. Methods Appl. Mech. Engrg. 199 (2010) 264–275. [5] Y. Bazilevs, V.M. Calo, J.A. Cottrell, J.A. 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. [6] 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. [7] 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.

22

Y. Zhu and X. Han / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112835

[8] X. Li, M.A. Scott, Analysis-suitable T-splines: characterization, refineability, and approximation, Math. Models Methods Appl. Sci. 24 (2014) 1141–1164. [9] L.B. Veiga, A. Buffa, D.C.G. Sangalli, Analysis-suitable T-splines are dual-compatible, Comput. Methods Appl. Mech. Engrg. 1 (2012) 249–252. [10] J.J. Zhang, X. Li, On the linear independence and partition of unity of arbitrary degree analysis-suitable T-splines, Commun. Math. Stat. 3 (2015) 353–364. [11] X. Li, J.J. Zhang, AS++ T-splines: LInear independence and approximation, Comput. Methods Appl. Mech. Engrg. 333 (2018) 462–474. [12] J.J. Zhang, X. Li, Local refinement for analysis-suitable++ T-splines, Comput. Methods Appl. Mech. Engrg. 342 (2018) 32–45. [13] H.M. Kang, F.L. Chen, J.S. Deng, Modified T-splines, Comput. Aided Geom. Design 30 (2013) 827–843. [14] X.D. Wei, Y.J. Zhang, L. Liu, T.J.R. Hughes, Truncated T-splines: Fundamentals and methods, Comput. Methods Appl. Mech. Engrg. 316 (2017) 349–372. [15] D. Forsey, R. Bartels, Hierarchical B-spline refinement, Comput. Graph. 22 (1988) 205–212. [16] C. Giannelli, B. Jüttler, H. Speleers, THB-Splines: The truncated basis for hierarchical splines, Comput. Aided Geom. Design 29 (2012) 485–498. [17] T. Dokken, T. Lyche, K.F. Pettersen, Polynomial splines over locally refined box-partitions, Comput. Aided Geom. Design 30 (2013) 331–356. [18] K.A. Johannessen, F. Remonato, T. Kvamsdal, On the similarities and differences between classical hierarchical, truncated hierarchical and LR B-splines, Comput. Methods Appl. Mech. Engrg. 291 (2015) 64–101. [19] J.S. Deng, F.L. Chen, Y.Y. Feng, Dimensions of spline spaces over T-meshes, J. Comput. Appl. Math. 194 (2006) 267–283. [20] J.S. Deng, F.L. Chen, X. Li, C.Q. Hu, W.H. Tong, Z.W. Yang, Y.Y. Feng, Polynomial splines over hierarchical T-meshes, Graph. Models 70 (2008) 76–86. [21] X. Li, J.S. Deng, F.L. Chen, Surface modeling with polynomial splines over hierarchical T-meshes, Vis. Comput. 23 (2007) 1027–1033. [22] L. Tian, F.L. Chen, Q. Du, Adaptive finite element methods for elliptic equations over hierarchical T-meshes, J. Comput. Appl. Math. 236 (2011) 878–891. [23] N. Nguyen-Thanh, H. Nguyen-Xuan, S.P.A. Bordas, T. Rabczuk, Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-dimensional elastic solids, Comput. Methods Appl. Mech. Engrg. 200 (2011) 1892–1908. [24] N. Nguyen-Thanh, J. Kiendl, H. Nguyen-Xuan, R. Wüchner, K.U. Bletzinger, Y. Bazilevs, T. Rabczuk, Rotation free isogeometric thin shell analysis using PHT-splines, Comput. Methods Appl. Mech. Engrg. 200 (2011) 3410–3424. [25] P. Wang, J.L. Xu, J.S. Deng, F.L. Chen, Adaptive isogeometrican analysis using rational PHT-splines, Comput. Aided Des. 43 (2011) 1438–1448. [26] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, Wiley, Chichester, 2009. [27] A. Embar, J. Dolbow, I. Harari, Imposing Dirichlet boundary conditions with nitsche’s method and spline-based finite elements, Int. J. Numer. Methods Engrg. 83 (2010) 877–898. [28] D.D. Wang, J.C. Xuan, An improved NURBS-based isogeometric analysis with enhanced treatment of essential boundary conditions, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2425–2436. [29] J.M. Escobar, J.M. Cascón, E. Rodríguez, R. Montenegro, A new approach to solid modeling with trivariate T-splines based on mesh optimization, Comput. Methods Appl. Mech. Engrg. 200 (2011) 3210–3222. [30] S.Y. Lee, G. Wolberg, S.Y. Shin, Scattered data interpolation with multilevel B-splines, IEEE Trans. Vis. Comput. Graph. 3 (1997) 228–244. [31] R.J. Zhang, W.Y. Ma, An efficient scheme for curve and surface construction based on a set of interpolatory basis functions, ACM Trans. Graph. 30 (2011) Article 10. [32] X.Y. Yuan, W.Y. Ma, Parametric mesh regularization for interpolatory shape design and isogeometric analysis over a mesh of arbitrary topology, Comput. Methods Appl. Mech. Engrg. 284 (2015) 906–942. [33] X.L. Han, A class of general quartic spline curves with shape parameters, Comput. Aided Geom. Design 28 (2011) 151–163. [34] M. Antonelli, C.V. Beccari, G. Casciola, A general framework for the construction of piecewise-polynomial local interpolants of minimum degree, Adv. Comput. Math. 40 (2014) 945–976. [35] H. Akima, A new method of interpolation and smoth curve fitting based on local procedures, J. Assoc. Comput. Mach. 17 (1970) 598–602. [36] M.S. Floater, Parameterization and smooth approximation of surface triangulations, Comput. Aided Geom. Design 14 (1997) 231–250.