T-Splines Based Isogeometric Topology Optimization with Arbitrarily Shaped Design Domains

: In this paper, a new isogeometric topology optimization (ITO) method is proposed by using T-splines based isogeometric analysis (IGA). The arbitrarily shaped design domains, directly obtained from CAD, are represented by a single T-spline surface which overcomes the topological limitations of Non-Uniform Rational B-Spline (NURBS). The coef ﬁ cients correlated with control points are directly used as design variables. Therefore, the T-spline basis functions applied for geometry description and calculation of structural response are simultaneously introduced to represent the density distribution. Several numerical examples show that the proposed approach leads to a coherent work ﬂ ow to handle design problems of complicated structures. The optimized results are free of checkerboard patterns without additional stabilization and ﬁ ltering techniques due to the proper-ties of T-splines, which also simpli ﬁ ed the post-processing. In addition, through performing local re ﬁ nement, we can easily achieve multiresolution optimization and in ﬁ ll optimization within the T-splines based framework. In general, the proposed method provides a possibility to design, analyze, and optimize engineering structures in a uniform model, which has the potential to improve design ef ﬁ ciency and reduce the cost of product development.


Introduction
Topology optimization (TO), which aims at finding the optimal distribution of materials in a prescribed design domain under certain constraints [1], has been gaining importance in recent years. Since the pioneering work of Bendsøe et al. [2], various impressive TO methods have been developed in the past three decades, such as solid isotropic material with penalization (SIMP) method [3,4], level set method (LSM) [5][6][7], evolutionary method [8] and moving morphable components (MMC) approach [9,10]. These methods have been applied to a wide range of applications [11][12][13][14][15][16] to obtain superior performance and cost-saving configurations of engineering structures, especially in the field of additive manufacturing (AM) [17][18][19] as TO may best exploit the design freedoms offered by AM.
In most of the topology optimization related issues, FEM has been used for the calculation of structure response and sensitivity analysis, which involves mesh generation procedure and results in heavy extra workload because mesh reparation is tiresome and unavoidable on most occasions. These tedious operation steps arising from different model description methods in modeling and optimization are expected to be surmounted. Additionally, the lower-degree continuity (C 0 ) of the unknown physical field in traditional FEM may lead to numerical instabilities [20] and affect the optimization accuracy. Thus, it becomes an old dilemma to establish a tight link between analysis models and accurate design models. Isogeometric analysis, proposed by Hughes et al. [21,22], has been proved to be a promising way out of this situation. IGA is a generalization of FEM which directly employs the basis functions used in representing the geometric model in CAD system to calculate the structural response. Because of its higher accuracy and compatibility with available FEM frameworks, IGA has been applied in a variety of fields [23][24][25][26].
In the last decade, there has been a growing interest in integrating IGA with topology optimization, termed as isogeometric topology optimization (ITO). The first work might be proposed by Seo et al. [27,28], where trimmed spline surfaces are used to treat the complex topology. Their optimization model that contains trimming curves and spline surface is consistent with CAD systems, which avoids additional model conversions. However, their sensitivity formulations may lead to a heavy computational burden. This kind of approach has also been applied in the optimization of shell structures [29]. A control points based SIMP method with material distribution approximated by NURBS basis functions was proposed by Hassani et al. [30] in 2012. Later, Dedè et al. [31] developed a phase field representation of the density distribution in TO with IGA applied to solve the phase field problems. Qian [32] presented a B-splines based method where the design space is restricted to the B-spline space. Their work proved that the splines based density distribution could provide an intrinsic filter. The efficiency [33] and robustness [34] have been discussed as well. IGA was also introduced into the LSM based topology optimization by Wang [35], where NURBS was used to represent the level set function. They also proposed an approach to achieve geometric constraints based on trimmed elements [36]. This kind of ITO method has been applied to solve stress problems [37], flexoelectric materials [38], and vibrating structures [39]. Besides, multi-resolution [40], multi-material [41,42], and multi-scale [43] topology optimization problems using IGA were also discussed. Recently, IGA was integrated into evolutionary based TO [44] and MMC based methods [45][46][47][48]. Gao et al. [49] proposed an ITO method for continuum structures using the Shepard function and NURBS basis functions. Based on their work, a design framework of auxetic metamaterials was developed [50]. To improve the efficiency of ITO, Wang [51] proposed a new framework including multilevel mesh, multigrid conjugate gradient method, and local-update strategy. More applications of ITO can be found in a recent review [52].
As we can see, most of the previous studies have focused on tensor product splines like B-spline and NURBS. The excellent mathematical and algorithmic properties, such as partition of unity, non-negativity and linear independence of NURBS benefit exceedingly the relevant methods. However, NURBS is severely limited by its tensor product topological structure. It is of great difficulty to represent complex engineering models such as airplanes or automobiles without using multiple patches and trimming features in the framework of NURBS [53,54], which brings the well-known problems of the gap and overlapping in CAD field. This limitation results in the fact that many available methods can just deal with topology optimization problems defined in simple geometries such as rectangular design domains. To overcome the defects of NURBS, T-spline was introduced into the CAD community in 2003 [55]. As a generalization of NURBS, T-splines possess a flexible topological structure by introducing T-juctions and extraordinary points. A trimmed NURBS surface can be converted to an untrimmed T-spline surface, and multiple NURBS patches are able to be merged into one single watertight T-spline patch [56]. In addition, T-splines have the ability to achieve sufficient local refinement while maintaining exact geometry [57]. All of these characteristics make T-splines an ideal technology for isogeometric analysis, which was described in detail in Bazilevs et al. [58,59]. We believe that the advantages of T-spline based isogeometric analysis can also be extended to be used in topology optimization. The topological flexibility makes T-splines based modeling especially attractive to be used for complex TO problem, so this is a worthy topic and has the probability to facilitate further integration of design and analysis. Therefore, in this paper, a new TO method based on the combination of isogeometric analysis and T-splines is proposed. The arbitrarily shaped design domain is represented by a single T-spline surface and the T-spline coefficients correlated with the control points are directly used as design variables that contribute to a continuous and smooth density distribution function (DDF). The Bézier extraction for unstructured T-splines is applied to simplify the computation and implement the isogeometric analysis. Hence, the same T-spline basis functions are used to represent not only the complex geometries but also the unknown field variables in analysis and the material distribution in topology optimization. Compared with NURBS based ITO, the proposed method is easier to handle design domains with complex geometric structures. The feature of local refinement of T-splines can obviously reduce the cost of computation. Moreover, the final result of the proposed method is described by a T-spline surface, which means that the representations of the input model and the output model keep consistent. To summarize, this approach can be applied to engineering structural design of arbitrarily shaped design domains and is likely to offer an alternative option for the development of new optimization tools.
The rest of the paper is organized as follows: In Section 2, we give a brief introduction to unstructured T-splines and Bézier extraction. Section 3 is dedicated to describing the new isogeometric topology optimization method, including T-splines based DDF, optimization formulation, sensitivity analysis, and the post-processing method. Section 4 shows three numerical examples, including complex structures, multiresolution, and infill optimization problems, which demonstrate the characteristics of the proposed method. Conclusions and future work are presented in Section 5.

T-splines Overview
This section provides a brief introduction of fundamental concepts underlying unstructured T-splines. More details are presented in Sederberg et al. [55,57,60]. In what follows, we will use the notation p for the polynomial degree, d s for the number of spatial dimensions, and d p for the number of parametric dimensions. We will focus on bi-cubic T-spline surfaces because of their predominance in industry.

The Unstructured T-mesh
A T-spline surface is constructed from a T-mesh which is a two-dimensional polygonal mesh that contains faces, edges and vertices. As an example shown in Fig. 1, each face in T-mesh is quadrilateral which is also referred to as element in the literature. Small circles identify the vertices of the T-mesh. Every vertex is associated with a weight x 2 R and a control point P 2 R d s , the gathering of which composes the control grid of T-spline surface.
In Fig. 1(a), T-junctions are marked with blue circles which are analogous to hanging nodes in conventional finite element method, while extraordinary points are denoted by red circles. We define extraordinary points as those interior vertices that are not T-junctions and whose valence does not equal four, namely, the number of edges that touch these interior vertices is not four. Then, an unstructured T-mesh can be defined as a T-mesh that contains extraordinary points, just as the T-mesh in Fig. 1(a).
The one-ring neighborhood of a vertex is the set of T-mesh faces that touch that vertex (yellow elements in Fig. 1(a)). The two-ring neighborhood involves the one-ring neighborhood and the T-mesh faces that touch the one-ring neighborhood (green elements in Fig. 1(a)). An n-ring neighborhood may be formed recursively in the same way [60].
We call the T-mesh faces and vertices contained in the two-ring neighborhood of an extraordinary point irregular faces and irregular vertices, respectively, while other T-mesh faces and vertices are called regular faces and regular vertices.

T-spline Blending Functions
To define T-spline blending functions, knot intervals (non-negative real numbers) should be assigned to each edge in the T-mesh. A valid knot interval configuration requires that the knot intervals on opposite sides of every element sum to the same value. Fig. 1(b) shows a valid knot interval configuration associated to the T-mesh in Fig. 1(a).
Then, each T-mesh vertex v A is associated with a T-spline blending function N A , where A 2 {1,…,n} and n is the global number of control points. Each N A is constructed according to the local grid structure and local knot intervals.
If v A is a regular vertex, N A can be defined as: where N 3 [Δξ A ](ξ) and N 3 [Δη A ](η) are cubic B-spline basis functions associated with the knot interval vectors Δξ A and Δη A , respectively. The local knot interval vectors can be obtained by the ray-intersection method proposed in [55]. The main process of the method is to march through the T-mesh in each topological direction, until p − 1 vertices or perpendicular edges are intersected. If a T-mesh boundary is crossed before p − 1 vertices or edges are encountered, zero-length knot intervals are appended to complete the knot interval vector, which creates an open knot vector structure along the boundary of the T-mesh. An example is shown in Fig. 2(a), the local knot interval vectors of v 1 and v 2 are given by :  In addition, Δξ A and Δη A also decide the support of N A , namely, the region that N A is non-zero. For example in Fig. 2(a) the support of N 1 is shown in pink and the support of N 2 in cyan, while these two blending functions are visualized in Fig. 2 However, if v A is an irregular vertex, it becomes difficult to define N A with acceptable continuity around extraordinary points. As the ray-intersection method can not be used, different studies have been put forward to handle this trouble [60][61][62]. Among these methods, Bézier extraction has been demonstrated as a promising approach in the application of isogeometric analysis. It is a great advantage that Bézier extraction provides a unified framework to calculate all the blending functions defined over an unstructured T-mesh, which means that we can utilize it to deal with not only the irregular vertex but also the regular vertex mentioned above.

Bézier Extraction
The aim of Bézier extraction is to generate the Bernstein-Bézier representation of the T-spline blending functions over each element. Given a vector N e ¼ N e a È É n e a¼1 , where N e a are the T-spline blending functions with support on the element e, n e is the number of these local blending functions and a is the local blending function index, there exists a linear operator C e such that N e ðn; gÞ ¼ C e Bðn; gÞ; (2)  where Λ(i, j) = (p + 1)(j − 1) + i, and B p i ðnÞ is the univariate Bernstein basis function defined over À1; þ1 ½ which can be calculated by Noted that the interval À1; þ1 ½ is chosen rather than the conventional unit interval 0; 1 ½ in the computer-aided geometric design (CAGD) community because the bi-unit interval can facilitate the Gaussian quadrature in isogeometric analysis.
Additionally, we use the IEN array [63] to store the relation between the e-th element local index a and the global index A, which can be stated as In Eq. (2), the linear operator C e is called the Bézier extraction operator of element e, which provides a simple representation of the T-splines blending functions. The computation of C e is the main process of Bézier extraction, which has been thoroughly explained in Scott et al. [60,64], so we will not go into details. However, it should be noted that the extraction procedure is different between regular and irregular elements. While the extraction of the regular element is based on knot insertion [64], the irregular element needs a linear interpolation scheme with an optimization process that assures G 1 continuity around the extraordinary points [60].
During the Bézier extraction, the elemental T-mesh should be generated, because some faces in the Tmesh are covered by some T-spline blending functions that can not be mapped to a single Bézier element. By adding all face extensions into the T-mesh, as shown in Fig. 3(a), all T-spline blending functions are guaranteed to be C ∞ within the faces of the elemental T-mesh. Meanwhile, each face corresponds to a Bézier element in physical space as shown in Fig. 3(c). And these Bézier elements make up a Bézier mesh in Fig. 3(b).
We want to note that Bézier extraction does not mean that we divide the T-spline surface into many Bézier surfaces separately, although the control points of these Bézier surfaces can be calculated easily using the Bézier extraction operator. We still make use of the T-spline control points and blending functions, which guarantee a higher continuity across elements. The essential lies in simplifying and unifying the representation and calculation of all the T-spline blending functions defined over an unstructured T-mesh.

The T-spline Geometrical Mapping
We define a mapping S e : e ! e from the parametric space n; g f g of the bi-unit Bézier element e to the physical space as shown in Fig. 3. Then given a control grid P i 2 e and the weight of each control point ω i , the T-spline geometrical mapping S e (ξ, η) can be given as: where w e ¼ x e a È É n e a¼1 is the element weight vector, W e = diag(w e ) is the diagonal weight matrix, and P e ¼ P e a È É n e a¼1 is the element control point matrix. We call R e n; g ð Þ ¼ R e a n; g ð Þ È É n e a¼1 the vector of element rational T-splines blending functions, which can be stated as Using the Bézier extraction Eqs. (2) and (7), we can calculate R e as R e n; g The derivatives of R e with respect to the local parametric coordinate ξ is @R e @n ¼ W e C e @ @n B n; g ð Þ W e n; g ð Þ ¼ W e C e 1 W e n; g ð Þ @B n; g ð Þ @n À B n; g ð Þ ½W e n; g ð Þ 2 @W e n; g ð Þ @n ! ; where W e n; g ð Þ ¼ w e ð Þ T C e B n; g ð Þ, and @R e @g can be obtained in the same way.

Analysis Suitability
Although we can construct free-form surface with arbitrary topology using T-splines in CAD, not every T-mesh is suitable for isogeometric analysis. A T-mesh without restriction may generate T-splines with linearly dependent blending functions, while linear independence is a requirement for successful analysis. Consequently, we have to utilize a mildly restricted subset of T-splines called analysis-suitable (AS) T-splines which was proposed by Li et al. in 2012 [65] to guarantee the linear independence. It was proven in [65] that a T-mesh without extraordinary points is analysis suitable if no horizontal and vertical T-junction extensions intersect with each other. In 2018, analysis-suitable++ (AS++) T-splines [66] was introduced to relax these restrictions, which can handle some special situations that the blending functions are linearly independent while some T-junction extensions intersect. For unstructured T-meshes that have extraordinary points, an analysis suitable T-spline is defined in [60] to be one whose T-mesh satisfies the condition in [65] and has no one-bay face extension spans an element in the three-ring neighborhood of an extraordinary point, and no extraordinary point lies within the three-ring neighborhood of another extraordinary point. To ensure the analysis suitability, we restrict ourselves within the set of analysis suitable T-spline.

Isogeometric Topology Optimization
In this section, we describe the new topology optimization method where T-splines are used to define the density distribution, and isogeometric analysis is applied to calculate the related structural response instead of the finite element method. The minimum compliance problem is considered only for the sake of numerical simplicity but without losing any generality.

T-spline Based Density Distribution
In order to construct the density distribution in topology optimization, a T-spline coefficient ρ i 2 [0, 1] is assigned to each control point of the design domain which is modeled by T-splines. In such case, the element density distribution v e n; g ð Þ can be stated as v e n; g ð Þ ¼ ðq e Þ T R e n; g ð Þ where q e ¼ q e a È É n e a¼1 is the element T-spline coefficient vector and R e n; g ð Þ is the blending function vector computed from Eq. (8).
This method can be understood as we add a new dimension to the 2D T-spline design domain. Therefore, it becomes a 3D T-spline surface that would be modified in each iteration, as shown in Fig. 4. Due to the regularity and non-negativity of T-spline rational blending functions, v e n; g ð Þ is always between 0 and 1, which meets the requirements of the SIMP method.
Following the SIMP interpolation approach, the Young's modulus E e (χ e (ξ, η)) according to the element density distribution is expressed by E e ðv e ðn; gÞÞ ¼ E min þ ½v e ðn; gÞ s ðE 0 À E min Þ; where E 0 is the solid material stiffness, E min is the minimal material stiffness used to avoid singularity in stiffness matrix, and s is the penalization power.

Isogeometric Analysis
For linear elastic problems, the element displacement field u e is approximated as u e n; g ð Þ ¼ q e ð Þ T R e n; g ð Þ ¼ X n e a¼1 R e a n; g ð Þq e a ; where q e a denotes the displacement at the IEN(a, e)-th control point. In the Galerkin's IGA method, the discrete equilibrium equation for elastic problems under static load may be expressed as where K is the stiffness matrix, U is the displacement vector, and F is the load vector associated with the control points. The stiffness matrix K is generated by assembling the local stiffness matrix K e . In the physical space Ω e , the element stiffness matrix K e is calculated by where B is the strain-displacement matrix, and D is the elastic tensor matrix. The integral for K e is pulled back onto the bi-unit Bézier element by where J is the Jacobian matrix of the T-spline geometrical mapping S e and can be expressed as J ¼ @x @n @y @n @x @g @y @g where the partial derivatives can be computed using Eqs. (6) and (9).
For two dimensional cases, the strain-displacement matrix B is formulated as ; where R e a is the a-th rational T-spline blending function of element e. Considering plane stress problems, the elastic tensor matrix D is given as where E e is the Young modulus corresponding to the density distribution and ν is the Poisson's ratio.
Substituting E e (χ e ) from Eq. (11) into Eq. (18), the elastic tensor, the elastic tensor matrix may be rewritten as where D 0 is the elastic tensor matrix for unit material.
The Gauss quadrature method is employed to achieve the numerical computation of the element stiffness matrix as follows where n gp is the number of the Gauss quadrature points (ξ i , η i ) and ω i is the corresponding quadrature weight. Note that the density value χ e is evaluated at each Gauss quadrature points to compute the Young modulus, which accounts for a more accurate approximation of the T-splines based density distribution than the element-wise constant way.

Topology Optimization Formulation
Given an arbitrarily shaped design domain Ω represented by T-splines, ρ i is assigned to each control point to generate the T-splines based density distribution v n; g ð Þ as in Section 3.1. Then the topology optimization formulation is developed based on χ, with the isogeometric analysis applied as in Section 3.2 to calculate the structural responses. Thus, for minimum compliance design problem, the formulation can be given as: Find : q i ði ¼ 1; 2; . . . ; nÞ Min : cðvÞ ¼ F T U S:t: : Z e v e n; g ð Þd e À c ⩽ 0 where ρ i is the T-spline coefficients which work as design variables, c is the objective function defined by the structural mean compliance, g is the global volume constraint; K is the stiffness matrix, U is the displacement vector, and F is the load vector; V 0 ¼ R d is the volume of the design domain, N e is the number of elements in Ω, and γ is the allowable volume fraction. The T-spline coefficient is limited to change between 0 and 1.

Sensitivity Analysis
To obtain the first-order derivative of the objective function c with respect to the design variables ρ i , the adjoint method [1,32] can be used. From Eq. (21), the sensitivity of compliance can be firstly given as where λ is the adjoint vector that can be computed from Here, the derivative of the stiffness matrix with respect to the design variables can be computed on an element level. We firstly substitute Eq. (19) into Eq. (15) to rewrite the element stiffness matrix as follows Then, the first-order derivative of the element stiffness matrix K e with respect to the design variables ρ i may be written as @K e @q i ¼ Z e sR e i ½v e sÀ1 ðE 0 À E min ÞB T D 0 BjJjd e : From Eq. (21), the first-order derivative of the volume constraint g with respect to the design variables ρ i can be obtained in a similar way as follows

Implementation
The method of moving asymptotes (MMA) [67] is used to solve the optimization formulation due to its superior characteristics and its wide application in structural optimization problems. Thus this continuous optimization problem is solved iteratively by a gradient-based optimization scheme. Fig. 5 shows the flowchart of the proposed topology optimization framework. As aforementioned, the input data contains a design domain constructed by T-splines and some essential parameters such as the prescribed volume fraction γ, the penalization power s and so on. Then, Bézier extraction is used to build the Bézier mesh and simplify the computation as described in Section 2.2. Next, the density distribution is constructed with the initialized design variables ρ i . At each iteration, the unknown structural responses are computed by the T-splines based IGA. The values and sensitivities of the objective and constraint functions are calculated to update the density distribution until satisfying the convergence criterion. Finally, the post-processing is applied to obtain the optimal T-splines model, which is described in the next section.

Post-processing
The aim of post-processing is to interpret the topology optimization results for further applications. It is of great difficulty to redesign or make use of the optimized structure if it can not be transformed into a CAD model [68][69][70]. Comparing to the conventional density-based topology optimization method that often results in a non-smooth discrete structure in the form of 0-1 elements, the proposed method has a more efficient and concise post-processing scheme due to the flexibility of T-splines.
It needs two steps to interpret the result of the proposed topology optimization method. 1) Detect the isocontour: In order to extract an appropriate structure from the optimized density distribution, a heuristic criterion is usually introduced in isogeometric topology optimization method, which can be mathematically expressed as [49] 0 ⩽ vðn; gÞ < v c void vðn; gÞ ¼ v c boundary v c < vðn; gÞ ⩽ 1 solid;

<
: (27) where χ c is a threshold and the boundary can be viewed as an isocontour of the density distribution.
We can see that Eq. (27) is similar to the level set representation of structures in LSM [5,71]. Additionally, it has been discussed in [49] that if the optimized densities are distributed near 0 or 1, it will be relatively suitable to set χ c = 0.5. The marching squares algorithm is a general method of detecting an isocontour. However, in the proposed method, there is an alternative way that is more suitable for CAD software. As shown in Fig. 6(a), the isocontor can be extracted by computing the intersection curves between a horizontal plane and the density distribution represented by a T-spline surface. The density value of the horizontal plane should be set to χ c .
2) Trim the T-spline design domain: To obtain the optimal structure, we can directly trim the T-spline design domain along the isocontour. Then the result will be represented by a single watertight Tspline surface, just like how we describe the arbitrarily shaped design domain. A detailed description of the T-splines trimming method is presented in [56]. Apart from this, some commercial CAD softwares can help to realize this operation on T-spline surface, such as Rhinoceros 3D and Autodesk Fusion 360, as shown in Fig. 6(b). In all examples, the design domain is represented by a bi-cubic T-spline surface. We design these arbitrarily shaped design domains using the commercial CAD software Rhinoceros 3D with the Autodesk T-spline plugin. The Young's modulus is considered as E 0 = 1 and E min = 1.0e − 9. The Poisson's ratio ν is set to be 0.3. The penalization power s is defined as 3. The imposed point load magnitude is assumed as 1. The 3 × 3 Gauss quadrature rule is used for each Bézier element. The convergence criterion is the relative difference of the objective function values between two iterations, which is set to be 1.0e − 4. The postprocessing threshold χ c is 0.5. Noted that no extra filtering schemes are contained in any example. The objective function values and total iteration number are denoted by "Obj" and "Iter", respectively.

Cantilever Beam with a Circular Hole
The first numerical example is an optimization problem of a cantilever beam with a circular hole. The structural design domain with the loads and boundary conditions is defined in Fig. 7(a) where L = 15, H = 10, and R = 3. As can be seen, the left side is fixed, and a concentrated force is vertically loaded at the low-right corner. The allowable volume fraction γ is set to be 40%, and the initial value of all design variables ρ i is equal to 0.4. Fig. 7(b) shows the T-spline surface that represents this structure. Noted that the T-spline model is a single watertight surface that contains 20 extraordinary points around the circular hole, while it requires to merge at least four patches of surfaces in NURBS based method as shown in Fig. 7(c).
The optimized results are presented in Fig. 8 including 2D view of the optimized density distribution in Fig. 8(a) and 3D view in Fig. 8(b). It can be seen that crisp boundaries are obtained due to the continuity of T-splines and the penalization method. Then, Fig. 8(c) shows the optimized structure after post-processing. Just like how we describe the geometry of the design domain, the final topology is represented by a single T-spline surface that can be designed and analyzed by IGA directly, as shown in Fig. 8(d). The volume fraction of the final T-spline surface is equal to 40.22%, which is just a little different from the prescribed volume fraction γ = 40%. In addition, the compliance value of the optimized structure is also recalculated. We get a value of 31.7746, which is lower than the optimized objective function 33.2217. Therefore, the post-processing scheme and the choice of χ c is appropriate. Finally, the convergence history of the objective function is displayed in Fig. 9 with some intermediate results of the density distribution. From the convergent curve, we can see that the proposed method gives a rapid and clear convergence. During the iterative process, it can be clearly seen that a three-dimensional T-spline surface is evolved gradually to control the topology of the two-dimensional structure.

Triangle Bracket
In this example, a design problem of a triangle bracket model is studied in detail to illustrate the ability of the proposed method to manage arbitrarily shaped design domains. Fig. 10(a) shows the complicated geometric structure of the triangle bracket design domain. As we can see, a horizontal force is applied, The allowable volume fraction γ is taken as 40%. It is a difficult task to model this complex design domain with tensor product NURBS patches without any surface trimming, while T-splines method can handle this easily, as shown in Fig. 10(b). The optimized density distribution is presented in Fig. 10(c), which is featured with sufficient continuity. And Fig. 10(d) displays the final T-spline surface. It can be seen that the optimized structure keeps the boundaries precisely. Additionally, the volume fraction of the final topology is equal to 40.06%, which demonstrates the accuracy of the post-processing scheme again.
Then, the influence of the initial design is discussed by applying three different initial density settings, as shown in Fig. 11. The first (left) column in Fig. 11 presents the three different initial density values that we investigate, namely, 0.2, 0.4 and 0.8. The second (middle) and the third (right) columns display the optimized density distribution in 2D and 3D views, respectively. As we can see, the results in the three cases are identical for the most part, while the optimized objective functions are close to each other as well. What is more, Fig. 12 shows the convergence histories of the three cases. It can be seen that all of the three convergent curves can reach similar convergence value, although different initial objective functions are given, and the beginning parts of the convergent curves are steep. Thus, to a certain extent, the proposed method is not sensitive to the initial densities.
The influence of the T-mesh that we utilized to construct the T-splines based design domain is also studied in this example, as shown in Fig. 13. It can be seen that we use three different T-meshes, which are generated by applying local refinement in different parts of Fig. 10(b). The optimized density distributions in 2D and 3D view are displayed in the second (middle) and the third (right) columns of Fig. 13 respectively. We can see that the refinement parts lead to more crisp structures in all of these three results. However, some small structures have appeared in Figs. 13(b) and 13(e), which means that the mesh dependency may occur in the proposed method. This phenomenon has been discussed in many works that focus on the isogeometric topology optimization based on NURBS or B-splines [32,49,72]. There are many ways to improve the performance of the proposed method and suppress the mesh dependency problems, such as controlling the orders of basis or applying the Shepard function [49]. In this example, we illustrate that performing multiresolution (or multilevel mesh) topology optimization [40,73] maybe one of the applicable ways to eliminate the mesh dependency from the proposed method. As can be seen from Fig. 14, the same T-splines model is utilized as the design domains (dark cyan surfaces), while three T-meshes with different refinement strategies are used to calculate the structure responses (yellow surfaces). The projection from the density distribution of the design domain to the Gauss quadrature points in the analysis domain can be obtained with high accuracy, because these Tspline surfaces precisely describe the same geometry, although they are built from distinct T-meshes. We can see that the optimized structures are almost the same, and the optimized objective functions are close to each other, which means that, in the proposed method, it is possible to refine the analysis model locally to achieve more accurate calculation without causing mesh dependency problems.

Multiple Loads Design Problem
In this section, a two-step optimization problem of a simply supported beam with multiple loads is studied to present the excellent flexibility of the proposed method. Fig. 15(a) displays the rectangular design domain, upon which five vertical loads with the same magnitude are applied individually. The structural sizes are defined as L = 100 and H = 50. Then, a T-spline surface with 2048 Bézier elements is given to represent this structure, as shown in Fig. 15(b). The objective function is the average of the five compliances corresponding to the five individual loads.
In the first step of the optimization process, the global volume fraction γ is set to be 50%. Figs. 16(a) and 16(b) shows the optimized density distribution in 3D and 2D view, respectively. Similar to the previous two examples, the density distribution has distinct interfaces between the solids and voids. The T-spline surface of the final topology in this step is shown in Fig. 16(c). We can see that ten internal holes appear after the first step optimization. It is noted that the volume fraction of Fig. 16(c) is 50.74%, which is very close to the prescribed volume fraction. Then, to demonstrate that the optimized structure can be modified or redesigned easily, in the second step an infill optimization is directly applied to the optimized result of the previous step, namely, the model in Fig. 16(c). The infill optimization that we realize is based on [74], the main idea of which is to impose constraints on local volume fractions in order to regulate the local material distribution and generate bone-like porous structures. Details can be found in Appendix A.  Fig. 17 shows the problem setup of the second step. As can be seen, the loads and boundary conditions remain unchanged, while the design domain becomes the T-spline surface obtained by the first step optimization. Note that the design domain has been refined many times from the T-spline surface in Fig. 16(c), because infill optimization requires considerable design variables to generate elaborate porous structures. The local volume fraction α and the local influence radius R e are taken as 65% and 2.0, Right (c, f, i): The optimized density distribution in 3D view with the optimized objective function and the total iterative step respectively, which means that in every local region with a radius of two the proportion of the solid area is restricted to no more than 65%.
The optimized density distribution is presented in Fig. 18(a). And Fig. 18(b) displays the final topology. As we can see, compared with the global volume constraint that generates large solid and empty parts, the infill optimization creates a structure dominated by crossing elongated sub-structures, which looks like the internal porous structure of a bone in nature. Fig. 19 gives the convergence history of the optimization process in the second step and some intermediate results. It can be seen that the density distribution evolves towards porous structures gradually, which demonstrates the ability of the T-splines based   Bézier elements, which is refined from the optimized structure of the first step method to describe very complicated material distribution. However, the convergence of the infill optimization is not as fast as the traditional problems of global volume constraints, with some fluctuations on the convergent curve.

Conclusion
In this paper, a new isogeometric topology optimization method is proposed for the integration of design, analysis and optimization framework. The main idea of this approach is to use the T-spline basis functions, which are employed for geometry description and interpolation of unknown field variables, as well as represent the density distribution. The isogeometric topology optimization formulation is developed for the minimization of the structural mean compliance. Several numerical examples demonstrate the feasibility efficiency and flexibility of the proposed method. Satisfactory results can be obtained from arbitrarily complex design domains without additional stabilization and filtering techniques. Moreover, the post-processing is simple and more data-adaptable compared with original density-based methods. It is proved that the topological flexibility of T-spline provides a solid basis for integrating design and analysis through IGA, which also gives topology optimization the probability to be However, the geometric modeling technique of T-spline is not mature enough. Thus the T-spline models with the industry-level complex design are very scarce. In the future, we will extend this approach to shell and 3D structures, as well as large scale industrial structures with the help of new development in geometry modeling technology and the intensive research in the sphere of isogeometric analysis. Additionally, the local refinement feature enables T-splines based method to support a non-uniform IGA mesh. It is of great value to fully study the effect of the IGA mesh distribution on the optimized designs, particularly when solving the stress-related problems [75]. Applications of T-splines in other types of ITO, such as LSM and MMC, also deserve further study.
Funding Statement: This research is supported by the Natural Science Foundation of China (Project Nos. 61972011 and 61572056).

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study. Under the framework of the T-splines based isogeometric topology optimization, the formulation of infill optimization can be given as: Find : q i ði ¼ 1; 2; . . . ; nÞ Min : cðvÞ ¼ F T U S:t: : where α is the local volume fraction and " q e measures the local material distribution around element e, which can be computed by N e is the set of surrounding elements whose centroid is closer than a given local influence radius R e to the centroid of element e, which can be stated as where x j ¼ S j ð0:5; 0:5Þ is the centroid of the j-th element.
As we can see, in Eq. (28), the local volume constraint " q e ⩽ a; 8e may bring several millions of constraints, which is difficult to solve. Thus, this constraint is firstly transformed to maxð" q e Þ ⩽ a. Then the p-norm function is used to approximate the max function as max 8e " q e ð Þ % k" qk p ¼ Finally the local volume constraint " g v ð Þ is stated as q p e ðvÞ " # 1 p À a ⩽ 0: The local volume constraint can be derived as follows @ " q e @v j ¼ @v j @q i ¼ R j i :