A Comprehensive Review of Isogeometric Topology Optimization: Methods, Applications and Prospects

Topology Optimization (TO) is a powerful numerical technique to determine the optimal material layout in a design domain, which has accepted considerable developments in recent years. The classic Finite Element Method (FEM) is applied to compute the unknown structural responses in TO. However, several numerical deficiencies of the FEM significantly influence the effectiveness and efficiency of TO. In order to eliminate the negative influence of the FEM on TO, IsoGeometric Analysis (IGA) has become a promising alternative due to its unique feature that the Computer-Aided Design (CAD) model and Computer-Aided Engineering (CAE) model can be unified into a same mathematical model. In the paper, the main intention is to provide a comprehensive overview for the developments of Isogeometric Topology Optimization (ITO) in methods and applications. Finally, some prospects for the developments of ITO in the future are also presented.


Introduction
Structural optimization [1] has attracted considerable attentions among researchers ranging from theoretical research to engineering applications, which aims to solve the optimal design of the load-carrying structures with the reasonable structural features, like the connectivity of holes, the shapes of boundaries. Overall speaking, structural optimization mainly contains three components as far as the design stage, presented in Figure 1, namely the conceptual design stage of Topology Optimization (TO), the basic design stage of shape optimization and the detailed design stage of size optimization. One of them, TO, has been identified as an important but with more challenges sub-discipline, and the main intention of TO is to seek for the optimal material layout with the expected structural performance in a design domain without the prior knowledge subject to several pre-defined constraints [2].
As we know, TO originates from a pioneering work [3] that discusses the frame-structures under the limits of economy of materials. Cheng and Olhoff [4,5] addressed the optimal design of solid elastic plates, which is considered as the seminar work for the structural optimization of continuum structures and attracts a wide of discussions in the last three decades. In 1988, Bendsøe and Kikuchi [6] used the homogenization approach to optimize the structural topology by gradually changing the sizes and orientations of holes in a design domain. After that, TO has accepted a myriad of discussions ranging from the developments of TO methods to the applications of different problems, and the details can refer to some comprehensive reviews of TO [7][8][9][10][11]. Up to now, there are several different topology optimization methods with the unique positive features which have been proposed in recent years.
The developed TO methods can be mainly divided into two branches as far as the representation model of the structural topology, including Material Description Models (MDMs) and Boundary Description Models (BDMs). In the first branch of TO methods, MDMs discrete the design domain to be a series of designable points or elements with the densities, namely the density-based TO methods. The density in each designable point or element determines the non/existence of material at the corresponding location in a design domain. This branch mainly contains the Solid Isotropic Material with Penalization (SIMP) method [12,13], and the Evolutionary Structural Optimization (ESO) method [14]. However, the second branch of TO methods uses the BDMs to display the structural topology, where a higher-dimensional function in an implicit or explicit form is constructed for the evolvement of topology in the design and structural boundaries are defined by the iso-contour/surface of the function. In this branch, the Level Set Method (LSM) [15][16][17], the phase field method [18,19], the recently proposed Moving Morphable Components/Voids (MMC/V) method [20][21][22][23] and the Bubble method [24,25] have been obtained considerable discussions. These developed TO methods have been also applied to address several different numerical problems, like the dynamic optimization [26][27][28], compliant mechanisms [29,30], stress problems [31][32][33], robust designs [34][35][36], materials design [37][38][39][40][41], concurrent topology optimization [42][43][44][45][46][47][48].
In the previously mentioned TO works, the classic Finite Element Method (FEM) [49] is applied to solve the unknown structural responses in numerical analysis. However, it is known that the FEM features several deficiencies in numerical analysis, like (1) the finite element mesh is just an approximant of the structural geometry, rather than the exact representation; (2) The neighboring finite elements have the loworder (C0) continuity of the structural responses, and the deficiency also exists in the higher-order finite elements; (3) The lower efficiency to gain a high quality of the finite element mesh. These drawbacks mainly stem from the use of different mathematical languages in geometric model and numerical analysis model: spline basis functions are used in the former whereas Lagrangian and Hermitian polynomials in the latter. Meanwhile, in TO, the optimized designs generally need the additional post-processing to meet the requirements of the practical engineering structures, so that the communication with CAD systems is compulsory. On the other side, these three deficiencies might cause the high possibility of the occurrence of numerical issues in TO. Recently, a promising and powerful alternative of the FEM, termed by the IsoGeometric Analysis (IGA), is proposed by Hughes and his co -workers [50,51] to perform the numerical analysis, which can completely remove the above limitations of FEM. In IGA, the core is that the same spline information including control points and spline basis functions is simultaneously applied into the representation of the structural geometry and solve the numerical analysis. The geometrical model and numerical analysis model are kept consistent in IGA. This such unification of the mathematical model in structural geometry and numerical analysis can offer benefits for the later optimization to resolve the above numerical issues occurred in TO.
Since the developments of IGA to eliminate the defects of the conventional FEM, several researchers have devoted to developing new TO methods and discussing their applications using IGA, rather the FEM. To the best knowledge of the authors, the first work introducing IGA into topology optimization might go back to Ref. [52], which discussed the shape optimization using IGA and its extension to the topological design. Later, an extensive work [53] used the trimmed spline surfaces to present structural boundaries and then proposed a novel Isogeometric Topology Optimization (ITO) framework based on TO and IGA, which opens up a new window for the development of TO in the future. After that, many research works have been performed to sufficiently consider the positive features of IGA into TO, which can develop more and more efficient and effective ITO methods for many numerical problems. Up to now, two publications present reviews for the IGA into structural optimization [54,55]. However, these two papers mostly focus on the descriptions about the introducing of IGA into shape optimization and its developments, and the discussions about the IGA into topology optimization are limited in these papers. Moreover, the considerations of IGA to replace the classic FEM in TO have obtained more and more attentions among many researchers in recent years. It is compulsory to provide an overview for the developments of ITO methods and their applications, which can provide more better research orientations and (2020) 33:87 suggestions for the newcomer in the field of TO or ITO, also other readers who have interests in this field. The rest of this paper is organized as follows: a brief description about the ITO methods in different types is presented in Section 2, and Section 3 provides the discussions about the applications of the ITO methods in different numerical problems. In Section 4, some prospects about the ITO in methods and applications are also presented. Finally, the paper ends with some concluded remarks in Section 5.

Isogeometric Topology Optimization (ITO) Methods
As already discussed in Introduction, Seo et al. [52,53] firstly implemented the ITO using the trimmed spline surfaces and IGA, where the trimmed surface analysis treats topologically complex spline surfaces using trimming information provided by CAD systems and it is also used for calculating structural response analysis and sensitivity calculation in TO. The spline surface and trimming curves are applied to represent the outer and inner boundaries of geometrical design models, in which the coordinates of control points of a spline surface and those of trimming curves work as design variables in TO.
In the design, this ITO framework deal with the inner front creation and inner front merging. When considering the complex structures in the optimization, the number of the trimming curves will increase, and a highly prohibitive computational cost might be caused. After that, the development of ITO starts to focus on how to construct a more efficient and effective ITO method based on the previous TO methods and IGA. Up to now, many different ITO methods have been developed. According to the classifications of TO methods already discussed in Introduction, we still divide discussions about ITO methods into two different branches, namely MDMs-based ITO methods and BDMs-based ITO methods. In the first branch using MDMs, the development of ITO methods strongly depends on the "density", namely the density-based ITO methods. As far as the second branch using the BDMs, the previous research works mostly develop ITO methods using the level set or MMC/V, namely the level set-based ITO methods and MMC/V-based ITO methods. Hence, we will provide the detailed discussions about the ITO methods in three different types, including the density-based, level set-based and MMC/V-based, respectively.

Density-Based
As we know, the homogenization approach is earlier used to realize the optimization of structural topology, which will introduces several numerical difficulties in the design. After that, several improvements are also discussed. One of them, the Solid Isotropic Material with Penalization (SIMP) method, can be viewed as a powerful alternative, which has accepted more and more attentions owing to its conceptual clarity and easy numerical implementation [12,13]. The basic intention of topology optimization to search the continuous material distribution is fully converted into seeking for the reasonable spatial arrangement of densities of finite elements. It is wellknown that some numerical artifacts are also occurred in the optimized solutions, like the checkerboards, "zig-zag" or wavy structural boundaries and mesh-dependency [56][57][58], and several works reveal that these issues mainly stem from the strong dependency on finite elements in SIMP method [59][60][61]. Hence, some alternative variants of SIMP are also developed to eliminate the numerical difficulties and produce the distinct material interface, like introducing the densities at elementary nodes [62][63][64]. A comprehensive review about the SIMP method can refer to [9,65]. Here, we provide a general mathematical model of the SIMP as far as the classic compliance minimization problem, given as: where c is the objective function, defined by the structural compliance, ρ is a vector containing a series of design variables, namely the element densities. ρ e denotes the e th element density, and p is the penalty parameter to enforce element densities to be 0 or 1. U is the global displacement field and K is the global stiffness matrix. u e is the element displacement, and k 0 is the element stiffness matrix. υ 0 is the elementary volume fraction, and V 0 is allowable material volume fraction. ρ min is the minimal value of design variables. N is the total number of element densities. Hence, in the SIMP method, the design aims to find a reasonable layout of ρ in the design domain with the expected structural compliance c, subject to the material volume fraction V 0 .
In 2011, Kumar and Parthasarathy [66] constructed B-spline finite elements for the density representation function and the displacement field in the design domain to eliminate the numerical artifacts of traditional elements, who reveal B-spline basis functions feature a smoothing effect to remove the mesh dependency, similarly to the density filtering schemes. Later, Hassani et al. [67] firstly developed an ITO method for structural , compliance problem, where densities are defined at control points and Non-Uniform Rational B-Spline (NURBS) basis functions are combined with the pre-defined densities at control points to develop the density distribution function for the representation of structural topology. As shown in Figure 2, we provide some numerical results.
We can easily find that although several numerical artifacts of SIMP can be successfully removed using the current ITO method, some new deficiencies are also shown in the optimized designs, like the blur and wavy structural boundaries. In the viewpoint of the authors, the current work opens up the combination of SIMP and IGA, which verify the feasibility of the introducing of IGA into SIMP. However, several new numerical artifacts are introduced. Meanwhile, the work directly employs the densities at control points to approximately represent the structural topology. It is suitable for the rectangular design domain, but it might introduce errors in the optimization for the curved structures. The main cause is that some parts of the control points are not at the curved design domain.
After that, Qian [68] developed a B-spline space for the topology optimization. In this work, an arbitrarily shaped design domain is embedded into a rectangular domain, which can sufficiently employ the tensor-product feature of B-splines to develop the density field for the representation of the structural topology in the design domain. The author reveals that the B-spline representation of the topology can offer an intrinsic filter for the topology optimization, which can effectively remove numerical artifacts and control minimal feature length in the optimized designs. Moreover, the B-spline space can decouple the representation of the density distribution from the finite element analysis, which avoids the re-meshing of the design domain in the multi-resolution. We also provide some numerical results in Ref. [68], as shown in Figure 3. As we can easily see, the structural features are similar to numerical results of SIMP, like the "zig-zag" or wavy boundaries. The main reason is that the final representation of the structural topology is still based on element densities which are defined by the B-spline density representation using control densities. The spatial distribution of element densities in the design domain has the intrinsic feature, namely "zig-zag". Meanwhile, the mapping from the densities at control points to element densities will increase the existence of intermediate densities. Hence, the structural boundaries of the optimized topology are still featured with a "zig-zag" or wavy shape, which still need the additional post-processing to smooth structural boundaries for the latter manufacturing.
Later, Gao et al. [69] constructed an enhanced density distribution function to develop a new ITO method. In the construction of the density distribution function, two steps are involved: (1) Smoothness: the Shepard function is firstly employed to improve the overall smoothness of the densities pre-defined at control points. (2) Continuity: the NURBS basis functions are linearly combined with the smoothed control densities to construct the density distribution function. In each optimization iteration, the density distribution function to represent the structural topology will be advanced. As shown in Figure 4, some numerical results are also given. As we can see, an enhanced density distribution function can offer more benefits for the optimization and the representation of the structural topology. However, it should be noted that the structural boundaries are represented by the iso-contour of the density distribution function with the iso-value (0.5) of the density. It originates from the (2020) 33:87 level set method, and the reasonability of the definition of the structural boundaries at the iso-contour/surface of the density distribution function is also discussed. We can easily find that the post-processing scheme is very simple, heuristic but efficient. However, it also introduces some errors in the evaluation of structural performance of the optimized designs. Lieu and Lee [70] developed a multiresolution scheme to topology optimization using the framework of IGA, where a variable parameter space is defined for the implementation of multiresolution TO using SIMP method. Then, they inherited the multiresolution ITO framework [70], and applied it to discuss the multi-material topology optimization problem [71], in which the alternating active-phase algorithm [72] for the multi-material topology optimization is directly used in the multiresolution ITO framework. Wang et al. [73] discussed the multiscale ITO for periodic lattice materials, in which the asymptotic homogenization is applied for the calculation of mechanical properties for lattice materials with uniform and graded relative density respectively. Taheri et al. [74] also studied the application of the ITO to the multimaterial topology optimization problem and the design of functionally graded structures, where the multi-material interpolation scheme proposed by Stegmann and Lund [75] to realize the discrete material optimization is directly used. Liu et al. [76] also addressed the stressconstrained topology optimization problem of plane stress and bending of thin plates using the ITO framework, where two stability transformation methods are developed to stabilize the optimization using the P-norm function for global stress constraint. Later, Gao et al. [77] proposed a NURBS-based Multi-Material Interpolation (N-MMI) model in the ITO method [69] to develop a Multi-material ITO (M-ITO) method. Then Gao et al. [78] employed the ITO method to study the design of auxetic metamaterials and the M-ITO method to discuss the optimization of auxetic composites, where a series of novel and interesting material microstructures with the auxetic property can be found. Xu et al. [79] also applied the ITO method to study the rational design of ultralightweight architected materials. The topology optimization of the spatially graded hierarchical structures is also discussed in the framework of ITO [80]. Xie et al. [81] also proposed a truncated hierarchical B-spline-based topology optimization to address topology optimization for both minimum compliance and compliant mechanism. Wang et al. [82] discussed the numerical efficiency of the ITO method and employed the multilevel mesh, MGCG and local-update strategy to improve the computational efficiency by mesh scale reduction, solving acceleration and design variables reduction. Zhao et al. [83] also addressed the T-Splines Based ITO method for the design domains with arbitrarily shape, where the arbitrarily shaped design domains is directly obtained from CAD and defined by a single T-spline surface. The T-spline can overcome the topological limitations of NURBS. However, it also introduces an important problem that how many control points should be arranged in the local structural features. The basic feature of TO is that we do not known the final optimized design without the prior knowledge. Hence, a uniform initial design is much better for the latter optimization, which can offer the equal opportunity for the advancement of each point in the design domain and avoid the occurrence of the local optimum design. However, when using T-splines to model the geometry and analysis, a non-uniform IGA mesh will occur and also a control lattice with nonuniform features will be utilized, which will introduce some numerical issues in the latter optimization.

Level Set-Based
It is known that Level Set Method (LSM) is numerical technique to track the interface and shape, which has been extensively used in many disciplines. The core of the LSM is to define a level set function with a higher-dimension to represent the structure, where the zero-level set is employed to represent the structural boundaries. The level set function with the negative values are applied to display the voids, and the solids in the design domain are represented by the level set function with the positive values, namely the implicit boundary representation model. Hence, the evolvement of the level set function can describe the advancing of the structural topology in the design domain.
As already discussed in Section 1, Sethian and Wiegmann [15] firstly employed the level set function to represent the structure topology and used structural stress to develop the evolving mechanism. After that, Wang et al. [16] innovatively developed the level-set topology optimization framework, where the upwind scheme and the finite difference method are utilized to solve the H-J PDEs to advance the structural topology. Allaire et al. [17] developed a level-set topology optimization method based on the classical shape derivatives in the level-set method for front propagation. Compared to MDMs, we can easily find that the level-set topology optimization is actually a shape optimization method but with a superior capability to implement the shape and topology optimization. The optimized topologies will have the smooth structural boundaries and distinct interfaces, and the LSM will feature several inherent physical merits: (1) a smooth and distinct boundary description for the optimized design, (2) the shape fidelity and higher topological flexibility during the optimization, (3) the shape and topology optimization are performed simultaneously and (4) a physical meaning solution of the H-J PDEs. The mathematical model of the level-set based TO method for the structural compliance problem can read as: where J is the objective function, defined by the structural compliance problem. u denotes the global displacement field in design domain, and Φ is the level set function with a higher dimension to represent the structural topology. D is the reference domain, and Ω is the design domain containing all admissible shapes. H is the Heaviside function which serves as a characteristic function. G is the volume constraint function. V 0 is the allowable material consumption. The elastic equilibrium equation is stated in the weak variational form, in which a is the bilinear energy function and l is the linear load function. υ is the virtual displacement field, which belongs to the kinematically admissible displacement space U. As shown in Figure 5, a 3D level set function with the corresponding 2D structural design domain is given.
In 2012, Shojaee et al. [84] discussed the composition of IGA with LSM to develop a level set-based ITO framework for the structural topology optimization, where the Radial Basis Function (RBF) is applied to parametrize the level set function. The corresponding numerical results are shown in Figure 6(a). In Ref. [84], the level set function is constructed by the RBF to show the topology, and IGA uses the NURBS basis functions to develop the analysis model. In actual, we can easily obtain that the geometric model and analysis model are not in an integrated mathematical language. , Later, Wang et al. [85] also proposed a parametrized level set-based ITO method using parametric level set method and IGA, where the same NURBS basis functions are used to parameterize the level set function and construct the solution space of numerical analysis. The geometric model and the analysis model of the structural topology can be unified, which coincides with the core of IGA. The numerical results of [85] are also presented in Figure 6(b). Then, Wang et al. [86] discussed the topology optimization for geometrically constrained design domains using the proposed level set-based ITO method, where the fast point-in-polygon algorithm and trimmed elements are utilized for ITO with the arbitrary geometric constraints. As shown in Figure 6(c), the corresponding numerical results are also given. Xia et al. [87] implemented Graphics Processing Units (GPU) parallel strategy for the level setbased ITO method to improve numerical efficiency. After that, Ghasemi et al. [88] also developed a level set-based ITO method but for the optimization of flexoelectric materials, where the NURBS-based IGA elements are successfully employed to model the flexoelectric effect in dielectrics and the energy conversion efficiency of flexoelectric micro and nanostructures is improved. Moreover, the point wise density mapping is directly used in the weak form of the governing equations and the adjoint sensitivity technique is applied to compute the derivative. Jahangiry et al. [89] also discussed the application of IGA in the structural level set topology optimization to develop a new level set-based ITO framework, where the control mesh is gradually updated in the optimization iterations, and then the authors also discussed the application of the new level set-based ITO framework in the topology optimization of the concentrated heat flow and uniformly distributed heat generation systems [90]. Lee et al. [91] also implemented the isogeometric topological shape optimization using dual evolution with boundary integral equation and level sets, where the implicit geometry  (2020) 33:87 using the level sets is transformed into the parametric NURBS curves by minimizing the difference of velocity fields in both representations. Xu et al. [92] employed the level set-based ITO method in Ref. [85] to discuss the design of vibrating structures to maximize the fundamental eigenfrequency and avoid resonance, and the related numerical results are shown in Figure 7(a). Yu et al. [93] also employed the level set-based ITO method in Ref. [85] to implement the multiscale topology optimization using the unified microstructural skeleton, where a prototype microstructure is defined to obtain a series of graded microstructures. Figure 7(b) shows the related numerical results. In Ref. [94], a level set-based ITO method was proposed for topology optimization to control the high-frequency electromagnetic wave propagation in a domain with periodic microstructures, where the high-frequency homogenization method is used to characterize the macroscopic highfrequency waves in periodic heterogeneous media. The corresponding numerical results are also given in Figure 7(c).

MMC/V-Based
Compared to the density-based and level set-based TO methods, MMC/V has implemented the topology optimization in an explicit and geometrical way. MMC/V can incorporate more geometry and mechanical information into topology optimization directly. Since the seminar work of MMC proposed by Guo et al. [20], it have been accepted more and more attentions in not only theoretical research but also engineering applications. Zhang et al. [95] developed a new MMC-based topology optimization method, where the ersatz material model is utilized through projecting the topological description functions of the components. Later, Guo et al. [21] studied the explicit structural topology optimization based on moving morphable components (MMC) with curved skeletons. In Refs. [22,23], the B-spline curves are used to describe the boundaries of moving morphable voids (MMVs) to develop the MMV-based topology optimization method. In 2017, Hou et al. [96] firstly proposed an MMC-based ITO method, where NURBS basis functions are applied to construct the NURBS patch to represent the geometries of structural components using explicit design parameters and the same functions are also applied into the latter IGA. As already indicated in Ref. [96], the proposed MMC-based ITO method can naturally inherit the explicitness of the MMC-based TO method, and also embraces the merits of IGA, such as a tighter link with Computer-Aided Design (CAD) and higher-order continuity of the basis functions. The numerical results are displayed in Figure 8(a). Xie et al. [97] also developed a new MMC-based ITO method based on R-functions and collocation schemes, in which the R-functions are used to construct the topology description functions to overcome the C1 discontinuity problem of the overlapping regions of components. As given in Ref. [97] to discuss the efficiency of the proposed method, the numerical results show that the current method can improve the convergence rate in a range of 17%-60% for different cases in both FEM and IGA frameworks. This proposed MMC-based ITO method was applied to the topology optimization for the symmetric structures using energy penalization method [98]. After that, Xie et al. [99] proposed a new MMC-based ITO method using a hierarchical B-spline which can implement the adaptive IGA to efficiently and accurately assess the structural performance. As far as the MMV-based ITO method, Zhang (2020) 33:87 et al. [100] proposed a new MMV-based ITO method, in which the MMV-based topology optimization framework is seamlessly integrated into IGA by using TSA (trimming surface analysis) technique. Comparatively speaking, the current MMV-based ITO method can flexibly control the structural geometry/topology. Meanwhile, it can also prevent the occurrence of self-intersection and jagged boundaries. The related numerical results are also shown in Figure 8(c). Later, Gai et al. [101] also studied the development of the MMV-based ITO method, where the closed B-spline boundary curves are utilized to model the MMVs to represent the structural topology. Du et al. [102] discussed the application of the MMC-based ITO method in the multiresolution topology optimization problem.

Other Types
Besides the previously mentioned works, the ITO methods are also developed based on other TO methods. Dedè et al. [103] proposed a phase field-based ITO method, where the optimal design can be obtained by the steady state of the phase transition described by the generalized Cahn-Hilliard equation. The numerical solutions are presented in Figure 9(a). Yin et al. [104] developed an ITO method based on the scheme of Bi-directional Evolutionary Structural Optimization (BESO), namely the BESO-based ITO method. Sahithi et al. [105] studied the evolutionary algorithms to realize the ITO of continuum Structures using the parallel computing, where the evolutionary optimization process and metaheuristics are used to optimize the layout of material in the design domain, and the related numerical results are shown in Figure 9(b).

Applications of ITO
In Section 2, we give a comprehensive review about the development of the ITO methods considering three components: the density-based ITO, level set-based ITO and MMC/V-based ITO. In the development of the ITO methods, the applications of the ITO methods are also involved into many numerical optimization problems, like the classic structural compliance problem with the single-material [67-69, 85, 89, 96, 97], the multi-material topology optimization problem [71,74,77,106], the trimmed spline surfaces [53,86,107], the functional graded structures [74,80].
In this section, we review the applications of the ITO in three important numerical optimization problems, including mechanical metamaterials, the splines used in IGA and the computational efficiency.

Mechanical Metamaterials
Mechanical metamaterials are a kind of artificial materials with counterintuitive mechanical properties that are obtained by the topology of their unit cell instead of the properties of each component [108]. Generally speaking, mechanical metamaterials are always associated with four elastic constants: Young's modulus, shear modulus, bulk modulus and Poisson's ratio. The corresponding subtypes of mechanical metamaterials mainly contains acoustic metamaterials, auxetic metamaterials, etc.
As already discussed in the definition of mechanical metamaterials, the effective macroscopic properties of materials strongly depend on the micro-architecture that are homogeneously arranged in the bulk material, rather than constituent properties of the base material. This feature of mechanical metamaterials can offer the high possibility for the applications of topology optimization to seek for a series of novel metamaterial microstructures with the promising macroscopic properties. Since the homogenization theory is developed to predict macroscopic effective properties [109], an inverse homogenization procedure is proposed for the optimization of a base unit cell with the negative Poisson ratio using topology optimization [110]. Later, this work is inspired and extended to the topology optimization of the rationally artificial materials with the extreme or novel properties [111], particularly for auxetic metamaterials with the Negative Possion's Ratio (NPRs) behavior.
The earlier work introducing the IGA into the design of mechanical metamaterials can go back to Ref. [112], in which the IGA-based shape optimization is developed for the design of smoothed petal auxetic structures via computational periodic homogenization. The authors also discussed the optimal form and size feature of planar isotropic petal-shaped auxetic structures with the tunable effective properties using the IGA-based shape optimization [113]. The IGA-based shape optimization for periodic material microstructures using the inverse homogenization was also studied in Ref. [114]. The introducing of IGA into topology optimization for the rational design of auxetic metamaterials can track to Ref. [78], which used the SIMPbased ITO method proposed in Ref. [69] and also numerically implemented the energy-based homogenization method to evaluate the effective macroscopic (2020) 33:87 properties using IGA with the imposing of the periodic boundary formulation on the base material unit cell. A reasonable ITO formulation for auxetic metamaterials with the re-entrant and chiral deformation mechanisms is developed, and several optimized design are shown in Figure 10(a). Later, the authors also discussed the computational design of auxetic composites via an IGA-based M-ITO method developed in Ref. [77], where an appropriate objective function with a weight parameter is also defined for the controlling of the generation of different deformation mechanisms with the re-entrant and chiral in auxetic composite microstructures [115]. The related numerical optimized microstructures with the auxetic are also shown in Figure 10(b). Later, Nguyen et al. [116] also discussed the design of auxetic metamaterials using the level setbased ITO method, where the reduced order model is utilized to reduce the computational degree of the linearly elastic equilibrium equation to improve the computational efficiency. Similarly, a series of novel and interesting auxetic microstructures in 2D and 3D, shown in Figure 10(c). Xu et al. [79] also utilized the density-based ITO method to discuss the rational design of ultra-lightweight architected materials with the extreme bulk modulus and extreme shear modulus, and a series of novel 3D ultra-lightweight architected material microstructures can be found. Nishi et al. [94] utilized the LSM-based ITO method to discuss the design of periodic microstructures in anisotropic metamaterials to control high-frequency electromagnetic wave, in which anisotropic metamaterials with the hyperbolic and bidirectional dispersion properties at the macroscale can be obtained.

Splines
In the development of the ITO method, a key in IGA is to the spline. In the earlier ITO works, the trimmed spline surfaces are employed to represent the structural topology. The outer and inner structural boundaries of the geometry are represented by a spline surface and trimming curves, in which design variables are the coordinates of control points of a spline surface and those of trimming curves [52]. This basic numerical technique is inherited in the later work [53,100,107], where the trimmed surface analysis is employed for the structural response analysis and sensitivity calculation in the optimization. A basic numerical scheme for the merging of the inner is shown in Figure 11(a). Later, the B-spline is employed in the construction of the geometrical model and B-spline basis functions are applied to develop the solution space in the IGA. Meanwhile, the B-spline-based IGA is introduced in the topology optimization. Qian [68] constructed a B-spline space for the topology optimization, where an arbitrarily shaped domain can be embedded into a rectangular domain modelled by the tensor-product B-splines. Some researchers studied the role of the B-spline in the topology optimization without using the IGA to solve the structural responses [117], where the free-form curve of closed B-splines is introduced as basic design primitives to realize topology optimization with small number of design variables. Then, the B-spline multi-parameterization method is proposed for topology optimization of thermoelastic structures [106]. After that, the hierarchical spline is applied into the development of the MMC-based ITO method, in which the adaptive IGA is implemented by the hierarchical B-spline to efficiently and accurately assess the structural performance [99]. Xie (2020) 33:87 et al. [81] developed a truncated hierarchical B-splinebased topology optimization. It should be indicated that sensitivity and density filters with a lower bound can be adaptively consistent with the hierarchical levels of active elements to remove the checkboard pattern and reduce the gray transition area. A basic illustration of the hierarchical B-spline is given in Figure 11(b). Comparatively speaking, NURBS, working as a mathematical model commonly used in computer graphics for generating and representing curves and surfaces, is also mostly employed in the development of the ITO method in three types. Wang et al. [85] developed the level setbased ITO method using NURBS, in which NURBS is firstly applied to parametrize the level set function to represent the structural topology and then construct the solution space in IGA to solve the unknown structural responses. Gao et al. [69] also employed NURBS to develop an enhanced density distribution function with the sufficient smoothness and continuity to represent the structural topology, and the same NURBS basis functions are also used in IGA. Hou et al. [96] used NURBS to construct the MMCs for the representation of the geometries of structural components (a subset of the design domain) with use of explicit design parameters, and the NURBS-based IGA is also applied to solve the structural responses. A basis description about the NURBS for the representation of the structural geometry is shown in Figure 11(c). Besides the above discussed splines, the T-spline is also used in the ITO method for the topology optimization, and the T-spline-based ITO method is developed to realize the optimization of design domain with arbitrary shapes [83] to eliminate the complexity of the multi-patch NURBS for the structural geometry. In actual, it will introduce an important problem that how many control points should be arranged in the representation of structural local features, which will have a significant effect on the latter topology optimization.

Computational Cost
Although computer has gained a great number of developments in recent years, the computational cost of topology optimization is still a prohibitive problem, especially for the common laptop. In order to improve the computational efficiency of the ITO in numerical implementations, several research works have been implemented in recent years. The most method is the use of multiresolution scheme in numerical calculation of the topology optimization [118]. In the multiresolution topology optimization, three distinct meshes are defined for the optimization: (1) a displacement mesh for the finite element analysis; (2) a design variable mesh for the optimization; and (3) a density or level set mesh to display the material distribution. The basic idea is that topology optimization can achieve the higher-resolution designs but with a lower computation cost as well. Lieu et al. [70] developed a multiresolution ITO method using SIMP to improve computational efficiency, and then applied it to address the multi-material topology optimization problem [71]. Du et al. [102] also utilized the multiresolution scheme in the MMV-based ITO method to reduce the computational cost. A simple illustration of multiresolution scheme is shown in Figure 12. Wang et al. [82] also improved the computational efficiency in three aspects: namely the mesh scale reduction, solving acceleration and design variables reduction, and the ITO method is developed using multilevel mesh, multigrid conjugate gradient method and local-update strategy. As already given in numerical results, the current proposed method can successfully reduce 37%-93% computational time compared to previous works. The GPU parallel strategy is also employed in the parameterized LSM-based ITO method to reduce the computational cost [87], where the parallel implementations are utilized in the initial design domain, IGA, sensitivity analysis and design variable update.

Prospects
In this Section, we will provide three main directions for the development of ITO in the future, including the Data-driven ITO, ITO for additive manufacturing and ITO considering the advantages of IGA in several problems. The details are given as follows.
(1) Data-driven ITO: It is known that the application of topology optimization for the complex engineering materials is very difficult due to the complexity. In recent years, the big data and machine learning have been becoming popularly in the field of computational mechanics, which can provide new windows for the topology optimization for complex problems. For example, the deep neural network is employed to approximate the field of variables to solve the boundary value equations in strong or weak forms [119]. In the work, the developed data-driven neural networks can efficiently reduce the computational costs. Meanwhile, the datadriven isogeometric shape optimization for auxetic microstructures is also studied [120]. Hence, in the future work of the ITO, the data-driven ITO method and its applications in several numerical problems will be the promising research topic.
(2) ITO for additively manufacturing: In recent years, additive manufacturing technique, a layer-by-layer manner to fabricate structures, has accepted great attentions and been becoming a powerful alternative to the conventional fabrication methods, like the machining and casting, due to its merits to manufacture the structures with specific features, like the cavity. Hence, additive manufacturing can offer the higher flexibility and efficiency for the fabrication of structures. The topology optimization design for additive manufacturing has proposed in recent years, and the comprehensive reviews for this topic can refer to Refs. [121,122]. IGA has the positive feature to unify Computer-Aided Design (CAD) model and Computer-Aided Engineering (CAE) into a same mathematical language, so that the ITO can offer more possibility for engineering structures from the conceptual design phase to the last manufacturing into an integrated process, if the development of the ITO can consider the additive manufacturing. This unification will significantly reduce the financial cost of the product design. Hence, in the future, the ITO for additive manufacturing will be also a hot research topic.
(3) The advantages of IGA in several problems: IGA has the compelling benefits in the field of shell and plate overall conventional approaches [123,124], because the smoothness of NURBS basis functions can offer a straightforward manner to construct the plate/shell elements, particularly for the thin shells and rotation-free. Meanwhile, the smoothness of NURBS basis functions can also offer more benefits for the analysis of fluids [125] and the fluid-structure interaction problems [126]. In addition, due to the ease of construction of the higherorder basis functions, IGA with more success can be utilized to solve PDEs with the forth-order (or higher) derivatives, for example the Hill-Cahnard equation [127]. Hence, in the future, the considerations of the ITO in the mentioned above numerical problems might be more meaningful for the development of this field.

Conclusions
In the current paper, we offer a comprehensive review for the Isogeometric Topology Optimization (ITO) in methods and applications. Firstly, we mainly divide the descriptions of ITO methods into three aspects, including the density-based ITO methods, level set-based ITO methods and MMC/V-based ITO methods. The corresponding discussion for each classification is clearly provided, and the development trajectory in each classification is also given. Secondly, the descriptions of the applications of ITO mainly focus on three components, namely the ITO for mechanical metamaterials, the splines in ITO and the computational cost of ITO. Finally, we also provide some prospects for the developments of the ITO methods and applications in the future, which contains the data-driven ITO to considerably reduce the computation cost, the ITO for additive manufacturing to consider the manufacturing problems into the initial conceptual design phase and the ITO considering the advantages of IGA in several problems.