{'fraction_non_alphanumeric': 0.05230284444194356, 'fraction_numerical': 0.013926792898742369, 'mean_word_length': 4.142361111111111, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Constructing Trivariate B-splines with Positive Jacobian by Pillow Operation and Geometric Iterative Fitting

The advent of isogeometric analysis has prompted a need for methods to generate Trivariate B-spline Solids (TBS) with positive Jacobian. However, it is difficult to guarantee a positive Jacobian of a TBS since the geometric pre-condition for ensuring the positive Jacobian is very complicated. In this paper, we propose a method for generating TBSs with guaranteed positive Jacobian. For the study, we used a tetrahedral (tet) mesh model and segmented it into sub-volumes using the pillow operation. Then, to reduce the difficulty in ensuring a positive Jacobian, we separately fitted the boundary curves and surfaces and the sub-volumes using a geometric iterative fitting algorithm. Finally, the smoothness between adjacent TBSs is improved. The experimental examples presented in this paper demonstrate the effectiveness and efficiency of the developed algorithm.


Introduction
While the traditional geometric design community focuses on the design of curves and surfaces, the advent of Isogeometric Analysis (IGA) [1] has made the development of methods for designing Trivariate B-spline Solids (TBSs) imperative. In IGA, a valid TBS should have a positive Jacobian value at every point in its domain. A negative Jacobian value at any points in the domain of the TBS can render the IGA invalid.
As is widely recognized, the generation of a valid TBS is difficult, owing to two key factors: The geometric condition for ensuring that a TBS is valid, i.e., the Jacobian value at every point is positive, is highly nonlinear. Hence, it is theoretically very difficult to guarantee the validity of a TBS. The state-of-the-art methods for generating TBS usually transform the validity problem into an optimization problem. However, the optimization is prone to fail owing to the high nonlinearity of the objective function. Therefore, to generate valid TBSs, the high nonlinearity of the objective function should be reduced and an efficient method should be developed for solving the optimization problem.
The regions with negative Jacobian are usually concentrated near the boundary curves where adjacent surfaces are smoothly stitched (refer to Fig. 4). To reduce negative Jacobian values, the input model should be segmented along the boundary curves using the pillow operation [2], a method originally developed for improving the quality of hexahedral meshes.
In this paper, we developed a method that can generate a TBS with a guaranteed positive Jacobian. Using a tetrahedral (tet) mesh with six surfaces segmented on its boundary mesh as the input, we first partition the tet mesh model into seven sub-volumes using the pillow operation. After each of them is parameterized into a cubic parameter domain, seven initial valid TBSs are constructed. Moreover, starting with the initial valid TBSs, the boundary curves, boundary surfaces and the TBSs are fitted by a geometric iterative fitting algorithm, known as the Geometric Feasible Direction (GFD) algorithm. In each iteration of the GFD algorithm, the movements of the control points are restricted inside a feasible region to ensure the validity. Finally, the smoothness between two adjacent TBSs is improved by the GFD algorithm. In this way, the validity of the generated TBSs is guaranteed with desirable smoothness between adjacent TBSs.
The structure of this paper is as follows: Section 1.1 provides an overview of the additional work performed to develop this study. Section 2 presents an overview of the proposed method. In Section 3, we introduce the validity conditions for B-spline curves, surfaces, and TBSs, and the geometric continuity between TBSs. In Section 4, the Geometric Feasible Direction (GFD) algorithm is developed. After elucidating the details of the developed method in Section 5, the experimental results are illustrated in Section 6. Lastly, Section 7 concludes the paper.

Related work
In this section, work related to TBS modeling and geometric iterative fitting is briefly reviewed.
TBS modeling: To analyze the arterial blood flow by IGA, a trivariate NURBS-solid modeling the artery was constructed using a skeleton-based method [3]. Following volume parameterization by a harmonic function, a cylinder-like trivariate B-spline solid was generated with a singular centric curve [4]. Moreover, trivariate B-spline solids with positive Jacobian values were produced from boundary representations using optimization-based approaches [5,6]. However, the optimization method may fail when the objective function is highly nonlinear. Based on the given boundary conditions and guiding curves, a NURBS solid was constructed to model a swept volume by the variational approach [7]. Moreover, trivariate T-spline solids are employed in IGA owing to their adaptive refinement capability [8]. To fill a genuszero triangular mesh, a T-spline solid was constructed after mesh untangling and smoothing [9]. Starting from a boundary surface triangulation with genus-zero, Zhang et al. developed a mapping-based method to generate rational trivariate solid Tsplines [10]. Furthermore, to fill a boundary triangulation with arbitrary genus, a polycube mapping-based algorithm for constructing T-spline solids was proposed in [11]. In [12], a volumetric T-spline was constructed for filling a B-rep model using Boolean operations, polycube mapping, and octree subdivision. The inputs for the aforementioned methods are triangular meshes. When the input was a genus-zero T-spline surface model, Zhang et al. constructed a solid T-spline whose boundary exactly conformed to the given T-spline surface model [13].
On the other hand, Catmull-Clark subdivision solids were developed to model the computational domain in IGA [14]. Additionally, other representations for spline solids include simplex splines [15] and polycube splines [16,17].
The methods described above usually generate a trivariate solid to fill a given B-rep model. However, the generation of a TBS by fitting a tet mesh model is much easier than by filling a B-rep model because it is very easy to produce a tet mesh using popular software such as, TetGen [18], and NetGen [19]. Hence, it is feasible to generate a TBS by fitting a tet mesh model. In Ref. [20], a tet mesh model is fitted by the geometric iterative method to generate a TBS. However, there are some regions with negative Jacobian values close to the boundary. In this experiment, a tet mesh is first segmented into seven subvolumes, each of which is fitted with a TBS by a geometric iterative fitting method. In this way, the generated TBSs are ensured to be valid, i.e. the Jacobian value at any point of the TBSs is positive.
In conclusion, the representations of a trivariate solid include B-spline solid, NURBS solid, T-spline solid, subdivision solid, etc. The input for a solid generation method is either a boundary representation, such as B-spline surface, a triangular mesh, or a volume mesh, such as a tet mesh. In our method, the input model is segmented to generate the valid TBS. While a tet meshes can be partitioned into sub-volumes, the boundary representations, as solids, are hard to be segmented owing to the loss of inner information. Therefore, tet meshes are taken as inputs for our method.
Geometric iterative fitting: Geometric Iterative Fitting (GIF), also called Progressive-Iterative Approximation (PIA), was first developed in [21,22]. The GIF method endows iterative methods with geometric meanings, facilitating the handling of geometric problems appearing in the field of geometric design. It was proved that the GIF method is convergent for Bspline fitting [23,24], NURBS fitting [25], T-spline fitting [26], subdivision surface fitting [27][28][29], as well as curve and surface fitting with totally positive basis [22]. Moreover, the GIF meth-ods have been employed in some applications such as reverse engineering [30,31], curve design [32], surface-surface intersection [33], etc. In Ref. [20], the GIF method was employed to fit a tet mesh model with a TBS, showing that the GIF method is capable of solving singular linear systems. In this paper, the GIF method is used to fit several sub-volumes with valid TBSs.

Algorithm overview
The TBS generation algorithm is demonstrated in Fig. 1. The input to this algorithm is a tet mesh with six boundary surfaces ( Fig. 3(a)). The tet mesh is first parameterized into a cubic parametric domain. Then, the cubic domain is segmented into seven sub-domains, as shown in the ( Fig. 3(b)). By mapping the subdomains into the input tet mesh, it is divided into seven subvolumes. After the initial TBSs are constructed, the boundary curves, boundary surfaces, and the seven TBSs are respectively fitted using the GFD algorithm. Finally, the smoothness between adjacent TBSs is improved. In subsequent sections, the details of the TBS generation algorithm will be elucidated.

Validity conditions and geometric continuity between TBSs
In this paper, what we want to generate is a composition of valid TBSs with as desirable as possible smoothness between adjacent TBSs. So, the validity conditions and geometric continuity definition between TBSs should be clarified.

Validity conditions
In this section, the validity conditions for a B-spline curve, surface, and TBSs are presented, respectively.

Definition 1. 1) A B-spline curve is valid, if its derivative
vector is nonzero at any parameter in its domain; 2) A B-spline surface is valid, if its normal vector is nonzero at any point in its parametric domain; 0 at any point in its parametric domain.
Given a B-spline curve with control point P i , i = 0, 1, · · · , m, and denoting the difference vectors as, the valid condition for the B-spline curve is, Proposition 1 (Validity condition for B-spline curves). A Bspline curve of degree d is valid if the apertures of the minimum circular cones enclosing the difference vectors {T i , T i+1 , · · · , T i+d−1 }, i = 0, 1, · · · , m − d, respectively, are all less than π 2 . The proof of Proposition 1 is evident.
Moreover, suppose we are given a B-spline surface of degree d u × d v with control points S i j , i = 0, 1, · · · , m, j = 0, 1, · · · , n, and denote the difference vectors as and, Let M I J , I = 0, 1, · · · , m − d u , J = 0, 1, · · · , n − d v be the subcontrol-polygon constituted by the control points, Moreover, suppose U I J and V I J are the unit axis vectors of the minimum circular cones C u I J and C v I J enclosing the difference vectors T u i j (1) and T v kl (2) of the sub-control-polygon M I J , starting from the apexes of cones, respectively. Then, a sufficient condition for the validity of a B-spline surface is presented as follows: Proposition 2 (Validity condition for B-spline surfaces). If T u i j · U I J > T u i j · V I J ≥ 0, and, T v kl · V I J > T v kl · U I J ≥ 0, where T u i j and T v kl are defined on each sub-control-polygon M I J , I = 0, 1, · · · , m − d u , J = 0, 1, · · · , n − d v , the B-spline surface is valid.
Proof: Because, we have, It means that, all of the vectors T u i j × T v kl of the sub-controlpolygon M I J are on the same side of the plane with U I J × V I J as the normal vector. It makes the normal vector of the given Bspline surface nonzero at any point in its domain. So, the given B-spline surface is valid. 2 Remark 1. Note that, the condition T u i j · U I J > T u i j · V I J ≥ 0, and, T v kl · V I J > T v kl · U I J ≥ 0, means that the two minimum circular cones C u I J and C v I J are separate. According to the definition of normal vector, the closer to orthogonality the two unit axis vectors U I J and V I J of the circular cones C u I J and C v I J , the better the validity of the B-spline surface.
Similarly, we can develop a sufficient condition for determining the validity of a TBS , and B d w k (w) are B-spline basis functions of degrees d u , d v , and d w , respectively, and, are control points. Denote the difference vectors as Moreover, letting be the sub-grid constituted by the control points we have,

Proposition 3 (Validity condition for TBSs).
If Moreover, because the local support property of TBS, the Jacobian value J(u, v, w) is determined by one of the sub-grid

Geometric continuity between TBSs
Given two adjacent TBSs (Fig. 2), and where, P i p , j p ,k p and Q i q , j q ,k q are control points, , and B k q (ω) are B-spline basis, the geometric continuity between the two TBSs P(u, v, w) and Q(µ, ν, ω) is defined as, Definition 2 (Geometric Continuity). The two TBSs P(u, v, w) and Q(µ, ν, ω) are G n geometric continuous along the common boundary surface P(u, v, 1) = Q(µ, ν, 0), if, 1) the two B-spline surfaces P(c, v, w) and Q(c, ν, ω) are G n geometric continuous along the common boundary curve P(c, v, 1) = Q(c, ν, 0), and, 2) the two B-spline surfaces P(u, d, w) and Q(µ, d, ω) are G n geometric continuous along the common boundary curve where c and d are arbitrary constants in their domains (refer to Fig. 2).

Geometric Feasible Direction algorithm
As pointed out in Section 2, the composition of TBSs is a step-by-step process, in the order of boundary curve fitting, boundary surface fitting, TBS fitting and smoothness improvement. The fitting problems and smoothness improvement in these steps are all modeled as constraint minimization problems, where P i = (x i , y i , z i ), i = 0, 1, · · · , m are the unknown control points, and E is an objective function. Moreover, the conditions for ensuring the validity of curves, surfaces, and TBSs are taken as the hard constraints (8). The specific minimization problems will be elucidated in subsequent sections.
To solve the constraint minimization problem (8), we developed the geometric feasible direction (GFD algorithm (Algorithm 1). GFD has clear geometric meanings, so it is a kind of geometric iterative fitting method. As listed in Algorithm 1, the inputs to GFD algorithm are an initial B-spline curve, a Bspline surface, or a TBS with control points arranged in a onedimensional sequence {P (0) i , i = 0, 1, · · · , m}, and a constrained minimization problem (8), including an objective function and the constraints.
Specifically, in the k th iteration, the gradient vector of the objective function E (8) is, where In the iteration, each control point P (k) i moves along the directional vector, 4 to produce the new control point P (k+1) i , i.e., where τ i ∈ [0, 1] is a weight. Initially, we set τ i = 1, i = 0, 1, · · · , m, and the direction D (k) in the k th iteration is composed of, Then, the weights τ i , i = 0, 1, · · · , m are adjusted one by one to make the direction D (k) (11) a feasible one, which satisfies the following two requirements: 1) The new B-spline curve, surface, or TBS, generated by replacing , satisfies the constraints of the minimization problem (8), To this end, we discretize τ i ∈ [0, 1] to {0, 1 n , 2 n , · · · , n n }, and select a weight as large as possible τ i , which satisfies the two aforementioned requirements 1) and 2). In our implementation, we take n = 20, and δ d = 0.3.
Moreover, let, After the feasible direction D (k) is figured out, we calculate the Armijo step α k (please refer to Ref. [35] for the calculation of the Armijo step). Then, the new control points are generated as, The GFD algorithm is presented in Algorithm 1, where the termination condition is taken as, The value of ε e will be specified in solving specific problems. The convergence analysis of the GFD algorithm is demonstrated in the Appendix.

Valid TBS generation by the geometric iterative method
In Section 2, the overview of the method for generating valid TBSs was presented. In this section, we will elucidate the details of the method.

Partition of the tet mesh model by pillow operation
As stated above, the input to our method is a tet mesh model with six surfaces segmented on its boundary mesh (Fig. 3(a)). The scaled Jacobian value of the generated TBS is heavily influenced by the segmentation of the boundary mesh. If two adjacent surfaces are C n (n ≥ 1) continuous along their common boundary, or even concave along the boundary, the Jacobian values at the points near the boundaries of the fitting TBS will be small or even negative (Figs. 4, 5(a)). This makes the generated TBS invalid. Fig. 4illustrates the distribution of the  scaled Jacobian values of the TBS fitting the tet mesh model in Fig. 3(a) by the method developed in Ref. [20]. It can be seen that, the small or negative Jacobian values mainly concentrate in the region (in blue) close to the boundaries of the surfaces, along which two adjacent surfaces are smoothly stitched. The quality of the generated TBS can be improved by the pillow operation [2], which was originally invented to improve the quality of a hexahedral mesh. Refer to Fig. 5(a), there is a quadrilateral ABCD, and the Jacobian value at vertex A is 0. After performing the pillow operation ( Fig. 5(b)), the Jacobian values at all of the vertices are greater than 0. In order to perform the pillow operation on the input tet mesh, the tet mesh model (Fig. 3(a)) is first parameterized into a cubic parameter domain Ω = [0, 1]×[0, 1]×[0, 1] by the volume parameterization method [20] (Fig. 3(b)). Then, a sub-domain is inserted in the parameter domain Ω. As illustrated in Fig. 3(b), the vertices of the cubes Ω and Ω c are denoted as v 0 , v 1 , · · · , v 7 , and v c 0 , v c 1 , · · · , v c 7 , respectively. Connecting v i to v c i , i = 0, 1, · · · , 7 generates six sub-domains Ω u , Ω d , Ω l , Ω r , Ω f , and Ω b . For example, the sub-domain Ω u is enclosed by the six faces Mapping the seven sub-domains into the original tet mesh model produces seven partitioned sub-volumes (Fig. 3(c)).

Construction of the initial TBSs
In Section 5.1, the input tet mesh model is partitioned into seven sub-volumes. Each of them are parameterized into the cubic parameter domain by the volume parameterization method developed in [20]. Note that the parameterization on the common boundary curves and common boundary surfaces of adjacent sub-volumes should conform with each other. Moreover, each cubic parameter domain is sampled into a (M + 1) × (N + 1) × (K + 1) grid. Similar to the parameterization, the grid on the common boundary curves and common boundary surfaces of adjacent sub-volumes should be the same. Mapping the grids into the corresponding sub-volumes leads to the control grids of the initial TBSs, whose knot vectors are uniformly distributed in [0, 1] × [0, 1] × [0, 1] with Bézier end conditions. Owing to the conformity of the parameterization and control grids on the common boundary curves and boundary surfaces, there is the unique control grid on each boundary curve or boundary surface.
The geometric iterative fitting method, which is employed to fit the sub-volumes (refer to Section 5.3), requires that the initial TBSs should be valid. If there is an initial TBS that is not valid, the grid sampling rule should be adjusted, or the initial TBS should be optimized to produce the valid initial TBS. In all of the examples we handled, the initial TBSs generated by the method aforementioned are all valid.

Geometric iterative fitting
The result generated by the developed method is a composition of seven valid TBSs. As stated above, the objective function for guaranteeing the validity of a TBS is highly nonlinear with a high number of unknowns, so the optimization is prone to fail. Even if we can find a solution, the computation for solving the optimization problem is complicated, owing to a significant number of unknowns. To reduce the difficulty in guaranteeing the validity of the TBSs, we solve this problem step by step, in the order of, As mentioned above, the first objective we want to reach is that the generated TBSs should be valid, that is, the Jacobian value at any point of each TBS should be greater than 0. This means that the boundary curves and boundary surfaces of TBSs should also be valid (refer to Section 3.1 for validity condition).

Boundary curve fitting
Suppose the data point sequence to be fitted is, By volume parameterization [20], each data point V j has been assigned a parameter t j ∈ [0, 1], j = 0, 1, · · · , M c . The data point sequence (13) will be fitted by a B-spline curve, where, P i are the control points, B k i (u) are the B-spline basis functions with degree d, and the knot vector of the B-spline curve is uniformly defined in [0, 1] with Bézier end conditions. The curve fitting should concern two aspects. One is the fitting error; and the other is the validity of the fitting curve. On one hand, the fitting error can be modeled as, On the other hand, denote α as the aperture of the minimum circular cone C enclosing the difference vectors, with unit axis vector, Evidently, if P (u) = 0 (14), the B-spline curve is not valid. On the contrary, the larger the norm P (u) , the better is the validity of the B-spline curve (14). With fixed lengths of the vectors T i in (15), the smaller the aperture α, the larger is the norm P (u) . Therefore, the improvement of validity of the fitting curve can be achieved by minimizing the following objective function, In conclusion, the fitting of the data point sequence (13) with a valid B-spline curve (14) can be formulated as the following minimization problem, 1) The control polygon of P(u) satisfies Proposition 1, 2) Compatibility condition, where λ c ∈ [0, 1] is a weight. Constraint 1) in the minimization problem (16) ensures that the fitted curve satisfies the validity condition for B-spline curves (Proposition 1). The values of λ c taken in our implementation are listed in Table 3. Moreover, the B-spline curves are constructed piece by piece. After one piece of B-spline curve is constructed, it is fixed. When constructing a new B-spline curve P(u), it should satisfy the compatibility condition 2) with the adjacent constructed B-spline curves, which belong to the same B-spline surface as P(u). That is, they should satisfy Proposition 2. Specifically, refer to Fig. 6(a), when constructing the new B-spline curve P(u) with control points P 00 , P 10 , · · · , P d u ,0 (red points), if there exists the piece of B-spline curve adjacent to P 00 , with control points P 00 , P 01 , · · · , P 0,d v (blue points), (the two pieces of B-spline curves belong to the same surface,) the difference vectors (refer to Eqs. (1) and (2)) T u 00 , T u 10 , · · · , T u d u −1,0 and T v 00 , T v 01 , · · · , T v 0,d v −1 should satisfy Proposition 2. The constrained minimization problem (16) is solved by the GFD algorithm (Algorithm 1), where the inputs are the initial B-spline curves constructed in Section 5.2, and the objective function and constraints of the minimization problem (16). In our implementation, the threshold ε e in the termination condition (12) is taken as 10 −4 .

Boundary surface fitting
After the boundaries of the seven sub-volumes are fitted by B-spline curves, they are fixed. The next task is to fit the boundary surfaces of the sub-volumes. Each of the boundary surfaces is a triangular mesh with four boundary curves. Suppose the inner vertices of such a triangular mesh are V i , i = 0, 1, · · · , M s , each of which is assigned a pair of parameter values (u i , v i ) by the volume parameterization [20]. Moreover, suppose the Bspline surface fitted to the triangular mesh is, (17) where, S i j are the control points, B d u i (u) and B d v j (v) are B-spline basis functions with degrees d u and d v , respectively, both defined in [0, 1] with Bézier end conditions. Similar to the curve fitting, the surface fitting also involves two aspects:the fitting precision, and the validity. The fitting precision is measured by the following energy function, Moreover, the validity of the B-spline surface should be improved as per the requirement. Suppose C u and C v are the minimum circular cones enclosing the difference vectors T u i j (1) and T v i j (2), respectively, with unit axis vectors T u and T v , taken as, The validity of a B-spline surface can be improved by two strategies. Firstly, as analyzed in Section 5.3.1, the apertures of the two cones should be as small as possible, which is achieved by minimizing the following two energy functions, Secondly, the two cones should be as close to perpendicularity as possible, which is formulated as the minimization of the following energy function (Eq. (19)), Therefore, the B-spline surface fitting the triangular mesh with validity guarantee can be generated by the following minimization problem, where λ s , µ s ∈ [0, 1] are weights, whose values can be found in Table 3, and Constraint 1) in Eq. (20) ensures that the B-spline surface S(u, v) (17) satisfies the validity condition for B-spline surfaces (Proposition 2). Similar to the curve case, the B-spline surfaces are constructed piece by piece. Hence, the construction of a new Bspline surface S(u, v) (the surface in light red in Fig. 6(b)), should satisfy the compatibility condition 2) (Eq. (20)) with existing B-spline surfaces (the two surfaces in light blue in Fig. 6(b)), which are adjacent to one of the corners of the new B-spline surface S(u, v), and belong to the same TBS with S(u, v). In other words, the corresponding sub-control-polygons of these surfaces adjacent to one corner of S(u, v) should satisfy Proposition 3. Specifically, refer to Fig. 6(b), supposing the degree of the TBS is d u × d v × d w , the three sub-control-polygons with size d u × d v , d v × d w , and d u × d w , respectively (refer to Eq. (3)), should satisfy the validity condition for TBSs, i.e., Proposition 3.
The constrained minimization problem (20) is solved by the GFD algorithm (Algorithm 1). The inputs to the GFD algorithm are the initial B-spline surfaces constructed in Section 5.2, the objective function and constraints of the minimization problem (20). According to Eq. (17), the control points of the B-spline surface are arranged in a one-dimensional sequence with lexicographic order, i.e., {S 00 , S 01 , · · · , S mn }. It should be pointed out that the threshold ε e in the termination condition (12) is set as 10 −5 .

TBS fitting
Now, it is the time to fit seven sub-volumes by TBSs. Similar to the boundary curve fitting and boundary surface fitting, the seven sub-volumes are also fitted one by one. Each sub-volume is a tetrahedral (tet) mesh, whose boundary triangular meshes were fitted with B-spline surfaces (refer to Section 5.3.2), and thus fixed in the TBS fitting. Suppose V n , n = 0, 1, · · · , M h are the inner vertices of a sub-volume mesh with parameter values (u n , v n , w n ), and the TBS fitting the sub-volume mesh vertices is, where H i jk are the control points, B d u i (u), B d v j (v), and B d w k (w) are the B-spline basis functions of degrees d u , d v , and d w , defined on the interval [0, 1], with Bézier end conditions, respectively.
Similar to the boundary curve and surface fitting, the generation of a TBS should involve the following factors: 1) The fitting precision to the tet mesh vertices, and, 2) the improvement of validity of each TBS.
On one hand, the fitting precision is modeled by the following formula, On the other hand, to improve the validity of a TBS, we define the difference vectors along the u−, v−, and w−directions, respectively, as, Suppose C u , C v , and C w are the minimum circular cones enclosing T u i jk , T v i jk , and T w i jk , with unit axis vectors T u , T v , and T w , taken as, respectively, Similar as the analysis in Section 5.3.2, the smaller the apertures of the cones, the better the validity of the TBS, which is modeled as the minimization of the following three energy functions, i.e., where, N u , N v , and N w are the number of the vectors T u i jk , T v i jk , and T w i jk , respectively. In addition, the closer the perpendicularities of the three cones are to each other, the better is the validity of the TBS, which can be formulated as the minimization of the following energy functions: Essentially, the minimization problem for generating the valid TBSs can be modeled as, 2) The boundary control points of H(u, v, w) (21) are fixed, where λ h , µ h ∈ [0, 1] are weights, and Constraint 1) makes the TBS H(u, v, w) satisfy the validity condition for TBSs (Proposition 3). The values of weights λ h , µ h ∈ [0, 1] are presented in Table 3.
The constrained minimization problem (24) is solved by the GFD algorithm (Algorithm 1), with the initial TBSs constructed in Section 5.2 as input. The control points of an input TBS (refer to Eq. (21)) are arranged in a one-dimensional sequence with lexicographic order, i.e., {H 000 , H 001 , · · · , H m u ,m v ,m w }. Moreover, we set ε e = 10 −6 in the termination condition (12).

Smoothness and fairness improvement
In this section, we will improve the smoothness between two adjacent TBSs, and the fairness of each TBS. Suppose H l (u, v, w), l = 0, 1, · · · , 6 is a TBS with control points H l,i jk .
Since the boundary paths are fixed before the TBS fitting, the two adjacent TBSs have reached G 0 continuity. According to Proposition 4 and Fig. 2, the two adjacent TBSs H l 1 (u, v, w) and H l 2 (u, v, w) are G 1 continuous, if the two unit vectors T l 1 ;i j and T l 2 ;i j , i.e., have the same direction. Then, the smoothness between a pair of adjacent TBSs can be improved by the minimization of the energy function, The smoothness between all of the adjacent TBSs can be improved by minimizing the energy, (1 − T l 1 ;i j · T l 2 ;i j ).
Meanwhile, the fairness of the TBSs is improved by minimizing the fairness energy, The fairness energy E h f air is a quadratic function of control points of the TBS H l . In our implementation, the composite trapezoidal rule [36] is utilized to calculate the integrals in the fairness energy E h f air . Therefore, the improvement of smoothness and fairness with validity guarantee can be modeled as the constrained minimization problem, where, λ f ∈ [0, 1] is a weight, listed in Table 3. Constraint 1) in (25) guarantees that each TBS satisfies the validity condition for TBSs (Proposition 3). The constrained minimization problem (25) is solved by the GFD algorithm (Algorithm 1). The input to the GFD algorithm is the composition of the seven TBSs, constructed in Section 5.3.3, with control points arranged in a one-dimensional sequence. The threshold ε e in the termination condition (12) is taken as 10 −5 .

Results
The algorithm developed in this paper is implemented using Visual Studio C++ 2010, and run on a PC with Intel Core2 Quad CPU Q9400 2.66 GHz and 4GB memory. In this paper, we demonstrated a method to generate six TBS models. In Fig. 7, the cut-away views of the TBS models are illustrated. It can be seen that the iso-parametric curves vary smoothly not only inside a single TBS, but between two adjacent TBSs as well. Moreover, in Fig. 8, the distribution of the scaled Jacobian values [37] of the TBSs is visualized with different colors.
The darker the red color, the higher the scaled Jacobian values. As shown in Fig. 8, the scaled Jacobian value of these models are all positive. Most region of the six TBS models is in red. Actually, in each of the six TBS models, the Jacobian values are larger than 0.5 in over 80% region. The statistics of TBSs generated by our method are listed in Table 1. The third column of Table 1 is the number of control points of TBSs, reading as m × n × k × l. It means that, the number of control points of the central TBS is m × n × k, those of the six surrounding TBSs are m × n × l, m × n × l, m × k × l, m × k × l, n × k × l, and n × k × l, respectively. The fourth column is the fitting error, which is defined as, where N i is the number of tet mesh vertices in the i th subvolume (i = 0, 1, · · · , 6), and L is the diagonal length of the bounding box of the whole model. The fitting errors for all of the six models are in the order of magnitude 10 −2 . Moreover, in the fifth column, the ratios of the volume of the region with scaled Jacobian in (0, 0.2] to the whole volume of the TBS are listed. It can be seen that the largest ratio (for the model Ball joint) is 1.96%, and the ratios for the other models are all below 0.5%, meaning that regions with scaled Jacobian in (0, 0.2] in generated TBSs are very small. In the last four columns, we list the time (in second) cost in curve, surface, TBS fitting, and quality improvement (i.e., smoothness and fairness improvement). The whole process for the six models takes time ranging from 11 min to 37 min.
Moreover, for comparison, the minimum scaled Jacobian values and average Jacobian values of TBSs generated by the method in Ref. [20] and our method, respectively, are presented in Table 2. The average Jacobian value is calculated as, where J(x, y, z) is the scaled Jacobian value [37] at (x, y, z). It can be seen from Table 2 that, while the minimum scaled Jacobian values of TBS models generated by the method in [20] are all negative, those of TBS models produced by our method are all positive. In addition, with the two methods, the average scaled Jacobian values are comparable.  (25)), in order to balance the energy functions well, and achieve desirable results. The weights employed in generating the six TBSs are presented in Table 3.
Limitations: One limitation of the developed method lies that, the segmentation manner, i.e., the pillow operation, is not desirable for any tet mesh models. For example, because the bottom of the Isis model shown in (Fig. 9) is slim, the quality of the TBS generated by our method is not very good, though the Jacobian value can be guaranteed positive. The fitting precision can just reach the order of magnitude 10 −1 , and the Jacobian values in 16.69% region lie in the interval (0, 0.2]. Therefore, more reasonable tet mesh segmentation method should be developed in the future.  2 Number of the control points of the TBSs generated by our method. 3 The ratio of the volume of the region with scaled Jacobian in (0,0.2] to the whole volume of the TBS. 4 Time (in second) cost in curve, surface, TBS fitting and quality improvement (smoothness and fairness improvement).

Conclusion
In this paper, we developed a method to generate a TBS by fitting a tet mesh model. The input to this method is a tet mesh model with six surfaces segmented on its boundary mesh. To improve the Jacobian values in the regions close to the boundaries, the tet mesh model was first partitioned into seven subvolumes using the pillow operation. Then, a geometric iterative fitting algorithm was developed to fit the boundary curves, boundary surfaces, and the sub-volumes, separately. Because the validity conditions are integrated into the geometric iterations, the Jacobian values of the fitting boundary curves, boundary surfaces, and TBSs are guaranteed to be positive. Therefore, the generated TBSs are valid for isogeometric analysis. Theorem 1. If ∇E(P 0 , P 1 , · · · , P m ) (8) is uniformly continuous on the region D ∩ L, and the objective function E(P 0 , P 1 , · · · , P m ) (8) is bounded, the GFD algorithm (Algorithm 1) is convergent.
Proof: Denote θ k as the angle between the feasible direction D (k) and −∇E (k) . In the iteration of the GFD algorithm, the following inequality holds, So, there exists some µ > 0, such that, Therefore, based on Lemma 1, together with that the objective function E(P 0 , P 1 , · · · , P m ) (8) is bounded, the GFD algorithm (Algorithm 1) is convergent. 2