On the Indispensability of Isogeometric Analysis in Topology Optimization for Smooth or Binary Designs

: Recently, isogeometric analysis (IGA), which uniﬁes the computer-aided design (CAD) model and the computer-aided engineering (CAE) model, has been adopted to develop the isogeometric topology optimization (ITO) framework. However, a critical study on the indispensability of IGA in topology optimization to take the place of the conventional ﬁnite element method (FEM) is still lacking. In the current work, two important problems are extensively discussed: (1) The lower numerical precision of the FEM resulting from the disuniﬁcation between the CAD and CAE models damages the effectiveness of the topology optimization, which suggests the indispensability of IGA in the replacement of the FEM in optimization; (2) a material penalization model is required to ensure the generation of a full loading-transmission path during optimization in classic density-based methods, which causes a greater overestimation of structural stiffness and also suggests the necessity of an ersatz material model. The current paper describes a promising ITO method with point-wise design to gain smooth or binary symmetrical topologies, for which an extended density distribution function (DDF) was constructed to describe the structural topology. Two benchmarks of the stiffness-maximization problem and compliant mechanism are studied in the context of the above issues. Finally, several topologically optimized designs with symmetry are obtained using the ITO method. approach the 0–1 density distribution, the degree of gray scale is equal to 0.096%, and the corresponding numerical error is equal to − 0.0028%. In conclusion, the developed ITO method with the extended DDF and ersatz material model considering the topology information has a superior ability to solve problems.


Introduction
For the last three decades, topology optimization has been regarded as an effective and efficient numerical tool to find the optimal material distribution with the expected structural performance in fixed-design domains, subject to several prescribed constraints [1]. Since the homogenization method was proposed by Bendsøe and Kikuchi in their seminal work [2], several different topology optimization methods have been proposed in recent years, such as SIMP (solid isotropic material penalization) [3,4], ESO (evolutionary structural optimization) [5,6], LSM (level-set method) [7][8][9], and MMC/Vs (moving morphable components/voids [10,11]). According to the topology description model, these methods can be mainly classified into two branches, namely material description models (MDMs) and boundary description models (BDMs). In BDMs, a higher-dimensional function is constructed to represent the structural topology, where the iso-surface/contour of the function is applied to present structural boundaries, as in the LSM and MMC/Vs. In the In actual, the NURBS basis functions utilized in the construction of the CAD (computeraided design) model can be also applied to construct the numerical response space in the CAE (computer-aided analysis) model, so that an integrated model can be developed in CAD and CAE to eliminate the numerical deficiencies in the traditional FEM (finite element method). This is the core of IGA (isogeometric analysis) [40,41]. Hence, the study of the utilization of IGA in topology optimization to eliminate the numerical issues of FEM (such as an approximant of structural geometry, the lower-order continuity of numerical responses, and poor numerical precision), has attracted many researchers [42][43][44][45]. Previous works studied the application of IGA into topology optimization methods, such as SIMP [35][36][37][38][39], LSM [46][47][48], and MMC/Vs [49][50][51], to study its effectiveness. A brief review of the combinations of IGA and topology optimization can be found in [52]. The aforementioned research works study the effectiveness of the application of IGA in topology optimization to solve several design problems. Although IGA presents several benefits compared with the conventional FEM in numerical analysis and also displays effectiveness in topology optimization, the indispensability of IGA in optimization to replace conventional FEM is still poorly understood.
Hence, the current work first aims to propose a benchmark ITO method of using MDMs, in which an extended DDF is constructed by a NURBS parameterization with a smoothing mechanism to control the smoothness degree and a threshold projection to govern the feature of binary property. The NURBS-based IGA then replaces the conventional FEM to solve structural responses with greater numerical precision; the Gauss quadrature method, considering the exact structural geometry information, is applied to evaluate the stiffness matrices. The material penalization model and the ersatz material model are both studied to show the corresponding features in the point-wise design with detailed topology information, rather than the previous element-wise design, featuring the linear interpolation of the density in finite elements. Finally, several numerical examples with stiffness-maximization and compliant mechanism designs are tested to show the effectiveness of the ITO and the indispensability of IGA in the optimization.

Smooth or 0-1 Designs in Topology Optimization
As shown in Figure 1, Ω is the structural design domain, and in Ω = Ω s ∪ Ω υ , Ω s denotes the structural design domain with the solid and Ω υ indicates the structural design domain with the void. Γ is the structural boundary, namely Γ = Ω s ∩ Ω υ . It is known that the main intention behind topology optimization is to find an appropriate material distribution Ω s with the optimized structural performance, and a classic topology optimization formulation with the volume constraint is given as: where J is an optimized objective function, such as the structural compliance. V max denotes the maximal value of the material consumption. Φ denotes a function that satisfies Φ = 1, i f Φ ∈ Ω s ; Φ = 0, i f Φ ∈ Ω υ to describe the material layout. In the classical SIMP method, Φ is transformed into a series of discrete values defined at finite elements, namely Φ e (e = 1, 2, . . . , N e ), and the initial discrete values of Φ are also relaxed into a continuous variation, namely {0, 1} → [0, 1] . It should be noted that the minimum of Φ e should be equal to a constant that is not equal to 0 to avoid numerical singularity. The corresponding mathematical model can be also written as in a discrete form:

A Brief Description
The classic SIMP method [3,4] has been viewed as a powerful and promising tool for topology optimization in recent years; it has been applied to many engineering optimization problems, a detailed review of which can be found in [53]. As has already been pointed out in many works, the critical feature of the SIMP method is that an isotropic material interpolation model for elastic property with respect to artificial densities can effectively convert the initial discrete design variables into a continuous version. In the latter, the three-field of design variables, topology variables, and physical variables was innovatively studied [53], and, consequently, the three-field SIMP method with the FEM, namely the FEM-based three-field SIMP, has been a benchmark method in recent years. A detailed mathematical formula can be given as: where 0 k denotes the solid-element-stiffness matrix, and e u is the element-displacement vector. U is the global displacement field, and F denotes the imposed load. K is the global stiffness matrix. ( ) V x and 0 V are the material volume fraction and the design domain volume, respectively. x is the vector of the design variables (i.e., artificial element densities e , 1,2, , e x e N =  ). In the FEM-based three-field SIMP, the vector of the design variables in the detailed implementations is involved in the three fields, including the design variables, topology variables and physical variables, details of which can be found in [54] with an 88-line MATLAB code to implement the above formulation, which can be called by a MATALB line top88(nelx, nely, volfrac, penal, rmin, ft) with six input parameters.

Several Critical Problems
As shown in Figure 2, three benchmarks, including cantilever beam, Michell-type structure, and MBB beam, were first optimized using the classic 88 MATLAB SIMP code [54]. The related parameters (nelx, nely, volfrac, penal, rmin, and ft) are listed below the

FEM-Based Three-Field SIMP Method 2.2.1. A Brief Description
The classic SIMP method [3,4] has been viewed as a powerful and promising tool for topology optimization in recent years; it has been applied to many engineering optimization problems, a detailed review of which can be found in [53]. As has already been pointed out in many works, the critical feature of the SIMP method is that an isotropic material interpolation model for elastic property with respect to artificial densities can effectively convert the initial discrete design variables into a continuous version. In the latter, the three-field of design variables, topology variables, and physical variables was innovatively studied [53], and, consequently, the three-field SIMP method with the FEM, namely the FEM-based three-field SIMP, has been a benchmark method in recent years. A detailed mathematical formula can be given as: where k 0 denotes the solid-element-stiffness matrix, and u e is the element-displacement vector. U is the global displacement field, and F denotes the imposed load. K is the global stiffness matrix. V(x) and V 0 are the material volume fraction and the design domain volume, respectively. x is the vector of the design variables (i.e., artificial element densities x e , e = 1, 2, . . . , N e ). In the FEM-based three-field SIMP, the vector of the design variables in the detailed implementations is involved in the three fields, including the design variables, topology variables and physical variables, details of which can be found in [54] with an 88-line MATLAB code to implement the above formulation, which can be called by a MATALB line top88 (nelx, nely, volfrac, penal, rmin, ft) with six input parameters.

Several Critical Problems
As shown in Figure 2, three benchmarks, including cantilever beam, Michell-type structure, and MBB beam, were first optimized using the classic 88 MATLAB SIMP code [54]. The related parameters (nelx, nely, volfrac, penal, rmin, and ft) are listed below the solutions. The final parameter ft for the filtering was fixed at 2 in all the numerical examples, and it is not displayed in Figure 2. The values nelx and nely denote the total number of the finite elements in two physical directions, respectively, and are also equal to the total number of design variables in two directions. The value of volfrac is the maximum value of the material volume fraction, and penal is the penalty parameter. As we can easily see from Figure 2, when a material penalization model with a penalty parameter equal to 3 is used, the optimized topologies with the full loading-transmission path can be achieved, and a number of intermediate densities are also produced. If the ersatz material model is set with a penalty parameter equal to 1, there is no complete loading-transmission path in the final topologies, and a large number of intermediate densities is produced. This mainly results from the ersatz material model, in which the elastic property and density in each finite element have a linear relationship, and it is not possible to enable the generation of the topology with the force-transmission path. mized solutions are provided in the second row of Figure 3. As can be observed, even if the penalization mechanism is not employed in the topology optimization, the threshold projection still has the ability to push the generation of the topology to some extent, such as the Michell-type structure and the MBB beam. However, according to the optimized results, the operation of the threshold projection as a filter is not sufficient to push the densities towards 0 or 1. The optimized solutions are different from the benchmark designs. Moreover, the optimization is not stable and cannot arrive at the convergent condition with a reasonable topology. Hence, Huang [29,30] defined a floating projection constraint in the optimization using the threshold projection that can completely force the densities towards 0 or 1.  We then removed the intermediate densities in the optimized designs obtained by the material penalization model, which are shown in the first row of Figure 2. In [55], several morphology-based black-and-white filters were applied to develop the three-field SIMP method; the greater effectiveness of its threshold projection has received significant attention [56]. As shown in the first row of Figure 3, the threshold projection is employed in the material penalization model, and the settings of the related parameters are consistent with [56]. We found that the densities in final results were almost equal to 0 or 1, and only an extremely limited number of intermediate densities were present. Meanwhile, if the ersatz material model and threshold projection are used, the corresponding optimized solutions are provided in the second row of Figure 3. As can be observed, even if the penalization mechanism is not employed in the topology optimization, the threshold projection still has the ability to push the generation of the topology to some extent, such as the Michell-type structure and the MBB beam. However, according to the optimized results, the operation of the threshold projection as a filter is not sufficient to push the densities towards 0 or 1. The optimized solutions are different from the benchmark designs. Moreover, the optimization is not stable and cannot arrive at the convergent condition with a reasonable topology. Hence, Huang [29,30] defined a floating projection constraint in the optimization using the threshold projection that can completely force the densities towards 0 or 1.
of the finite elements in two physical directions, respectively, and are also equal to the total number of design variables in two directions. The value of volfrac is the maximum value of the material volume fraction, and penal is the penalty parameter. As we can easily see from Figure 2, when a material penalization model with a penalty parameter equal to 3 is used, the optimized topologies with the full loading-transmission path can be achieved, and a number of intermediate densities are also produced. If the ersatz material model is set with a penalty parameter equal to 1, there is no complete loading-transmission path in the final topologies, and a large number of intermediate densities is produced. This mainly results from the ersatz material model, in which the elastic property and density in each finite element have a linear relationship, and it is not possible to enable the generation of the topology with the force-transmission path.
We then removed the intermediate densities in the optimized designs obtained by the material penalization model, which are shown in the first row of Figure 2. In [55], several morphology-based black-and-white filters were applied to develop the three-field SIMP method; the greater effectiveness of its threshold projection has received significant attention [56]. As shown in the first row of Figure 3, the threshold projection is employed in the material penalization model, and the settings of the related parameters are consistent with [56]. We found that the densities in final results were almost equal to 0 or 1, and only an extremely limited number of intermediate densities were present. Meanwhile, if the ersatz material model and threshold projection are used, the corresponding optimized solutions are provided in the second row of Figure 3. As can be observed, even if the penalization mechanism is not employed in the topology optimization, the threshold projection still has the ability to push the generation of the topology to some extent, such as the Michell-type structure and the MBB beam. However, according to the optimized results, the operation of the threshold projection as a filter is not sufficient to push the densities towards 0 or 1. The optimized solutions are different from the benchmark designs. Moreover, the optimization is not stable and cannot arrive at the convergent condition with a reasonable topology. Hence, Huang [29,30] defined a floating projection constraint in the optimization using the threshold projection that can completely force the densities towards 0 or 1.  According to the above numerical results [55][56][57] and recent studies [27][28][29][30][31], on some important problems, the following conclusions can be drawn.

•
As shown in Figure 4, if the design domain uses a denser mesh, the numerical precision is improved and the optimization with the ersatz material model and threshold projection can find relatively appropriate designs. This reveals that numerical precision has a significant influence on the ability of threshold projections to generate structural topologies. However, increasing the number of finite elements results in a prohibitive computational cost, Although the optimization can find solutions to some extent, the densities dramatically change in the optimization and cannot converge in a stable process. Hence, the use of IGA to unify CAD and CAE in an integrated model with greater numerical precision at the same finite elements can offer a powerful alternative to FEM for solving unknown responses [40,41]. mechanism in the density distribution.

•
It is imperative to construct a more effective DDF to represent the structural topology with smoothness by using NURBS approximation mechanism rather than the Lagrange interpolation. Furthermore, the DDF can combine with the threshold projection to extensively remove the intermediate densities in the transition area, and the optimizer can avoid the discrete characteristic in the topology description using the DDF. Additionally, during the optimization, the DDF with the implicit representation of structural boundaries is applied to describe the structural topology, and the DDF at the Gauss quadrature points is utilized in IGA to evaluate the material elastic properties. Hence, the current ITO can eliminate the considerable dependence on elementary densities in previous works.

An Extended DDF to Implicitly Represent Structural Topologies
In the optimization, the DDF is employed to represent the structural topology, and structural boundaries are immersed into the iso-surface/contour of the DDF. As already discussed in [38], the definition of the DDF should maintain three conditions: (1) Nonnegativity; (2) strict binding within 0-1; and (3) non-negative derivatives. Moreover, the smoothness and continuity of the DDF have a significant effect on the optimized solutions. The development of the smoothness mechanism using the Shepard function and maintaining higher-order continuity using NURBS basis functions are imperative. However, although the optimized designs have smooth structural boundaries and distinct material

•
Previous works [29,30] revealed that the ersatz material model is more suitable for optimization to obtain a smooth design in the density framework and the material penalization model can further benefit the optimization to obtain a pure 0-1 design. Meanwhile, the threshold projection is introduced to define the floating constraint to push the design variables towards 0 or 1. In fact, the material penalization model, with its overestimation of structural compliance, leads to numerical deviations in the structural stiffness of the final topologies, and a unified topology optimization based on the ersatz material model, using the density to produce smooth designs and pure 0-1 designs is more meaningful. Moreover, the floating constraint increases the numerical difficulties of the optimization. It is imperative to define the projection mechanism in the density distribution. • It is imperative to construct a more effective DDF to represent the structural topology with smoothness by using NURBS approximation mechanism rather than the Lagrange interpolation. Furthermore, the DDF can combine with the threshold projection to extensively remove the intermediate densities in the transition area, and the optimizer can avoid the discrete characteristic in the topology description using the DDF. Additionally, during the optimization, the DDF with the implicit representation of structural boundaries is applied to describe the structural topology, and the DDF at the Gauss quadrature points is utilized in IGA to evaluate the material elastic properties. Hence, the current ITO can eliminate the considerable dependence on elementary densities in previous works.

An Extended DDF to Implicitly Represent Structural Topologies
In the optimization, the DDF is employed to represent the structural topology, and structural boundaries are immersed into the iso-surface/contour of the DDF. As already discussed in [38], the definition of the DDF should maintain three conditions: (1) Nonnegativity; (2) strict binding within 0-1; and (3) non-negative derivatives. Moreover, the smoothness and continuity of the DDF have a significant effect on the optimized solutions. The development of the smoothness mechanism using the Shepard function and maintaining higher-order continuity using NURBS basis functions are imperative. However, although the optimized designs have smooth structural boundaries and distinct material interfaces, the transition area from 0 to 1 in the designs has a number of intermediate densities [38,39,52], which enhance the numerical errors due to the material penalization model in the evaluation of material elastic properties. In the current work, we first focus on developing an extended DDF by introducing threshold projection, which can reduce the sities [38,39,52], which enhance the numerical errors due to the material penalization model in the evaluation of material elastic properties. In the current work, we first focus on developing an extended DDF by introducing threshold projection, which can reduce the number of intermediate densities in the optimization. The construction of the extended DDF involves three components: (1) Smoothness; (2) continuity; and (3) binary property. The detailed steps are given as follows, shown in Figure 5:

•
Step 1: Introducing the design variables at the control points, namely the controldesign variables: Step 2: Using the smoothing mechanism to improve the overall smoothness of the control-design variables: where φ  denotes the smoothed control-design variables. ψ is the Shepard function, which is defined by C 4 compactly supported radial basis functions ϕ .

•
Step 4: Using NURBS to construct a density-response surface for the design domain: where Φ is a NURBS density-response-surface with a higher-dimension, whole-design domain, and it can be viewed as a density distribution function representing the material distribution. Compared to the initial DDF in [38,39], the current DDF can be viewed as an extended but powerful function, which can describe not only the smooth structural topology that can be fabricated, but also control the intermediate densities to extensively reduce the numerical errors in the evaluation of the structural performance.

•
Step 1: Introducing the design variables at the control points, namely the controldesign variables: • Step 2: Using the smoothing mechanism to improve the overall smoothness of the control-design variables: where φ denotes the smoothed control-design variables. ψ is the Shepard function, which is defined by C 4 compactly supported radial basis functions ϕ.

•
Step 3: Introducing threshold projection [55][56][57] to reduce the intermediate densities in the transition area:φ where β and η are the related parameters in threshold projection.

•
Step 4: Using NURBS to construct a density-response surface for the design domain: where Φ is a NURBS density-response-surface with a higher-dimension, whole-design domain, and it can be viewed as a density distribution function representing the material distribution. Compared to the initial DDF in [38,39], the current DDF can be viewed as an extended but powerful function, which can describe not only the smooth structural topology that can be fabricated, but also control the intermediate densities to extensively reduce the numerical errors in the evaluation of the structural performance. Furthermore, NURBS is a response surface for the design domain. In the construction of the extended DDF, the densities at the control points work as control coefficients, which means that the DDF can be also viewed as a density-response surface with a higherdimension structural design domain. In the representation of the structural topology, the iso-contour/surface of the DDF is applied to represent the structural boundaries. The DDF with higher values than the iso-contour represents the solids in the design domain, and the lower values than the iso-value in the DDF describe the voids in the structural design domain, shown in Figure 6. The detailed mathematical model is expressed as follows: where Φ iso is the value of the iso-contour/surface of the DDF, which is equal to 0.5 in the latter numerical examples. However, it should be noted that although the immersed representation of the structural boundaries using the DDF inherits from the level-set implicit description model, the optimization of the DDF has an intrinsic difference compared to LSM. In the current ITO using the DDF, the main intention is to optimize the DDF by advancing the control-design variables, rather than deriving the moving and merging of the iso-contour/surface of the DDF.
of the extended DDF, the densities at the control points work as control coefficients, which means that the DDF can be also viewed as a density-response surface with a higher-dimension structural design domain. In the representation of the structural topology, the iso-contour/surface of the DDF is applied to represent the structural boundaries. The DDF with higher values than the iso-contour represents the solids in the design domain, and the lower values than the iso-value in the DDF describe the voids in the structural design domain, shown in Figure 6. The detailed mathematical model is expressed as follows: where iso Φ is the value of the iso-contour/surface of the DDF, which is equal to 0.5 in the latter numerical examples. However, it should be noted that although the immersed representation of the structural boundaries using the DDF inherits from the level-set implicit description model, the optimization of the DDF has an intrinsic difference compared to LSM. In the current ITO using the DDF, the main intention is to optimize the DDF by advancing the control-design variables, rather than deriving the moving and merging of the iso-contour/surface of the DDF.

NURBS-Based IGA
In IGA, the same NURBS basis functions used in the construction of the structural geometry are still applied to develop the numerical response space, which can ensure the integration of the CAD and CAE models. A detailed derivation of linearly elastic problems can be found in [38]. In the case of the Galerkin IGA formulation, the system stiffness matrix is achieved by assembling all the IGA element-stiffness matrices. In the physical space, the calculation of the IGA element-stiffness matrix can be stated as: where e K is the IGA element stiffness matrix, e Ω is the IGA element domain, and e B indicates the strain-displacement matrix. In the calculation, the iso-parametric formulation is employed, in which two mappings should be defined, namely the mapping of : e e Ω → Ω X from the parametric space to the physical space and an affine mapping of : from the bi-unit parent element to the parametric element. The integration in Equation (9) needs to be pulled back first into the parametric element and then into the bi-unit parent element. It is involved the inverses of the mappings X and Y , and the detailed form is given by:

NURBS-Based IGA
In IGA, the same NURBS basis functions used in the construction of the structural geometry are still applied to develop the numerical response space, which can ensure the integration of the CAD and CAE models. A detailed derivation of linearly elastic problems can be found in [38]. In the case of the Galerkin IGA formulation, the system stiffness matrix is achieved by assembling all the IGA element-stiffness matrices. In the physical space, the calculation of the IGA element-stiffness matrix can be stated as: where K e is the IGA element stiffness matrix, Ω e is the IGA element domain, and B e indicates the strain-displacement matrix. In the calculation, the iso-parametric formulation is employed, in which two mappings should be defined, namely the mapping of X :Ω e → Ω e from the parametric space to the physical space and an affine mapping of Y : Ω e →Ω e from the bi-unit parent element to the parametric element. The integration in Equation (9) needs to be pulled back first into the parametric element and then into the bi-unit parent element. It is involved the inverses of the mappings X and Y, and the detailed form is given by: where J 1 and J 2 are the Jacobi matrices of the two mappings, respectively. The Gauss quadrature method is adopted to calculate Equation (10), and the detailed computation is expressed by: where w i and w j are the corresponding quadrature weights. D 0 is material constitutive elastic tensor matrix of solid. Φ ξ i , η j denotes the density at the Gauss quadrature point ξ i , η j . In Figure 7, we present the comparison of the FEM and IGA in terms of the calculation of the element-stiffness matrix. As can be observed, in the FEM, a parametric element without the consideration of the structural topology information is utilized in the calculation of all the finite element-stiffness matrices. In the FEM, if the element densities are equal to 1, all the finite-element matrices are the same. However, in IGA, each elementstiffness matrix has its own parametric element considering the topology information. The stiffness matrices of all the solid elements are different. Hence, the calculation of the IGA element-stiffness matrices with structural topology information is performed pointto-point, which can offer higher numerical precision compared to the FEM. Moreover, it is noted that the current ersatz material model is considered for each point with the topological information, rather than the conventional element-wise design, which features a linear approximation, using the density.
where 1 J and 2 J are the Jacobi matrices of the two mappings, respectively. The Gauss quadrature method is adopted to calculate Equation (10), and the detailed computation is expressed by: where i w and j w are the corresponding quadrature weights. 0 D is material constitutive elastic tensor matrix of solid. ( ) denotes the density at the Gauss quadrature Figure 7, we present the comparison of the FEM and IGA in terms of the calculation of the element-stiffness matrix. As can be observed, in the FEM, a parametric element without the consideration of the structural topology information is utilized in the calculation of all the finite element-stiffness matrices. In the FEM, if the element densities are equal to 1, all the finite-element matrices are the same. However, in IGA, each elementstiffness matrix has its own parametric element considering the topology information.
The stiffness matrices of all the solid elements are different. Hence, the calculation of the IGA element-stiffness matrices with structural topology information is performed pointto-point, which can offer higher numerical precision compared to the FEM. Moreover, it is noted that the current ersatz material model is considered for each point with the topological information, rather than the conventional element-wise design, which features a linear approximation, using the density.

Comparisons between the ITO and the FEM-Based Three-Field SIMP
As shown in Figure 8, two flowcharts of the ITO method and the FEM-based threefield SIMP are provided. The distinct differences mainly contain the following points: (1) The extended DDF in the ITO method offers a promising approach to topology representation using a smooth and continuous function, rather than a series of discrete-element densities in the three-field SIMP method; (2) the same NURBS basis functions adopted in the extended DDF and numerical analysis can effectively ensure their integration in the topology optimization. However, the linear function in the construction of the three-field

Comparisons between the ITO and the FEM-Based Three-Field SIMP
As shown in Figure 8, two flowcharts of the ITO method and the FEM-based three-field SIMP are provided. The distinct differences mainly contain the following points: (1) The extended DDF in the ITO method offers a promising approach to topology representation using a smooth and continuous function, rather than a series of discrete-element densities in the three-field SIMP method; (2) the same NURBS basis functions adopted in the extended DDF and numerical analysis can effectively ensure their integration in the topology optimization. However, the linear function in the construction of the three-field element densities and the shape function in the numerical analysis using FEM are completely different, which can result in numerical deviations; (3) in the numerical analysis, the geometry information of the design domain and the topology information of the DDF can be exactly considered in the IGA, whereas only the topology information stemming from the discrete element densities is utilized in the FEM. The difference is discussed in the last paragraph of Section 3.2, as shown in Figure 7. The development of the current ITO method is based on the FEM-based three-field SIMP method. For example, the filtering scheme and a threshold projection in the FEM-based three-field SIMP method are inherited in the construction of the DDF, where the only difference is the NURBS parameterization to ensure the features of the continuity. Overall, the application of IGA in the optimization can improve the numerical precision. However, the effectiveness of the higher numerical precision under the same finite elements in IGA and the indispensability of IGA are still not applied in topology optimization.
pletely different, which can result in numerical deviations; (3) in the numerical analysis, the geometry information of the design domain and the topology information of the DDF can be exactly considered in the IGA, whereas only the topology information stemming from the discrete element densities is utilized in the FEM. The difference is discussed in the last paragraph of Section 3.2, as shown in Figure 7. The development of the current ITO method is based on the FEM-based three-field SIMP method. For example, the filtering scheme and a threshold projection in the FEM-based three-field SIMP method are inherited in the construction of the DDF, where the only difference is the NURBS parameterization to ensure the features of the continuity. Overall, the application of IGA in the optimization can improve the numerical precision. However, the effectiveness of the higher numerical precision under the same finite elements in IGA and the indispensability of IGA are still not applied in topology optimization.

The Problems of Compliance-Minimization and the Compliant Mechanism
The optimization problems of structural compliance and the compliant mechanism are discussed to show the effectiveness of ITO with the extended DDF and ersatz material model.

Compliance-Minimization Formulation
As already discussed in Section 3, the control densities work as design variables, and the optimizer focuses on finding the optimized DDF with the expected structural performance. As far as the structural stiffness problem is concerned, the objective function is structural compliance, and the detailed formula is given as:

The Problems of Compliance-Minimization and the Compliant Mechanism
The optimization problems of structural compliance and the compliant mechanism are discussed to show the effectiveness of ITO with the extended DDF and ersatz material model.

Compliance-Minimization Formulation
As already discussed in Section 3, the control densities work as design variables, and the optimizer focuses on finding the optimized DDF with the expected structural performance. As far as the structural stiffness problem is concerned, the objective function is structural compliance, and the detailed formula is given as: where φ i,j denotes the (i, j) th control-design variable. J is the objective function, which is defined by the mean structural compliance. G corresponds to the volume constraint, where V max denotes the maximum material consumption. Φ is the DDF to represent the structural topology. a denotes the bilinear energy function, and l is the linear load function. u is the displacement field in structural design domain, and g is the prescribed displacement vector at the Dirichlet boundary Γ D . δu is the virtual displacement field belonging to the Sobolev space H 1 (Ω). The variational weak form of the equilibrium-state equation can be established, and the detailed form is given as: where f is the body force and h is the boundary traction on the Neumann boundary Γ N . Meanwhile, the sensitivities of the objective function and constraint functions with respect to the control-design variables are required in the gradient-based optimization algorithm to derive the advancement of the structural topology. The details can be found in [38,39]. The first-order derivative of structural compliance with respect to the controldesign variables can be expressed as: Similarly, the derivative of the volume constraint with respect to the control-design variables is given by:

Compliant Mechanism Design Formulation
As far as the optimization of the compliant mechanism is concerned, the objective function should correspond to the maximum output displacement, and the material consumption works as the constraint function. The final form of the mathematical model can be stated as: Find : φ i,j (i = 1, 2, . . . , n; j = 1, 2, . . . , m) where u out is the output displacement. L is a column vector; the value of L at the output position is equal to 1, and the other values are equal to 0. The first-order derivative of the output displacement with respect to the control-design variables can be given as follows. First, a general form of the derivative of the objective function with respect to the control-design variables is given as: According to the elastic equilibrium equation, a new form of Equation (17) is given as: Hence, we gain a detailed form that can be expressed by: The derivative of the volume constraint with the control-design variables remains consistent with Equation (15). A known gradient-based mathematical optimization algorithm, namely the OC (optimality criteria) method, is applied to solve the optimization problems in the current work.

Numerical Examples
In this section, several 2D and 3D numerical examples are performed to show the effectiveness of the ITO method. Three benchmark examples, namely cantilever beam, Michell-type structure, and MBB beam, are considered to maximize the structural stiffness. First, the importance of the threshold projection in the extended DDF is sufficiently discussed to present the indispensability of the extended DDF compared to the initial DDF in [38,39]. Second, the effectiveness and indispensability of the IGA in replacing the conventional FEM in topology optimization are extensively discussed. Third, the problem of stiffness overestimation using the material penalization model is studied, and the effectiveness of the use of the ITO method with the ersatz material model to achieve smooth or binary topologies is also addressed. Next, three 3D numerical examples are also tested using the ITO method with the ersatz material model. Finally, the example of a classic force-inventor design is also studied using the ITO method with the ersatz material model to achieve smooth or binary designs. In all the examples, the Young's moduli of the solid and void materials are defined as 1 and 1 × 10 −9 to avoid the numerical singularity, respectively. The Poisson's ratio is set as 0.3. The imposed point load magnitude is defined as 1, except for the special definition. Gauss quadrature points with 2 × 2 (2D) or 2 × 2 × 2 (3D) are chosen in each IGA element. In the threshold projection, the parameter η is set as 0.5, β increases within the optimization, and the increasing criterion is listed in Algorithm 1. It is noted that only linearly elastic materials are considered. The initial values of the control-design variables in all the examples are set at 1. The radius in the Shepard function r m is equal to 2 in the numerical examples. The convergent condition is that the L ∞ norm of the difference of control-design variables between two consecutive iterations is less than 1% within the maximum 400 steps. In the case of three benchmarks, the corresponding boundary and load conditions are shown in Figure 9, where the structural sizes are denoted by L and H.

The Effectiveness of the Extended DDF
In this subsection, the key intention is to study the effectiveness of the extended DDF and show the necessity of the initial DDF. The material penalization model is used in this subsection, and the penalty parameter is set as 3. The detailed design parameters of the three benchmarks are defined in Table 1, including the structural sizes, the degrees of the NURBS basis functions in parametric directions, the number of control points or controldesign variables, the number of IGA elements, and the maximum material consumption.
As shown in Figure 10, the optimized cantilever beam solutions using the ITO method with the initial DDF and the extended DDF are provided, including the optimized layouts of the control-design variables, the DDF, and the corresponding topologies. The optimized designs of the Michell structure and MBB beam are presented in Figures 11 and

The Effectiveness of the Extended DDF
In this subsection, the key intention is to study the effectiveness of the extended DDF and show the necessity of the initial DDF. The material penalization model is used in this subsection, and the penalty parameter is set as 3. The detailed design parameters of the three benchmarks are defined in Table 1, including the structural sizes, the degrees of the NURBS basis functions in parametric directions, the number of control points or controldesign variables, the number of IGA elements, and the maximum material consumption. As shown in Figure 10, the optimized cantilever beam solutions using the ITO method with the initial DDF and the extended DDF are provided, including the optimized layouts of the control-design variables, the DDF, and the corresponding topologies. The optimized designs of the Michell structure and MBB beam are presented in Figures 11 and 12, respectively. As we can easily find, the optimized values of the control-design variables in    Figure 10. The optimized designs of cantilever beam: (a1-c1) The ITO method using the extended DDF; (a2-c2) the ITO method using the initial DDF. Figure 11. The optimized designs of the Michell structure: (a1-c1) The ITO method using the extended DDF; (a2-c2) the ITO method using the initial DDF. Figure 11. The optimized designs of the Michell structure: (a1-c1) The ITO method using the extended DDF; (a2-c2) the ITO method using the initial DDF. Figure 11. The optimized designs of the Michell structure: (a1-c1) The ITO method using the extended DDF; (a2-c2) the ITO method using the initial DDF. Figure 12. The optimized designs of the MBB beam: (a1-c1) The ITO method using the extended DDF; (a2-c2) the ITO method using the initial DDF. Figure 12. The optimized designs of the MBB beam: (a1-c1) The ITO method using the extended DDF; (a2-c2) the ITO method using the initial DDF.
4, x FOR PEER REVIEW 15 of 28 Figure 13. The data layout of the density values of the control-design variables in the optimized designs presented: (a) the corresponding layout in Figure 10(a1); (b) the corresponding layout in Figure 11(a1); (c) the corresponding layout in Figure 12(a1); (d) the corresponding layout in Figure  10(a2); (e) the corresponding layout in Figure 11(a2); (f) the corresponding layout in Figure 12(a2).

The Indispensability of the IGA in Replacing the FEM in Topology Optimization
In this subsection, the key intention is to discuss the indispensability of the IGA in replacing the FEM in topology optimization. First, the FEM-based three-field SIMP method is applied to optimize the three examples above with the same design parameters, namely the number of finite elements, the maximum volume fraction, and the same increasing criterion of the parameters in the threshold projection. In Figure 14, the optimized topologies of the three examples are provided. The topologies in the first row are optimized by the FEM-based three-field SIMP method with the material penalization model, and the second row corresponds to the designs optimized by the FEM-based three-field SIMP method with the ersatz material model. The related parameters (nelx, nely, volfrac, penal, and rmin) are defined below the optimized topologies. We can easily infer that the FEM-based three-field SIMP method with the ersatz material model cannot achieve a reasonable topology with a full loading-transmission path, and the optimization cannot arrive at the convergent condition.  Figure 11(a1); (c) the corresponding layout in Figure 12(a1); (d) the corresponding layout in Figure 10(a2); (e) the corresponding layout in Figure 11(a2); (f) the corresponding layout in Figure 12(a2).

The Indispensability of the IGA in Replacing the FEM in Topology Optimization
In this subsection, the key intention is to discuss the indispensability of the IGA in replacing the FEM in topology optimization. First, the FEM-based three-field SIMP method is applied to optimize the three examples above with the same design parameters, namely the number of finite elements, the maximum volume fraction, and the same increasing criterion of the parameters in the threshold projection. In Figure 14, the optimized topologies of the three examples are provided. The topologies in the first row are optimized by the FEM-based three-field SIMP method with the material penalization model, and the second row corresponds to the designs optimized by the FEM-based three-field SIMP method with the ersatz material model. The related parameters (nelx, nely, volfrac, penal, and rmin) are defined below the optimized topologies. We can easily infer that the FEMbased three-field SIMP method with the ersatz material model cannot achieve a reasonable topology with a full loading-transmission path, and the optimization cannot arrive at the convergent condition. and the second row corresponds to the designs optimized by the FEM-based three-field SIMP method with the ersatz material model. The related parameters (nelx, nely, volfrac, penal, and rmin) are defined below the optimized topologies. We can easily infer that the FEM-based three-field SIMP method with the ersatz material model cannot achieve a reasonable topology with a full loading-transmission path, and the optimization cannot arrive at the convergent condition. Secondly, the proposed ITO method is employed to optimize the three benchmark examples; the related design parameters are also listed in Table 1. A critical difference is  Table 1. A critical difference is that the ersatz material model is used in the current subsection. As already discussed in Section 3, there is no penalization mechanism in the ersatz material model. As shown in Figure 15, three topologies optimized by the ITO method with the ersatz material model are presented, which are similar to the optimized designs shown in Figures 10-12. Meanwhile, the convergent histories of the structural compliance of the three benchmark examples are shown in Figure 16, plotted in red. It can be easily seen that the iterative curves of the three benchmarks are all characterized by a clear and smooth variation during the optimization. Additionally, some slight fluctuations emerge in the convergent process, such as the iteration 29, 70 in the iterative curve. This results from the increase in the parameter Beta in the threshold projection during the optimization, which intends to push the control-design variables toward 0 or 1. Hence, we may conclude that the extended DDF with the threshold projection in the current ITO method tends to gradually approach the binary layout with a reasonable topology when the parameter Beta of the threshold projection increases during the optimization. that the ersatz material model is used in the current subsection. As already discussed in Section 3, there is no penalization mechanism in the ersatz material model. As shown in Figure 15, three topologies optimized by the ITO method with the ersatz material model are presented, which are similar to the optimized designs shown in Figures 10-12. Meanwhile, the convergent histories of the structural compliance of the three benchmark examples are shown in Figure 16, plotted in red. It can be easily seen that the iterative curves of the three benchmarks are all characterized by a clear and smooth variation during the optimization. Additionally, some slight fluctuations emerge in the convergent process, such as the iteration 29, 70 in the iterative curve. This results from the increase in the parameter Beta in the threshold projection during the optimization, which intends to push the control-design variables toward 0 or 1. Hence, we may conclude that the extended DDF with the threshold projection in the current ITO method tends to gradually approach the binary layout with a reasonable topology when the parameter Beta of the threshold projection increases during the optimization. Additionally, several intermediate cantilever-beam topologies are presented in Figure 17. We can easily see that the earlier iterations with the increase in the parameter Beta facilitate the generation of the key structural members in the design domain, until a full loading-transmission path is generated in the final topology. In the latter iterations, the increase in the parameter Beta aims to push all the control-design variables towards 0 or 1 and form a nearly pure 0-1 distribution of the extended DDF. However, the FEM-based three-filed SIMP method with the ersatz material model has no ability to seek a reasonable topology in a stable process, as presented in Figure 14(a2-c2). With the increase in the finite elements in the FEM to improve the numerical precision, it can enhance the stability of the optimization to find several ineffective designs, as shown in Figure 3. Hence, the higher numerical precision of IGA in replacing the FEM can maintain the effectiveness of the ITO method with the ersatz material model to find an appropriate design. Hence, it is indispensable to develop the framework of the ITO to seek the optimal material distribution in the design domain.   Figure 17. We can easily see that the earlier iterations with the increase in the parameter Beta facilitate the generation of the key structural members in the design domain, until a full loadingtransmission path is generated in the final topology. In the latter iterations, the increase in the parameter Beta aims to push all the control-design variables towards 0 or 1 and form a nearly pure 0-1 distribution of the extended DDF. However, the FEM-based three-filed SIMP method with the ersatz material model has no ability to seek a reasonable topology in a stable process, as presented in Figure 14(a2-c2). With the increase in the finite elements in the FEM to improve the numerical precision, it can enhance the stability of the optimization to find several ineffective designs, as shown in Figure 3. Hence, the higher numerical precision of IGA in replacing the FEM can maintain the effectiveness of the ITO method with the ersatz material model to find an appropriate design. Hence, it is indispensable to develop the framework of the ITO to seek the optimal material distribution in the design domain.
the ITO method with the ersatz material model to find an appropriate design. Hence, it is indispensable to develop the framework of the ITO to seek the optimal material distribution in the design domain.

The Effectiveness of the Ersatz Material Model Compared to the Material Penalization Model in the ITO
As already described in [29,30], the use of the material penalization model to evaluate elastic properties can lead to the critical problem of the overestimation of the structural stiffness in the optimization. In Section 5.1, the ITO method with the material penalization model is employed to optimize three benchmark examples, and the ITO method with the ersatz material model is used in Section 5.1 to discuss the indispensability of IGA. Based on several numerical results in [29,30], the problem of stiffness overestimation always emerges with the intermediate densities. The development of the extended DDF aims to reduce a nearly binary topology in the ITO method. This means that there is no difference between the material penalization model and the ersatz material model when a binary design is formed in the optimization. However, while a binary design can completely reduce the intermediate values of the density, it also introduces several difficulties later in the manufacturing process. Hence, it is compulsory to retain a limited number of interme-

The Effectiveness of the Ersatz Material Model Compared to the Material Penalization Model in the ITO
As already described in [29,30], the use of the material penalization model to evaluate elastic properties can lead to the critical problem of the overestimation of the structural stiffness in the optimization. In Section 5.1, the ITO method with the material penalization model is employed to optimize three benchmark examples, and the ITO method with the ersatz material model is used in Section 5.1 to discuss the indispensability of IGA. Based on several numerical results in [29,30], the problem of stiffness overestimation always emerges with the intermediate densities. The development of the extended DDF aims to reduce a nearly binary topology in the ITO method. This means that there is no difference between the material penalization model and the ersatz material model when a binary design is formed in the optimization. However, while a binary design can completely reduce the intermediate values of the density, it also introduces several difficulties later in the manufacturing process. Hence, it is compulsory to retain a limited number of intermediate densities to develop a small transition area that can construct smooth boundaries from solids to voids using the iso-contour of the extended DDF. The index used to present the gray scale should be defined to rationally choose a criterion that can determine the smooth or binary designs in the optimization. In a similar manner to [55,56], the detailed mathematical equation is stated as: where ϑ is the indicator parameter of the degree of the gray scale in the optimized design. As already described in Section 4, the developed DDF is applied to not only represent structural topology, but also calculate the related material elastic properties. However, as described in [29][30][31]38], the structural compliance of the final topology represented by the iso-value of the DDF is not equal to the optimized objective function calculated by the DDF; this difference is due to the direct step process of the density in the representation of solids and voids. Hence, we should re-analysis the optimized topology and calculate the numerical error. A symbol J 0_1 is applied to denote the re-calculated structural compliance according to the optimized topology, and the detailed mathematical model of the numerical error can be expressed by: where ∆ denotes the numerical error. In Figures 15 and 18, the optimized topologies of the three benchmark examples are provided. The former presents the topologies optimized by the ITO method with the ersatz material model, and the latter gives the results obtained by the ITO method with the material penalization model. The corresponding designs have slight differences. Meanwhile, the related objective functions in the three benchmarks, namely the optimized objective function calculated by the DDF and the re-calculated structural compliance using the optimized topology, their numerical error, and the total iterative step, are provided in Table 2. As can easily be seen, the optimized objective functions of three benchmarks using the ITO with the material penalization model are all larger than the values using the ersatz material model. Meanwhile, the numerical results obtained using the material penalization model have larger objective functions of the extended DDF than the re-calculated structural compliance, with a deviation of almost 10%, which suggests the overestimation of the material penalization model in the evaluation of structural compliance.
where Δ denotes the numerical error. In Figures 15 and 18, the optimized topologies of the three benchmark examples are provided. The former presents the topologies optimized by the ITO method with the ersatz material model, and the latter gives the results obtained by the ITO method with the material penalization model. The corresponding designs have slight differences. Meanwhile, the related objective functions in the three benchmarks, namely the optimized objective function calculated by the DDF and the re-calculated structural compliance using the optimized topology, their numerical error, and the total iterative step, are provided in Table 2. As can easily be seen, the optimized objective functions of three benchmarks using the ITO with the material penalization model are all larger than the values using the ersatz material model. Meanwhile, the numerical results obtained using the material penalization model have larger objective functions of the extended DDF than the re-calculated structural compliance, with a deviation of almost 10%, which suggests the overestimation of the material penalization model in the evaluation of structural compliance.    Figure 18a are shown in Figure 19. It can be seen that the earlier iterations quickly push the control-design variables toward 0 or 1 due to the material penalization model. After the generation of a full loading-transmission path in the design domain, the initial appropriate topology is formed. In the latter iterations, the role of the threshold projection is to work as a step function to gradually reduce the degree of gray scale, as clearly shown by the convergent curves of the degree of the gray scale in Figure 20. Hence, in the ITO method with the material penalization model, the penalization mechanism plays a key role in generating the topology, whereas the threshold projection in the extended DDF affords this function in the ITO method with the ersatz material model.  Figure 21. As can be easily observed, if the parameter Beta in the threshold projection increases, the degree of gray scale gradually decreases and reaches 0.092% (almost equal to 0). Meanwhile, the re-analysis of structural compliance is performed after the generation of the complete loading-transmission path in the optimized topology. As can be observed, the numerical error also increases until the convergent condition is reached, and the final value is equal to −1.8%. This reveals one important feature in the optimization using the ITO method with the ersatz material model: the structural compliance calculated from the optimized topology is lower than the objective function evaluated from the DDF in the optimization.   14, x FOR PEER REVIEW 19 of 28 control-design variables. Even if the values of control-design variables approach 0 or 1, some intermediate densities are still introduced at the Gauss quadrature points because the values of the NURBS basis functions range from 0 to 1. Similarly, the degree of gray scale and numerical error in the designs of the three benchmarks optimized by the ITO method with the ersatz material model are shown in Figure 21. As can be easily observed, if the parameter Beta in the threshold projection increases, the degree of gray scale gradually decreases and reaches 0.092% (almost equal to 0). Meanwhile, the re-analysis of structural compliance is performed after the generation of the complete loading-transmission path in the optimized topology. As can be observed, the numerical error also increases until the convergent condition is reached, and the final value is equal to −1.8%. This reveals one important feature in the optimization using the ITO method with the ersatz material model: the structural compliance calculated from the optimized topology is lower than the objective function evaluated from the DDF in the optimization.   In Figure 20, the iterative histories of the degree of gray scale and the numerical error in the ITO with the material penalization model for the three benchmark examples are provided. With the increase in the parameter Beta in the threshold projection, the numerical error gradually becomes reduces, and finally reaches 8.46% of the cantilever beam, even if the degree of gray scale is extremely small, equal to 0.073%. This difference in numerical error mainly results from the fact that the objective function is calculated by the DDF at the Gauss quadrature points. As defined in Section 3, the values of the extended DDF at the Gauss quadrature points are calculated by the NRUBS basis functions and control-design variables. Even if the values of control-design variables approach 0 or 1, some intermediate densities are still introduced at the Gauss quadrature points because the values of the NURBS basis functions range from 0 to 1. Similarly, the degree of gray scale and numerical error in the designs of the three benchmarks optimized by the ITO method with the ersatz material model are shown in Figure 21. As can be easily observed, if the parameter Beta in the threshold projection increases, the degree of gray scale gradually decreases and reaches 0.092% (almost equal to 0). Meanwhile, the re-analysis of structural compliance is performed after the generation of the complete loading-transmission path in the optimized topology. As can be observed, the numerical error also increases until the convergent condition is reached, and the final value is equal to −1.8%. This reveals one important feature in the optimization using the ITO method with the ersatz material model: the structural compliance calculated from the optimized topology is lower than the objective function evaluated from the DDF in the optimization. Additionally, the other two cantilever-beam numerical designs are shown in Figure  22, where the 30th design in Figure 22a is consistent with the final topology in [38]. As can be seen, the numerical error is significantly larger, equal to 28.9%. This introduces a large numerical deviation. Even if the degree of gray scale decreases to 5.8%, as in Figure 22b, the corresponding numerical error, equal to 15.2%, is still much larger than the value of the error in Figure 23b, which is equal to −0.92%. In Figure 23, three of the numerical designs of the cantilever beam are given. In Figure 23b, the degree of gray scale is equal to 6.8%, the optimized design features smooth structural boundaries, and the final numerical error is extremely small, equal to −0.92%. Hence, an optimized topology with smooth boundaries but an extremely small numerical error can be obtained, which can perfectly meet the requirements of practical engineering using the current ITO method. Hence, it can be concluded that the ITO method with the ersatz material model has superior capabilities. The limited value of the degree of gray scale is 5%, based on the numerical results. When the value is larger than 5%, the optimized topologies have smooth boundaries, and vice versa. Additionally, the other two cantilever-beam numerical designs are shown in Figure 22, where the 30th design in Figure 22a is consistent with the final topology in [38]. As can be seen, the numerical error is significantly larger, equal to 28.9%. This introduces a large numerical deviation. Even if the degree of gray scale decreases to 5.8%, as in Figure 22b, the corresponding numerical error, equal to 15.2%, is still much larger than the value of the error in Figure 23b, which is equal to −0.92%. In Figure 23, three of the numerical designs of the cantilever beam are given. In Figure 23b, the degree of gray scale is equal to 6.8%, the optimized design features smooth structural boundaries, and the final numerical error is extremely small, equal to −0.92%. Hence, an optimized topology with smooth boundaries but an extremely small numerical error can be obtained, which can perfectly meet the requirements of practical engineering using the current ITO method. Hence, it can be concluded that the ITO method with the ersatz material model has superior capabilities. The limited value of the degree of gray scale is 5%, based on the numerical results. When the value is larger than 5%, the optimized topologies have smooth boundaries, and vice versa. boundaries but an extremely small numerical error can be obtained, which can perfectly meet the requirements of practical engineering using the current ITO method. Hence, it can be concluded that the ITO method with the ersatz material model has superior capabilities. The limited value of the degree of gray scale is 5%, based on the numerical results. When the value is larger than 5%, the optimized topologies have smooth boundaries, and vice versa.

Discussion of the Influence of the IGA Mesh
In the current subsection, the key focus is to address the influence of the IGA mesh on the effectiveness of the threshold projection in the ITO method, and also to discuss the effect of the IGA mesh on the numerical error using the immersed representation model of the extended DDF. In this subsection, the ersatz material model is used in the ITO method. The other design parameters remain consistent with the second row of Table 1 and with the change in the IGA mesh. In the current example, three different IGA meshes, namely 100 × 50, 160 × 80 and 200 × 100, are considered for the later optimization of the cantilever beam.
In Figure 24, three smooth cantilever beam designs are given based on the settings in Section 5.3, and three binary designs for the cantilever beam are presented in Figure 25. As can be easily seen in Figures 24 and 25, the ITO method with the ersatz material model can effectively optimize reasonable topologies with the complete full-loading path in the

Discussion of the Influence of the IGA Mesh
In the current subsection, the key focus is to address the influence of the IGA mesh on the effectiveness of the threshold projection in the ITO method, and also to discuss the effect of the IGA mesh on the numerical error using the immersed representation model of the extended DDF. In this subsection, the ersatz material model is used in the ITO method. The other design parameters remain consistent with the second row of Table 1 and with the change in the IGA mesh. In the current example, three different IGA meshes, namely 100 × 50, 160 × 80 and 200 × 100, are considered for the later optimization of the cantilever beam.
In Figure 24, three smooth cantilever beam designs are given based on the settings in Section 5.3, and three binary designs for the cantilever beam are presented in Figure 25. As can be easily seen in Figures 24 and 25, the ITO method with the ersatz material model can effectively optimize reasonable topologies with the complete full-loading path in the design domain, even if the IGA mesh is coarse in the design, with dimensions of 100 × 50. Hence, we can state that the coarser or denser IGA mesh has no influence on the ability of the threshold projection to remove intermediate densities and force the densities to be 0 or 1. In the FEM-based three-field SIMP method, the finite element mesh also has no effect on the threshold projection, which has an effect on the ability of the SIMP method to seek the optimized topologies. Furthermore, the related numerical errors of the structural compliance are also presented in Figures 24 and 25. We can easily infer that the numerical error gradually reduces with the increase in the IGA mesh, and this feature is present in the optimization of both smooth and binary designs. The main reason for this is that the denser IGA mesh can effectively enhance the numerical precision in the analysis and lower the evaluation errors in the structural compliance, which originates from the immersed representation model of the extended DDF.

Stiffness-Maximization Designs in 3D
In this subsection, the main intention is to discuss the utility of the proposed ITO method with the extended DDF and ersatz material model in the optimization of 3D structures. In Figure 26, three different 3D structures with loads and boundary conditions are

Stiffness-Maximization Designs in 3D
In this subsection, the main intention is to discuss the utility of the proposed ITO method with the extended DDF and ersatz material model in the optimization of 3D structures. In Figure 26, three different 3D structures with loads and boundary conditions are

Stiffness-Maximization Designs in 3D
In this subsection, the main intention is to discuss the utility of the proposed ITO method with the extended DDF and ersatz material model in the optimization of 3D structures. In Figure 26, three different 3D structures with loads and boundary conditions are defined, including Michell structure, L-type structure, and Bridge-type structure. The related optimization parameters of the three benchmarks are provided in Table 3, namely, the orders of the NURBS basis functions in three parametric directions, the number of control points or control-design variables, the number of IGA elements, and the maximum material consumption.
x FOR PEER REVIEW 23 of 28   Figure 27, the optimized designs of the Michell structure, L-type structure, and Bridge-type structure are provided, including the smooth structural topologies that can be manufactured and the nearly pure 0-1 designs, which show the effectiveness of the current ITO method with the extended DDF. The convergent histories of the structural compliance, the parameter Beat, the degree of gray scale, and the numerical error of the Michell structure are provided in Figure 28. As can easily be observed, the iterative trajectories feature smooth variations and a stable convergence process. Some of the intermediate topologies of the Michell structure are presented in Figure 29. The role of the threshold projection in the earlier iterations focused on generating a complete loading-transmission path. In the latter iterations, the increase in the parameter Beta in the threshold projection gradually transformed the smooth design into a wavy topology with a 0-1 density layout.  In Figure 27, the optimized designs of the Michell structure, L-type structure, and Bridge-type structure are provided, including the smooth structural topologies that can be manufactured and the nearly pure 0-1 designs, which show the effectiveness of the current ITO method with the extended DDF. The convergent histories of the structural compliance, the parameter Beat, the degree of gray scale, and the numerical error of the Michell structure are provided in Figure 28. As can easily be observed, the iterative trajectories feature smooth variations and a stable convergence process. Some of the intermediate topologies of the Michell structure are presented in Figure 29. The role of the threshold projection in the earlier iterations focused on generating a complete loading-transmission path. In the latter iterations, the increase in the parameter Beta in the threshold projection gradually transformed the smooth design into a wavy topology with a 0-1 density layout.

The Design of a Force Inventor
In this subsection, the key intention is to demonstrate the effectiveness and utility of the currently developed ITO method in the optimization problem of the compliant mechanism. A classic force-inverter design problem is considered here, and the corresponding loads and boundary conditions are presented in Figure 30. In Table 4, the related optimization parameters are defined, including the orders of the NURBS basis functions, the structural sizes, the number of control points, the number of IGA elements, and the maximum material volume fraction. Meanwhile, in the calculation of the IGA stiffness matrix, the input and output should add the spring stiffness, which should be equal to 0.1.
Michell structure are provided in Figure 28. As can easily be observed, the iterative trajectories feature smooth variations and a stable convergence process. Some of the intermediate topologies of the Michell structure are presented in Figure 29. The role of the threshold projection in the earlier iterations focused on generating a complete loading-transmission path. In the latter iterations, the increase in the parameter Beta in the threshold projection gradually transformed the smooth design into a wavy topology with a 0-1 density layout.

The Design of a Force Inventor
In this subsection, the key intention is to demonstrate the effectiveness and utility of the currently developed ITO method in the optimization problem of the compliant mechanism. A classic force-inverter design problem is considered here, and the corresponding

The Design of a Force Inventor
In this subsection, the key intention is to demonstrate the effectiveness and utility of the currently developed ITO method in the optimization problem of the compliant mechanism. A classic force-inverter design problem is considered here, and the corresponding imum material volume fraction. Meanwhile, in the calculation of the IGA stiffness matrix, the input and output should add the spring stiffness, which should be equal to 0.1.  Figure 30. Force-inverter design with boundary and loading conditions. Figure 30. Force-inverter design with boundary and loading conditions. In Figure 31, the optimized structural topologies with the smooth and 0-1 feature are both provided, including the distribution of the control-design variables, the DDF, and the corresponding topology. As can easily be observed, there is a limited number of intermediate densities in the control-design variables, which can be applied to construct a transition area from 0 to 1 in the DDF that features smooth boundaries. In the smooth design of the force inverter, the degree of gray scale is equal to 7.04%, and the value of the numerical error is equal to −4.1%. In the final 0-1 design of the force inverter, the control-design variables almost approach the 0-1 density distribution, the degree of gray scale is equal to 0.096%, and the corresponding numerical error is equal to −0.0028%. In conclusion, the developed ITO method with the extended DDF and ersatz material model considering the topology information has a superior ability to solve problems.
14, x FOR PEER REVIEW 25 of 28 In Figure 31, the optimized structural topologies with the smooth and 0-1 feature are both provided, including the distribution of the control-design variables, the DDF, and the corresponding topology. As can easily be observed, there is a limited number of intermediate densities in the control-design variables, which can be applied to construct a transition area from 0 to 1 in the DDF that features smooth boundaries. In the smooth design of the force inverter, the degree of gray scale is equal to 7.04%, and the value of the numerical error is equal to −4.1%. In the final 0-1 design of the force inverter, the controldesign variables almost approach the 0-1 density distribution, the degree of gray scale is equal to 0.096%, and the corresponding numerical error is equal to -0.0028%. In conclusion, the developed ITO method with the extended DDF and ersatz material model considering the topology information has a superior ability to solve problems.

Concluding Remarks
The main intention of the current work was to address three important problems in previous works on topology optimization: (1) The necessity of using isogeometric analysis to replace the finite element method in optimization; (2) the overestimation problem of structural stiffness in density-based methods using the material penalization model; and (3) the segmented representations of smooth and binary topologies. Hence, a promising method of isogeometric topology optimization (ITO), with powerful effectiveness and high efficiency, is proposed. The main conclusions obtained based on several 2D and 3D numerical examples are as follows: • IGA with higher numerical precision than the conventional FEM plays a significant

Concluding Remarks
The main intention of the current work was to address three important problems in previous works on topology optimization: (1) The necessity of using isogeometric analysis to replace the finite element method in optimization; (2) the overestimation problem of structural stiffness in density-based methods using the material penalization model; and (3) the segmented representations of smooth and binary topologies. Hence, a promising method of isogeometric topology optimization (ITO), with powerful effectiveness and high efficiency, is proposed. The main conclusions obtained based on several 2D and 3D numerical examples are as follows: • IGA with higher numerical precision than the conventional FEM plays a significant role in optimization by maintaining the form of a full loading-transmission path in the topology.

•
The extended DDF with the ersatz material model can be more beneficial by lowering the overestimation of the structural stiffness resulting from the material penalization model.

•
The development of the extended DDF to unify the representations of binary and smooth designs can offer more benefits for optimization.

•
In the ITO method with the ersatz material model, the threshold projection in the extended DDF can effectively support the generation of an appropriate topology.