Concurrent optimization of structural topology and infill properties with a CBF-based level set method

In this paper, a parametric level-set-based topology optimization framework is proposed to concurrently optimize the structural topology at the macroscale and the effective infill properties at the micro/meso scale. The concurrent optimization is achieved by a computational framework combining a new parametric level set approach with mathematical programming. Within the proposed framework, both the structural boundary evolution and the effective infill property optimization can be driven by mathematical programming, which is more advantageous compared with the conventional partial differential equation-driven level set approach. Moreover, the proposed approach will be more efficient in handling nonlinear problems with multiple constraints. Instead of using radial basis functions (RBF), in this paper, we propose to construct a new type of cardinal basis functions (CBF) for the level set function parameterization. The proposed CBF parameterization ensures an explicit impose of the lower and upper bounds of the design variables. This overcomes the intrinsic disadvantage of the conventional RBF-based parametric level set method, where the lower and upper bounds of the design variables oftentimes have to be set by trial and error. A variational distance regularization method is utilized in this research to regularize the level set function to be a desired distanceregularized shape. With the distance information embedded in the level set model, the wrapping boundary layer and the interior infill region can be naturally defined. The isotropic infill achieved via the mesoscale topology optimization is conformally fit into the wrapping boundary layer using the shape-preserving conformal mapping method, which leads to a hierarchical physical structure with optimized overall topology and effective infill properties. The proposed method is expected to provide a timely solution to the increasing demand for multiscale and multifunctional structure design.


Introduction
Structures containing cellular infills or micro/meso architectures can possess fine-tuned properties and extra functionalities with a low density [1][2][3][4]. With the help of modern additive manufacturing technologies [5][6][7], the multiscale structures are foreseeing a wide range of engineering applications [8][9][10]. However, how to identify the exact structural layout at different scales to maximize the performance is a challenging topic.
At micro/meso scale, topology optimization of metamaterials or lattice structures has been well studied since Sigmund's pioneering work [11], where the material is optimally distributed inside the design domain to achieve desired effective material properties. The homogenizationbased topology optimization method [12][13][14], the solid isotropic material with penalization (SIMP) method [15,16], the level set method (LSM) [17][18][19], and the evolutionary method [20] are most commonly used approaches in this field. With an inverse homogenization method, Sigmund designed truss-like base cells with different material properties [21,22]. With the SIMP method, Sigmund and Torquato [23] elaborated the three-phased microstructure designs with extreme thermoelastic properties. As for the conventional LSM, Mei and Wang [24] designed multimaterial microstructures utilizing the "color" level sets. Vogiatzis et al. [25] employed a reconciled level set approach to design multimaterial microstructures with negative Poisson's ratio (NPR). As an extension to the conventional level set approach, Wang et al. [26] employed a numerically robust parametric level set method (PLSM), in order to design a series of microstructures with different Poisson's ratios. In this systematic computational design framework, the homogenization was used for calculating the material property and PLSM was used for updating the structural shape and topology. With the evolutionary approach, a series of papers have been published for designing microstructures with extreme properties [27][28][29]. A comprehensive review for the metamaterial design using evolutionary structural optimization is provided in Ref. [30].
Instead of merely designing the micro/meso structures, recent years have witnessed more and more researches on optimizing both the structure and its constructing material properties in a concurrent manner [31]. The SIMP-based concurrent structural design was reported by Rodrigues et al. [32] and Coelho et al. [33], for designing the overall structural geometry and its infill metamaterials. Liu et al. [34] uncoupled the macro and the micro densities in the SIMP method, to carry out a concurrent topology optimization of truss-like material. Deng and Chen [35] studied the concurrent topology optimization by designing structures with multiple porous metamaterials under radon field loading uncertainty. To achieve a high-performance shell-infill structure, Wu et al. [36,37] proposed a concurrent SIMP-based topology optimization with a local volume constraint to design bone-like structures. The SIMP-based concurrent approach has also been extended to designs utilizing thermal elastic materials [38]. Within the level set framework, Sivapuram et al. [39] concurrently optimized the overall topology and the metamaterial properties with a decomposition formulation for an easy algorithm parallelization. Wang et al. [40] reported a concurrent multiscale design with a shape metamorphosis technology, aiming at ensuring the connectivity of the adjacent microstructures. Li et al. [41] proposed an integrated multiscale structure design, where SIMP is used for macroscale structure design and a radial basis functions (RBF)-based PLSM is used for microscale structure design. Li et al. [42] further developed the multiscale structure design based on the concurrent PLSM approach, in which the connectivity between the adjacent microstructure patches were ensured by virtual kinematic connectors. Via the evolutionary approach, Xia and Breitkopf [10] designed the overall structural topology and local material property concurrently, using the multilevel finite element approach (FE 2 ) to handle the material nonlinearity. Readers interested in evolutionary approaches to concurrent designs of multiscale structures can be referred to Refs. [43][44][45][46] for further details. In this paper we employ the level set approach and improve it to make it suitable for multiscale topology optimization.
With an implicit design representation, the designs represented by the level set model has clear design boundaries. Mathematically, the level set model can possess extra geometrical information such as curvatures, and can handle topological changes in a robust way. In the conventional level-set-based topology optimization, the Hamilton-Jacobi partial differential equation (PDE) describes the dynamics of the boundary motion [47]. The design velocity field derived from shape sensitivity analysis [19,48,49] drives the design boundary to its optimum location. However, the conventional LSM has several numerical limitations [50][51][52]. Besides, introducing multiple constraints to the design can be challenging when the Lagrange multiplier method is used. When utilized for multiscale structure design, introducing extra design variables can be difficult for the conventional level set approach. To increase the numerical efficiency and robustness of conventional LSM, Wang et al. [53] proposed PLSM with the RBF [54]. By utilizing the method of moving asymptotes (MMA) [55] as the optimizer, numerical efficiency can be improved and design constraints can be imposed in a straightforward way. The PLSM provides the flexibility in parameterizing the structural topology as well as the material properties, which enables the concurrent topology optimization for multiscale structures. Apart from using RBF, the R-function and B-spline are also used in the level set parameterization scheme. The readers can be referred to Refs. [56][57][58] for further details.
The RBF-based PLSM approximates the level set function as the weighted summation of RBF kernel functions. As design variables, those kernel function weights are updated during the optimization. However, those kernel function weights have no exact physical meaning, and their lower and upper bounds cannot be explicitly defined. In practice, those bounds are set by trial and error. To get the explicit bounds for the design variables, a cardinal basis function (CBF) was constructed based on the RBF partition of unity collocation method [59]. The CBF satisfies the Cardinal properties which equals to 1 at the selected node and equals to 0 elsewhere. This means when CBF is used as the kernel function, the corresponding weights have the physical meaning of being the heights of the level set function. Moreover, a distance regularization energy functional [60,61] can be introduced to regularize the level set function to be a distance function in a zone along the boundary and flat surfaces elsewhere. The regularized distance function can ensure an accurate material property interpolation from the level set model to the physics model. The flat surfaces can impose preferred numerical stability and can help the creation of new holes [62]. Besides, the signed distance function can preserve the distance information directly. This means with one level set function, a shell-infill multiscale structure can be designed [63][64][65]. The shell layer of this structure can function as a protective coating against loading or extreme exterior conditions [66] while the infills can provide the designed functionalities.
Another important issue of multiscale structure optimization is to consider its manufacturability. The proposed multiscale structure is generated by mapping the optimal infill into its wrapping shell. This process will be tricky on the boundary areas, since the infill cells will be rotated or deformed. A rotated or deformed anisotropic infill cell cannot preserve its designed property. To address this issue, an isotropic infill micro/meso structure is designed to be filled into the interior infill region. The isotropy of the infill structure is ensured by imposing the isotropy constraint [67][68][69][70] to the topology optimization process. As for the mapping, a local shape-preserving conformal mapping is employed in this research. With the conformal mapping, the material properties of the designed infill structure can be mathematically preserved. The final mapped structure can be sent to 3D printer directly to be manufactured with a high fidelity of its designed material property.
To sum up, in this research, the concurrent multiscale structure design is carried out via a variational CBF-based PLSM. The infill properties are treated as design variables to be updated together with the structural design simultaneously. A unified isotropic infill metamaterial is designed to lower the computational cost for the prefabrication of 3D printing. With a modified Heaviside material property interpolation, from macroscale, the void, the wrapping shell, and the infill material region will be naturally discriminated. By constructing a user-defined distance regularization energy functional, the distance regularization effect can be directly controlled. A second stage optimization is directly carried out to find out the isotropic infill metamaterial layout. With both the overall geometry and the infill structure design in hand, the local shape-preserving conformal mapping is performed to integrate them together. To make the mapping adaptive to complicated geometries, instead of using the conventional Ricci flow with four control points [71], in this paper, a multi-control-point conformal mapping is proposed so that it is much flexible in controlling area distortions, especially for irregular domains.
The remaining paper is organized as follows. In Section 2, the strain energy method is introduced to calculate the effective material properties of a unit cell. In Section 3, the level set representation, its parameterization, and the distance regularization are elaborated. The general problem setting for the proposed concurrent optimization is formulated in Section 4. The detailed concurrent optimization for given problems are detailed in Section 5 with numerical results. The conformal mapping method is introduced in Section 6, together with mapped multiscale structure results. At last, this work is summarized in Section 7.
2 Prediction of the effective material properties based on the strain energy method For a given unit cell, noting its homogenized elasticity tensor as C H ijkl , a constitutive law can be formulated based on the Hooke's law: where ij and ε kl are the effective stress and strain tensor, respectively. Under a linear elastic problem setting, the stress and the strain tensors of the homogeneous medium should be equal to the average stress and strain in the microstructure, as ij ¼ 1 where V is the volume of the unit cell.
With a 2D plane-stress assumption, Eq. (1) can be expanded to the form of Eq. (2): On the other hand, another equilibrium equation can be built based on the fact that the strain energy stored inside both the homogeneous medium, noted as U H , and inside the unit cell, noted as U, is the same. With four special strain field as shown in Fig. 1: Four strain energies can be calculated accordingly as U 1 , U 2 , U 3 and U 4 . For example, with the first loading case, as listed in Eq. (3): For a unit cell with volume of 1, C H 1111 ¼ 2U 1 . Similarly, based on different boundary conditions, the elastic tensor can be re-written as 3 Level set representation, parameterization and the distance-regularized level set function

Level set implicit boundary representation
The level-set based implicit boundary representation [72] can be expressed as where F represents the level set function. Within the design domain D, the structure region is denoted by W and the design boundary is represented as G.
x signifies an arbitrary point in the design domain D and t is a pseudo time for the dynamic shape optimization process. This implicit boundary representation is illustrated in Fig. 2.
The dynamic boundary evolution is governed by the Hamilton-Jacobi PDE as where the V n is the normal design velocity field.

Parameterization of the level set function using a kernel function
With a given kernel function, the level set function can be written as the weighted summation form: In this study, j , representing ðx j Þ for simplicity, is the weight on the jth node. Ψ j ðxÞ (j = 1,2,..,n) is the newly constructed CBF, with the Kronecker delta property as The readers can be referred to Ref. [62] for the complete construction procedure of the CBF kernel function. A comparison is made in Fig. 3 between a regular Gaussian RBF and the CBF. One thing needs to be mentioned is that with the analytical parametric level set function expression, the optimization results generated via PLSM can be manufactured via the high-resolution layer-image-based continuous liquid interface production (CLIP) 3D printing technology [63,73]. By bridging the PLSM design methodology and the CLIP printing process, the optimal structure can be manufactured at a lower prefabrication computational cost.

Distance-regularized level set function
To maintain a distance-regularized level set evolution, a distance regularization energy functional is introduced in the current study. Generally, for a given level set function Φ, the distance regularization energy functional RðjrΦjÞ can be expressed as where PjrΦj is the regularization energy potential.
To obtain a distance-regularized level set function, Li et al. [60] proposed the "double-well" regularization potential which takes the form: The "double-well" potential PðjrΦjÞ, and its diffusive rate D ¼ d p ðjrΦjÞ are plotted in Fig. 4, where the diffusive rate has the form: It can be seen from Fig. 4(b) that this "double-well" distance regularization effect will be affected by the "selecting point" at 0.5. With the minimization of the energy R, when the jrΦj is below the selecting point 0.5, with a positive diffusive rate, jrΦj tends to be driven to 0. Oppositely, when jrΦj is above the selecting point 0.5 but below 1, with a negative diffusive rate, jrΦj tends to be driven to 1. When jrΦj is above 1, with a positive diffusive rate, jrΦj tends to be driven to 1. This will regulate the level set function Φ to a desired distanceregularized shape.
Easily, this selecting point at 0.5 can be modified based on needs. One of the constructing procedures was recently reported in Ref. [61]. By following a similar procedure and choose 0.25 as the "selecting point", a new energy potential P n used in this research can be built as As can be seen from Fig. 5, when the jrΦj is below the selecting point 0.25, with a positive diffusive rate, jrΦj tends to be driven to 0 and when jrΦj is above the selecting point 0.25 but below 1, with a negative diffusive rate, jrΦj tends to be driven to 1. When jrΦj is above 1, with a positive diffusive rate, jrΦj tends to be driven to 1. With a wider range for jrΦj to be driven to 1, the corresponding level set function will be relatively more likely to be driven to a signed distance function near the design boundary. In this case, the shell region for the multiscale structure will be well maintained.
To provide an intuitive process of imposing the distance regularization energy functional, the whole process is visualized in Fig. 6. By importing a black/white image (TO is short for topology optimization) as a binary initial level set function and minimizing the distance regularization energy according to it, the final distance-regularized level set function can be achieved, as shown in Fig. 6(c). It can be noticed that the overall "TO" shape is preserved. As can be seen from the zoom-in view in Figs. 6(d) and 6(e), the transition area of the level set function is regularized into a signed distance function while the rest areas are kept flat.
The signed distance function can be used to discriminate the wrapping shell region from the interior infill region, which can provide an accurate mapping from the level set model to the physics model.

General concurrent optimizing settings for multiscale structure design
To formulate the proposed multiscale concurrent structural topology optimization in a general way, the level set nodal height , together with the infill material properties (Young's modulus E and the Poisson's ratio v) are set to be the design variables. Within the distance regularized parametric level set framework: where the total objective function J tot contains the distance regularization energy functional R with weight w.
To discriminate a wrapping shell region and the infill region from the macroscale, the corresponding material interpolation from the level set model to the physics model needs to be modified. Conventionally, the material property interpolation is done by directly getting the Heaviside function value of the level set function. This process will generate a binary field in order to separate the material region from the void region, as shown in Figs. 7(a) and 7(b). In the proposed material interpolation scheme, a second Heaviside function is calculated within the material region at a given transition width Δ. In this case, the design domain will be separated into three phases, which can be used to represent the void, the wrapping boundary layer, and the interior infill region of the multiscale structure, as shown in Fig. 7(d).
With the gradient-based optimizer MMA, the optimization process requires the calculation of sensitivities. The  sensitivities with respect to will be derived based on different optimization problems. The sensitivities for the material properties can be analytically calculated via a forward finite difference scheme. Generally, the derivative of a given function f at Point x can be defined by the limit: After the concurrent topology optimization process, the optimal macroscale structural topology and the optimal infill material properties can be obtained. To find out the actual layout of the infill structure, a second optimization can be set up in a least square manner. An extra isotropic condition is also included as an optimization constraint to ensure the isotropy of the infill structure.

Minimize
: aðu,vÞ ¼ lðvÞ, where C H ijkl is the component of the homogenized elasticity tensor and their optimal value targets, got from the concurrent optimization, are noted as C * ijkl , HðΦÞ is the Heaviside function of the level set Φ, V const is the design volume constraint, u is the displacement field, and v is the test function. The weight of the distance regularization energy functional R is noted as w # . The bilinear energy form a(u,v) and the linear load form l(u) inside the computational domain D, with p representing the body force and τ representing the boundary traction, are described as The structural isotropy [68] is ensured when Essentially, the isotropy is ensured when the current C H 1212 hits the targeted C * 1212 calculated from Eq. (18). One more thing that needs to be mentioned is the monitoring of the isotropy during the optimization. Ideally, when this least square optimization comes to its convergence, the selected elastic tensor entries will hit the targeted values. However, in practice, this condition may not be satisfied 100%. But even without hitting all the selected targets, the isotropy may still be valid as long as the condition in Eq. (18) is hold. An isotropy polar plot index [70] is utilized in this research to visualize the isotropy of the infill structure throughout its optimization process, for an easy isotropy monitoring. The details are introduced as follows: For a unit cell with a rotation of angle , as shown in Fig. 8, the original elastic tensor will be transformed as shown in Eq. (19) and the corresponding expressions are listed in Eq. (20).
By plotting the normalized value of each expression in Eq. (20) with changing from 0°to 360°, an ideal isotropic material will generate a circle with the radius of 1. If the material is not isotropic, the plot will deviate away from the standard circle. This is illustrated in Fig. 9 where the red dashed line circle with a radius of 1 represents a plot of an isotropic material while the blue curve represents an anisotropic one. By plotting the current material with blue Fig. 8 The coordinates change of a unit cell curve to be compared with the stand red circle throughout the optimization process, the isotropy of the structure can be monitored.
Another detail that needs to be mentioned is how to determine the bounds for the infill structure material properties. The bounds should be chosen reasonably to get a valid design. However, since the infill structure only occupies a portion of the total volume of the design domain, the overall Young's modulus in each direction would be lower than the constructing material itself [74]. In the current research, the bounds are set by using Ref. [27] as a reference, where microstructures were designed to reach maximum material properties.
At last, a systematic flowchart is drawn in Fig. 10 for a more intuitive representation of the entire proposed concurrent topology optimization process.

5.1
The concurrent design of a Michell-type structure and its infill material property In this section, a two-by-one design domain is used for the designing of a multiscale Michell-type structure with fixed-fixed supports and multiple loads. The design domain is shown in Fig. 11. The forces F 1 = 1 and F 2 = 0.5 are applied at the bottom edge. The overall structure optimization statement can be made as where E ijkl ε ij ðuÞε kl ðvÞ represents the strain energy density inside the design domain. The sensitivity analysis can be Fig. 9 The isotropy polar plot Fig. 10 The flow chart for the concurrent topology optimization process Long JIANG et al. Concurrent optimization of structural topology and infill properties taken from Ref. [62], which takes the form of Eq. (22). In these equations, a smoothed Dirac delta function is represented as δðΦÞ. As can be noted, the shape derivative is extended to the whole design domain naturally. In this example, the design domain is discretized into 100Â50 quadrilateral elements. The Young's modulus for the shell material is set to be 1 and the void to be 0.001. As for the infill region, the range for the E varies from 0.01 to 0.165 and the range for the Poisson's ratio varies from 0.3 to 0.4. The minimum mean compliance optimization is carried out under the total volume constraint of 60%. The optimization convergence history, together with the level set evolution and the design evolution are shown in Fig. 12. The optimal infill material properties E is 0.165 with infill material Poisson's ratio 0.3. The overall objective function is minimized to 12.4199 and the final total volume is 0.60246.
j ¼ 1,2,:::,n: With the optimal infill material property targets, the infill structure optimization can be carried out following Eq. (16). The optimization convergence history, the evolution-history for the level set function, the design evolution and the final design polar isotropy plot are shown in Fig. 13. It can be seen from Fig. 13(a) that with the least square objective function minimized, the selected targets are hit. The isotropy of the infill structure can be verified intuitively from Fig. 13(b) where the blue-color polar plot of the current infill structure is almost overlapping with the reference red unit circle at the end of the optimization process.

Concurrent designing of an NPR structure and its infill material property
To achieve a structure with desired Poisson's ratio, a common approach is to formulate a least square problem to design a structure with a specified elastic tensor. However, from the authors' experience, trying to hit all targets simultaneously requires fine adjustment of parameters inside the optimization algorithm, which is tricky in most Fig. 11 The boundary condition for designing the Michell-type structure with fixed-fixed supports and multiple loads Fig. 12 Convergence history and the design evolution for the macroscale Michell-type structure with fixed-fixed supports and multiple loads. Upper: The level set function evolution; lower: The corresponding design evolution of the times. To decrease the number of adjustable weights, NPR structure design can be formulated directly as aðu,vÞ ¼ lðvÞ, where v t is the Poisson's ratio target. To ensure an overall orthotropic structure behavior, C 1111 and C 2222 are constrained to be equal. Considering that the homogenized elastic tensor entries can be isolated out by the aforementioned strain energy method, the design sensitivities can be borrowed from Eq. (22). In this example, a 1Â1 domain is discretized into 50Â50 quadrilateral elements. The Young's modulus for the void region is 0.001 and 1 for the wrapping shell, respectively. For the infill structure, the Poisson's ratio is given a range from 0.4 to 0.5. The range of the Young's modulus of the infill structure is from 0.01 to 0.13. The total volume is constrained within the range from 49% to 51%. The overall Poisson's ratio target is set to be -0.3. The design evolution and the convergence history are shown in Fig. 14. The final total volume ratio is 0.5098 and the final overall Poisson's ratio is -0.25841. For the infill structure, the final Young's modulus is 0.13 and Poisson's ratio is 0.4.
With the infill material properties achieved, the secondstage topology optimization is formulated in the least square manner. The design domain and the material properties for void and material are kept the same. The design evolution and the convergence history for the infill structure are shown in Fig.15.
Furthermore, the NPR performance for the entire multiscale structure is verified in Fig. 16. A 0.3 downwards displacement is added on a 4Â4 multiscale NPR structure array. The boundary conditions are shown in Fig. 16(a) and the pre/post-deformation simulation results for the multiscale structure are shown in Figs. 16(b) and 16(c), Fig. 13 Convergence history and the isotropy of the infill structure design for the Michell-type structure. (a) The convergence history and design evolution for the infill structure: The level set function evolution (upper) and the corresponding design evolution (lower); (b) the isotropic polar plot of the corresponding infill structure throughout the optimization process: Reference standard circle (red) and isotropic polar plot of the current infill structure (blue); (c) the elastic tensor of the corresponding infill structure The corresponding design evolution Fig. 15 Convergence history and the isotropy of the infill structure design for the NPR structure. (a) The convergence history and design evolution for the infill structure: The level set function evolution (upper) and the corresponding design evolution (lower); (b) the isotropic polar plot of the infill structure throughout the optimization process: Reference standard circle (red) and isotropic polar plot of the final infill structure (blue); (c) the elastic tensor of the corresponding infill structure respectively. A shrink from the right side of the structure array can be observed from Fig. 16(c), which verifies the NPR performance.
6 Multi-control-point conformal mapping, the mapped optimization results and its manufacturing With the achieved optimal overall structure and the optimal infill structure, in this section, these two parts are integrated together via the local shape-preserving conformal mapping [75]. The shape preserving effect can be visualized in Fig. 17, as the small circles on a human face will remain their circular shape after the conformal mapping.
Firstly, the concept of conformal mapping is introduced here in smooth settings. Let ω ¼ f ðzÞ : ℂ↕ ↓ ℂ be a complex function on the plane. Denote where i is the unit imaginary root. Then f is said to be conformal if ∂f ∂z ¼ 0: (25) In this research, the local shape-preserving character can ensure that the infill structure will maintain its designed material properties. This can be a great advantage compared with putting the infill structure into the design in a periodic manner and trim the ones on the design boundary. In that case, the elements cut by the design boundary might result in overhanging structure, isolated islands and some other unwanted results. Those defects will sabotage the designed infill structural properties and bring challenges to the manufacturing process.
Under discrete settings, conformal mappings can be computed by discrete Ricci flow method [76][77][78]. Some more discrete Ricci flow algorithms regarding efficiency and adaptivity improvements are reported in Refs. [79,80]. For further conformal mapping algorithms, the readers are referred to Ref. [81] for more information.
Given a triangular mesh Σ ¼ ðV ,F,EÞ, a face element is denoted with corner vertex v i , v j , and v k by f ijk , and the angle between rays are denoted as ij where ∂Σ is the boundary of mesh Σ. Now the discrete Ricci flow is defined as follows. Given a circle packing metric to Σ, i.e., on each vertex v i a positive real number is defined as γ i , then the edge length between vertices v i and v j is l ij ¼ γ i þ γ j . With those parameters, all angles can be calculated in Σ. Denote u i = logγ i , then the discrete Ricci flow is defined as where k ¼ K 1 ,K 2 ,:::,K n À Á T is the user defined target curvature.
In our case, a conformal mapping from an irregular planar region to a polygonal region needs to be found such that the inner angles are either π=2 or 3π=2. The polygonal region is filled with regular infill structures, and then they are mapped back with the inverse of the computed conformal mapping, which is also a conformal mapping that preserves local shapes. To realize this, on the boundary of input mesh ∂Σ, multiple control points W :¼ fw 1 ,w 2 ,:::,w k g ∂Σ can be selected based on the need. Then the target curvature on each vertex are defined by where K i ¼ π=2 is chosen if the target polygonal region has an outward right angle and K i ¼ -π=2 is chosen if the target polygonal region has an inward right angle at Point v i . Compared to Ricci flow method with only four control points, as shown in Fig. 18, which will eventually map the input region to a rectangular region, the proposed multicontrol-point method provides more flexibility with the benefit of lower area distortion. The proposed mapping process is illustrated in Fig. 19 using the mapping of the Mitchell-type structure as an example. Considering the symmetry, only the right half of the structure is taken into the mapping procedure and the final structure, as shown in Fig. 19(i) is achieved by a simple mirroring technique. Following the proposed mapping process, both the mapped Mitchell-type structure and the NPR structure are shown in Figs. 20 and 21, respectively. It is worth noting that after the mapping, the infill and the shell are integrated together. However, the actual material used for the final multiscale structure should be less than its designed total volume ratio since the mapped infill structure is porous. Practically, the actual volume ratio for the entire multiscale structure should be the macroscale structural volume ratio times the infill micro/meso scale structural volume ratio.
With the mapped structural design, the final manufacturing process can be carried out. The mapped result is sent to a fused deposition modeling (FDM) 3D printer and manufactured by using polylactic acid (PLA) filament. The printed result is shown in Fig. 22.

Conclusions
In this paper, a concurrent CBF-based PLSM topology optimization is proposed for the designing of multiscale structures. By using the newly constructed CBF, the explicit design variable bounds can be passed to the optimizer, which is a significant advantage over the RBFbased PLSM. Driven by MMA, the proposed approach can efficiently handle multiple constraints. Coupled with mathematical programming, the multiscale structural topology optimization can be directly carried out in a concurrent manner. The distance-regularized level set function is used to discriminate the shell and the infill region from the macro scale. From the micro/meso scale, the isotropic infill structure will bring advantages to the structural mapping/manufacturing process. With the help of the local shape-preserving conformal mapping, the infill material properties can be mathematically preserved. By introducing multiple control points, the current proposed conformal mapping can be more flexible and adaptive in Fig. 18 The conventional conformal mapping Ricci flow method with four control points Fig. 19 The multi-control-point conformal mapping process. (a) The shell-infill structure optimization result; (b) half of the original design; (c) the infill region; (d) the meshed infill region with multiple control points (red) and one central point (green); (e) the zoom-in view of the triangular mesh; (f) the isotropic infill unit cell; (g) the mapped infill structure; (h) half of the shell-infill mapped structure; (i) the final shell-infill multiscale structure Fig. 20 The conformal mapping result for the Mitchell-type structure handling complex geometries. The mapped multiscale structure can be directly sent to 3D printers for manufacturing, to conclude the design-mapping-manufacturing multiscale structure topology optimization process.