A level set shape metamorphosis with mechanical constraints for geometrically graded microstructures

Functional performance benefits brought by graded materials or a variation of microstructures have been established. However, the challenge of connecting the adjacent microstructures remains open. Without connected microstructures, the design is only theoretical and cannot be realized in practice. This paper presents a shape metamorphosis approach to address this challenge. The approach formulates an optimization problem to minimize the difference between two shapes, where the shapes generated during the optimization iterations form the graded microstructures. A stiffness-based constraint is proposed to eliminate potential breakage or shrinkages in generating the intermediate shapes, which are used to build a transition zone for the two shapes to be connected. This also ensures that the optimal functionality is not compromised significantly. A series of numerical studies demonstrate that the proposed approach is capable of producing smoothly connected microstructures for multiscale optimization.


Introduction
The rapid advance of additive manufacturing (AM) offers freedom to engineering design owing to its capability of the unprecedented level of complexities, length scales, and material compositions (Rosen 2007) thereby opening up a new design freedom at miscroscale that was not easily accessible with traditional manufacturing technologies. The AM processes are classified into seven categories: VAT photopolymerization, material jetting, binder jetting, material extrusion, powder bed fusion, sheet lamination, and directed energy deposition ( variety of AM techniques have been used to manufacture structures made of architected cellular materials, which are usually computationally designed and optimized. For instance, unit cells of polymer have been manufactured by stereolithography in the category of VAT photopymerization (Tumbleston et al. 2015); the selective laser melting in the category of powder bed fusion has been used to manufacture cellular materials based on metal (Cansizoglu et al. 2008); electrophoretic deposition has been used to make ceramicbased cellular materials (Jang et al. 2013). With the feasibility of manufacturing small length-scale features in lattice and cellular microstructures demonstrated, interests are growing in designing such multiscale features for novel metamaterial and multifunctionalities.
Topology optimization has been widely accepted to be an ideal design method for additive manufacturing (Liu et al. 2018). It offers a systematic and mathematically driven framework for exploring these newly accessible design spaces and enables integrating manufacturing constraints into the design stage. In the past three decades, various topology optimization approaches have been applied to design material microstructures to achieve specific properties. In general, the material design is achieved by solving an optimization problem based on the homogenization theory. Comprehensive reviews on this topic can be found in, e.g., Schaedler and Carter (2016) and Osanov and Guest (2016). Material microstructures have been designed to achieve (a) extremal properties, e.g., maximal bulk or shear moduli (Huang et al. 2011;Neves et al. 2000), negative Poisson's ratio (Sigmund 2000), extremal thermal expansion (Sigmund and Torquato 1997); (b) targeted properties, e.g., effective Young's moduli (Yang et al. 2013), Poisson's ratio (Sigmund 1995); (c) multifunctionality, e.g., simultaneous transportation of heat and electricity (Torquato et al. 2002), and stiffness and conductivity (Challis et al. 2008).
To further exploit the enhancement of engineering performance brought by tailoring material microstructure, spatially varying material microstructures have been investigated. Multiscale topology optimization is one of the strategies getting increasingly explored in recent years to simultaneously design material distributions in both structural-and material-scale domains. Rodrigues et al. (2002) proposed a hierarchical design to simultaneously solve two subproblems related to the structure and material designs using the density-based approach. Recently, several multiscale optimization schemes have been proposed using various optimization techniques, e.g., the level set method (Wang et al. 2016;Sivapuram et al. 2016), the bidirectional evolutionary procedure (Yan et al. 2014), the deshomogenization method pioneered in Pantz and Trabelsi (2008), then elaborated in Allaire et al. (2018) and Groen and Sigmund (2018), and their applications range from the classic compliance minimization problem to thermoelastic (Deng et al. 2013;Yan et al. 2016), natural frequency , and frequency response (Vicente et al. 2016). This technique is increasingly receiving attention in recent years to design structures with optimal distribution of graded microstructures in both micro-and macroscales to maximize the benefits Wang et al. 2017).
One consequence of using the homogenization theory is that the connectivity between adjacent microstructures may not be guaranteed due to the assumption of the periodic boundary condition. This means the spatially varying microstructures may not actually be connected, resulting in fictitious configuration that cannot be held together and the discontinuity at the interfaces can give rise to nonoptimal functionalities. They are impossible to fabricate and do not possess the desired functionalities (Vermaak et al. 2014;Liu and Ma 2016). Recent literature aims to address this connectivity issue by designing functionally graded cellular materials (FGCMs). FGCMs are to fulfill a specific property or functionality, and the connectivity is ensured by grading the microstructures rather than changing the topologies of the microstructures. Recently, Zhu et al. (2019) extended the classical asymptotic homogenization theory to design graded infill structures with the use of the coordinate pertubation technique (Liu et al. 2017). As the topology does not fundamentally vary in the structural domain, the graded microstructures are naturally connected. Zhou and Li (2008) proposed a group of three methods, including (a) kinematical connective constraint, (b) pseudo load, and (c) unified formulation with non-linear diffusion, in order to preserve the connectivity between adjacent microstructures in the optimal design of twophase FGCM with a gradual change in targeted elastic properties. The first two approaches require prescribing connective areas, while the third approach ensures the connectivity implicitly. Given the freedom of the third method, it was further developed for 3D cases in Radman et al. (2013). Recently, the first approach was adopted by  to design gradient cellular structures with multifunctional considerations for microstructures. Panetta et al. (2015) use a family of tileable parameterized 3D small-scale structures to vary material properties across an object to achieve graded material structures. An innovative scheme was developed recently by Du et al. (2018) through imposing a connectivity index constraint in the standard multiscale topology optimization problem. However, it can be expensive in terms of computational cost to satisfy additional constraints especially when there are many microstructures.
Another strategy is to connect microstructures through post-processing. Distance field-based shape interpolation approach is one of the popularly used techniques to manipulate the geometries via linear or power interpolation techniques (de Buhan et al. 2016;Cramer et al. 2016;Wang et al. 2017;Dapogny et al. 2017). The interpolated geometries are then used to build a transition zone to connect the two distinct microstructures as shown in Fig. 1. This shape interpolation technique has been developed originally for morphing between images in computer graphics community (Raya and Udupa 1990). Comparing with the connectivity constraint-based approaches, the postprocessing approaches are conceptually simple and can be more robust, although they may result in a degradation of the structural performance of the design since the postprocessing stage is not taken into account during the optimization. Hence, the post-processing-based techniques are the focus of the present study. One challenge identified in the interpolation approaches is that the connectivity may not be guaranteed, if the two shapes to be connected are significantly different from each other. Figure 2 illustrates this via an example where transition microstructures are constructed by shape interpolation to connect two microstructures in Fig. 2a and c. The gray area represents void, while the black area is solid material. Three intermediate microstructures from the shape-based linear interpolation method are given in Fig. 2b. It shows that the intermediate shapes do not extend to the edges of the unit cell. Therefore, these topologies do not connect. A closer examination of the two topologies in Fig. 2a and c reveals that the outer edges of the topology in Fig. 2c contract to fit the corresponding regions of Fig. 2a during the morphing process. In addition, the second intermediate shape in Fig. 2 has four disconnected parts. To avoid these issues, Wang et al. (2017) conducted the shape interpolations between a fully solid material and a cellular material to generate a family of intermediate shapes, and microstructures without breakages were selected to build the transition zone. This breakage and shrinkage phenomenon was not reported in Cramer et al. (2016) mainly due to the fact that the interpolation was only applied to microstructures with different volume fractions but the same topology. In multiscale topology optimization, distinct cellular structures are expected to be connected, and more sophisticated shape interpolation approaches are needed.
In computer graphics, shape interpolation belongs to one of the intensively studied topics, namely shape morphing, which is a process where a shape continuously changes from one into another. Various shape morphing methods have been developed for 2D images and 3D objects. For image morphing, interpolation approaches based on a scene or an object are the two main categories. Scene-based approaches use the image intensity or feature to conduct linear or non-linear interpolation between intensity values that are correlated on their matrix locations within respective images (Grevera and Udupa 1996;Frakes et al. 2008;Penney et al. 2004). Object-based approaches operate on image features or contours (Beier and Neely 1992; Lee and Lin 2002;Raya and Udupa 1990). Depending on the way to represent objects, 3D morphing techniques are broadly divided into boundary or surface representation approaches and volumetric methods. Surface-based methods (Shlafman et al. 2002;Kanai et al. 2000) continuously change one polygonal surface into another, while the volumertic approaches (Breen and Whitaker 2001;Kilian et al. 2007;von Tycowicz et al. 2015) modify the voxel values of a volume data set in order to smoothly transform the object defined by the volume from a source shape into a target shape. Although several methods exist, they are not directly applicable to generate a set of intermediate shapes with the preservation of connectivity for adjacent shapes. This is because morphing has been used mainly in visualization for visual effects in movies and animation. Therefore, they do not concern about the spatial connections between adjacent shapes. Furthermore, many of these methods rely on a substantial human interaction. Since the aim of our method is not to make it look visually natural to human eyes but for physics-based functionalities, such a level of reliance on human interaction can be challenging, unlike for the purpose of visual effects. This paper, therefore, adopts a morphing method as a basis and builds on it to formulate a tailored post-processing method for multiscale topology optimization. Among the morphing methods, a shape metamorphosis method developed by Breen and Whitaker (2001) is found to be suitable for our purpose. It formulates the metamorphosis process as an optimization problem to maximize the similarity between the intermediate object and the target object iteratively. By doing so, it enables to smoothly deform the source object into the target object. The objects are represented by signed distance fields, and the optimization problem is solved by the level set method. This method thus naturally integrates with our level set topology optimization.
Inspired by level set shape metamorphosis, this paper introduces a new computational method for generating a series of microstructures that smoothly connect two adjacent microstructures for multiscale topology optimization.
The key idea here is to eliminate the shrinkage and breakage problems in the intermediate geometries by imposing a mechanical constraint to naturally produce an engineering design with mechanical integrity. The detailed formulations of the new method are introduced in Section 2, and its performance is evaluated through numerical studies in Section 3.

Level set shape metamorphosis with constraint on mechanical property
Mathematically, a shape or an object can be implicitly represented by a level set function, φ, as: where x is an arbitrary point in the domain Ω. Signed distance is a common choice for defining φ in level set structural topology optimization and shape metamorphosis.
Consider two objects, Ω A and Ω B and their boundaries denoted as A = ∂Ω A and B = ∂Ω B , respectively. Using the level set function, the similarity of these two objects or the extent of Ω A overlapping with Ω B can be measured by where H (·) is the Heaviside function, and φ B is the level set function of the object B. This integral achieves its maximum if the two objects are the same or completely overlaps. (Breen and Whitaker 2001) used this property to develop a shape metamorphosis method to transform an object to a target one Ω T , which is formulated as an optimization problem where the similarity or the overlap area of two objects is maximized, that is expressed as: where Ω i is the intermediate shape, φ is the level set function for Ω i , and φ T is the level set function for the target shape. However, the intermediate shapes may suffer from the shrinkage and breakage problems if the two objects to be transformed are distinct from each other. In order to prevent these problems, a mechanical constraint is imposed on (3). This enforces a requirement for the intermediate shapes generated in the optimization process to transfer the mechanical loads between the boundary conditions, hence the intermediate shapes are well-connected. In addition, this optimization problem will be solved by the level set methods (Breen and Whitaker 2001;Osher and Fedkiw 2003) that makes it straightforward to be embedded into the level set structural topology optimization (Dunning and Kim 2015;Sivapuram et al. 2016). The constrained shape metamorphosis is thus formulated as follows: where g ≥ g * denotes the stiffness-based constraint with g * denoting the lower bound.
As the objective of the proposed method is to address the microstructural connectivity issue in the multiscale topology optimization, the objects to be connected are thus microstructures from the inverse homogenization method, which uses the homogenization method to link the underlying material microstructure and constituents with the effective elastic properties to design the material topology (e.g., Sigmund and Torquato (1997)).
The effective elastic tensor for a microstructure, usually chosen as a periodic unit cell, Y , can be calculated as: where |Y | is the volume of the cell in a periodic domain Y , D pqrs is the elastic tensor of the base material, d is the spatial dimension,ε pq is the prescribed strain field, and ε (ij ) pq is the corresponding strain field induced byε pq and is solved by the the finite element model with the periodic boundary conditions at microscale: where v is the virtual displacement field belonging to the kinematically admissible displacement space χ. We refer the readers to Guedes and Kikuchi (1990) and Perić et al. (2011) for more details on the numerical implementation of (6).
Here, the inverse homogenization method was implemented in the level set topology optimization context. It combines the homogenization method to obtain the effective properties with the level set method to achieve shape evolutions and topological changes of the microstructure until the desired material properties are obtained. The readers are referred to, e.g., Sivapuram et al. (2016) and Du et al. (2018) for more details on the implementation in the level set method. A brief description of the inverse homogenization method to design the material microstructure is given as follows for the completeness. To design the material for a given property as a function of the effective elastic tensor D H , it is formulated as an optimization problem to properly distribute the solid material within a unit cell Y subject to a given material volume, α|Y |. It can be mathematically stated as follows: The stiffness-based constraint is based on the mechanical properties, including bulk modulus, shear modulus, and strain energy, that are derived from the effective elastic tensor calculated by homogenization methods. Without loss of generality, three types of constraints are introduced and investigated in the 2D context. The constitutive law between macroscale stress,σ , and strain,ε, can be characterized by the generalized Hooke's Law as: With this at hand, the constraints are established as: (i) Bulk modulus, κ: (ii) Shear modulus, μ: (iii) Total strain energy, : whereε i is the ith unit strain, andσ i is the corresponding stress.

Sensitivity analysis
The constrained shape metamorphosis problem of (4) is solved with the level set topology optimization method, which solves the Hamilton-Jacobian equation to move the boundary of a shape (Sethian and Wiegmann 2000;Wang et al. 2003;Allaire et al. 2004). It is expressed as the following advection equation: with V n denoting the normal inward velocity of structural boundary and t a pseudo time. It requires the firstorder derivatives of the objective function and of the constraints with respect to the design variables to form the velocity function, V n , to ensure that the structure evolves toward a maximum of J . The contribution of each shape sensitivity term to the velocity is determined by solving the optimization problem. The shape sensitivity for the objective function, which measures the similarity between the two microstructures to be connected in terms of their overlapping areas, is straightforwardly expressed as: For the constraints, they are functions of the effective elastic tensor D H of a microstructure. Hence, the essence is to calculate the shape derivative of D H with respect to Ω. The shape sensitivity of D H is well-known and written as (see, e.g., Sivapuram et al. 2016): With the shape sensitivity of D H at hand, shape sensitivities of the proposed constraints are straightforward to calculate through simple operations. For instance, the shape sensitivity of the bulk modulus-based constraint in (9) is the sum of corresponding four elements of ∂D H ∂Ω in (14).

Generation of intermediate microstructures in the transition zone
Every microstructure from the morphing process can be the candidate to build the transition zone as shown in Fig. 1, but it may not be necessary to use all of them. In order to select key microstructures to form the transition zone, we develop an algorithm to measure the connectivity at the interface boundary between each pair of adjacent microstructures, which are denoted as the left-and the righthand-side cells or microstructures as shown in Fig. 3. For this, the connectivity of two microstructures at the interface is quantitatively evaluated by an index, β, defined as: where φ(x) is the implicit function to describe the microstructure, e.g., signed distance function, and the subscripts l and r refer to the right-hand side and the lefthand side cells, respectively. If two microstructures are fully connected, β = 1.0, while β = 0 indicates completely disconnected microstructures (see Fig. 3). It should be noted that the denominator of (15) will always be positive in the common design for the minimum compliance. Based on this connectivity index, a selection scheme, which aims at using the least number of shapes to build a transition zone while maintaining a desirable level of connectivity, is developed to draw the intermediate microstructures. It includes the following steps: (1) Start the first microstructure generated from the metamorphosis; (2) put the first microstructure and the second microstructure as a pair and calculate their connectivity index; (3) if the connectivity index is larger than a given threshold of connectivity index, β 0 , then form the first microstructure with the next one and calculate their connectivity index; (4) otherwise, the first microstructure should connect with the last one, e.g., the ith microstructure, with the connectivity index over β 0 . (5) The procedure is then repeated to find the optimal shape to connect the ith shape from the rest. (6) The selection algorithm stops until the last shape is connected. (7) Finally, the selected shapes are used to build a transition zone between the two microstructures to be connected. It should be noted that the number of selected microstructures may vary depending on the imposed connectivity index. We will discuss this point in the numerical studies.

Numerical implementation
The implementation of the above method involves two parts, including finite element analysis (FEA) and level set method (LSM). In structural topology optimization, the area-fraction weighted fixed grid approach (AFG) has been employed in FEA and LSM. A design domain is discretized into a finite element model with four nodes bilinear elements, and a structure is expressed by the area fraction of material in an element, where the non-structural regions within the fixed mesh are filled with a fictitious weak material. As previously mentioned, a structure is described by a level set function. Hence, the design domain is also discretized into a level set mesh with bilinear elements for conducting the numerical analysis for the level set method. The resolution of the finite element model can be different from the level set mesh. Here, the same resolution is applied to the two meshes. In a level set domain, each node will have a level set value. The level set element is cut depending on the level set values, and the area fraction for the element can be calculated to link with the finite element model for conducting the structural analysis. Hence, (4) is now ready to be implemented in a standard level set optimization method. Sensitivities are determined numerically from the stress and strain fields using the fixed grid FEM, and the optimization problem is solved to find the boundary velocities to evolve the structure using (12).

Connecting microstructures from inverse homogenization
This section presents numerical investigations of the proposed method on microstructures designed for specific mechanical properties. The focus of the investigations are in two folds, i.e., the preservation of mechanical functionality and the effectiveness of the transition zone connecting two microstructures. The microstructures selected to be connected have distinct topologies which suffer from the aforementioned breakage or shrinkage issues when using the standard shape interpolation techniques. Two sets of microstructures are designed to demonstrate the robustness of the proposed approach. In the first set, denoted as Experiment 1, the two distinct microstructures to be connected are generated from the inverse homogenization method for the maximum bulk modulus and the maximum shear modulus as shown in Fig. 4a and b for a given material usage in terms of area fractions of 25% and 40%, respectively. It should be noted that they are actually local optimal solutions for the desired properties as the problems are solved by the gradient-based scheme. The base material is isotropic with Young's modulus E = 1.0 and Poisson's ratio ν = 0.3. They are found through solving an inverse homogenization problem using the level The microstructure in Fig. 4a has the bulk and shear moduli of κ = 0.0753 and μ = 0.0028, while those properties for the microstructure in Fig. 4b are κ = 0.1332 and μ = 0.1093. In the second set, denoted as Experiment 2, microstructures with maximum values on shear modulus and D H 1112 , respectively, are to be connected. Figure 4c gives the optimum microstructure with the maximum value on D H 1112 under an area fraction of 25% having κ = 0.0632 and μ = 0.0627.

Investigation of transition microstructures
For the first set of microstructures in Experiment 1, three formulations are studied. Two constraint formulations based on shear and bulk moduli are naturally constructed considering the mechanical functionalities of two microstructures. In addition, the constraint based on strain energy in (11) is also considered as it is a more generic and global measure. The results from each of these constraints together with the one without imposing a constraint are compared. Since the number of iterations required for the metamorphosis subject to different constraint varies, the history curves are presented with a normalized number of iterations in the abscissa axis. Figure 5a shows the history of the objective function value, i.e., the similarity between the current and final microstructures, and Fig. 5b the area fraction of the microstructure at each iteration when morphing from Fig. 4b to a. It can be seen from Fig. 5a that the objective function, J , monotonically increases, which indicates that the overlapping area between the current and target microstructure evolves smoothly. It can also be read from the morphing process given in Fig. 5b that the intermediate microstructures start from the source one and gradually gain materials to form microstructures that are the union of the two microstructures to be connected, and the redundant materials comparing with the starting one are then gradually removed. This trend is consistently seen in all the cases subjected to different constraints as given in Table 1. Now, let us examine how mechanical properties evolve in the morphing process. Figure 6a presents the history  Fig. 6b.
To further investigate the proposed method, the second set of microstructures in Experiment 2 is examined. As the two microstructures to be connected have the maximum properties on D H 1112 and μ, respectively, two constraints based on these two properties are imposed on the metamorphosis. Again, another metamorphosis without a constraint on the mechanical property is given as a reference. Figure 7 shows the history of the similarity function between the generated microstructure and the one to be connected and the evolution of the intermediate microstructures in terms of their area fractions. Similar evolution processes as in the first set of microstructures can be observed. Among the three constraint-based formulations, the strain energy-based one shows slightly a different history in the change of area fractions of the intermediate microstructures in Fig. 7b. They gain more material before removing the redundant regions, which implies the strain energy-based constraint may preserve the structural integrity better. This is reflected in Fig. 8 where both moduli are maintained high when the strain energy constraints are used. The differences between with and without a constraint on the mechanical property are evident also. Both bulk and shear moduli for the intermediate microstructures immediately drop to near zero when a constraint on the mechanical property is not considered, while these two elastic properties are around or above the lower values when a mechanical constraint is considered.

Forming a transition zone with the selection scheme
The proposed selection scheme is used to select shapes to form a transition zone by using a targeted connectivity index, β 0 . To illustrate the proposed selection scheme, two key factors including the targeted connectivity index and the applied constraint are investigated. The construction of a transition zone for the microstructures with maximum shear modulus and maximum bulk modulus in Experiment 1 is investigated first. To illustrate the influence of the target connectivity index on the transition zone, four target connectivity indices, β 0 , of 0.85, 0.8, 0.75, and 0.7 have been considered to illustrate a decreasing extent of connection in selecting shapes. Table 1 lists the built transition zone for the four β 0 values. It can be seen that the number of key intermediate microstructures required to form a transition zone reduces with the decreasing target connection index. For instance, 32 microstructures have been selected to form the transition zone for a targeted index of 0.85 in Table 1 with the use of the bulk constraint, while only 14 microstructures are selected for the target connection index of 0.7. All the three cases of constraints in Table 1 show that the pattern of transition zone does not have significant changes with the change in the target connectivity index. It demonstrates that the metamorphosis process with mechanical constraints reliably produces a set of geometrically smooth and connected microstructures. By comparing the transition zones under the same target connectivity index in Table 1, the effect of the imposed constraint on the formation of transition zones can be observed. The transition zones with the target connectivity index of 0.85 generated from the three constraints are compared. Only 16 microstructures have been selected for the shear modulus constraint. The bulk modulus constraint performs similarly to the compliance constraint in terms of the number of microstructures used to build the transition zone. A total of 36 microstructures are required to build the transition zone when using compliance constraint, while the constraint on bulk modulus requires 32. The constraints on bulk modulus or compliance are actually against the morphing process from the microstructure with the maximum bulk modulus to the microstructure with the maximum shear modulus, and they lead to the fact that more microstructures are required to build the transition zone. Although the constraint on compliance is outperformed by the other two constraints in this study, this is likely to be problem dependent.
Similar studies are carried out for the second set of microstructures in Experiment 2, which have maximum values on shear modulus and D H 1112 as shown in Fig. 4c. The influences of the target index and the type of constraint on the formation of a transition zone are given in Table 2. Among the three constraints, the constraint on shear modulus requires the least number of microstructures to form the transition zone, while the constraint on D H 1112 needs more microstructures than the other two. It can be seen that the selected microstructures grade from the shape with the maximum property on D H 1112 and the other with the maximum property on μ, where the shear resistance gradually improves with the growth of diagonal components. The constraint on shear modulus is more efficient due to the morphing process is toward the microstructure with dominant shear resistance, while the constraint based on D H 1112 is opposite to the morphing process.

Applications in multiscale design
We apply the proposed method to solutions from multiscale topology optimization. Two examples, bi-clamped and Lshaped beams, are presented here. In the bi-clamped beam example, the distribution of material is only considered at the microscale, while both micro-and macroscale topologies are designed for the L-shape beam.

Bi-clamped beam with topological design microstructures
A two-end clamped beam with a length of 2 m and a height of 1 m as shown in Fig. 9 is studied. Two concentrated forces of 1 N located at 0.7 m away from the two vertical edges are applied. The base material is linear isotropic with Young's modulus of 1 and Poisson's ratio of 0.3. The optimization problem is to minimize the overall compliance under a volume constraint of 40% on the microscale design domain; only material microstructures are thus optimized to achieve the objective by using the multiscale optimization approach introduced by Sivapuram et al. (2016). The optimization problem is stated in (16). The design domain is discretized into a mesh of 200 × 100 quadrilateral elements and divided into five regions as shown in Fig. 9. A unit cell is assigned to each region with a mesh of 100 × 100 linear quadrilateral elements. Due to the symmetry of boundary and loading conditions, three microstructures are actually optimized as Minimize : It can be seen from Fig. 10 that the solutions are not well connected at the interface. According to (15), the connectivities are 0.4 and 0.2 for interface between microstructures 1 and 2, the one between 2 and 3, respectively. The proposed approach is thus applied to build connections for them. Three constraints are selected in consideration of the underlying mechanical behaviors of the microstructures. The solution for regions 1 and 5 featured prominent diagonal struts for sustaining shear stress. Thick vertical components in microstructure for regions 2 and 4 and thick horizontal components in the solution for microstructure 3 show that the microstructures are governed by bending stresses. Hence, the shear and bulk modulus-based constraints together with the compliance constraint are considered. A set of intermediate shapes is generated when applying the proposed metamorphosis approach with one of the three constraints. Figure 11 shows the histories of bulk and shear moduli for the intermediate shapes that are candidates to connect microstructures 1 and 2. It shows that all of the three constraints are effective in generating the intermediate shapes without losing the mechanical functionalities. Similar phenomena have been observed for microstructures 2 and 3.
These intermediate shapes are then used to build transition zones for the two interfaces through the selection scheme based on the connectivity index from (15). Here, a targeted connectivity index of 0.85 is applied, and the resulting transition zones are shown in Figs. 12 and 13. For the interfaces between microstructures 1 and 2 or 4 and 5, all of the constraints can generate well-connected transition zones as shown in Fig. 12. The constraint based on shear modulus in (10) leads to thicker diagonal components to satisfy the constraint on shear modulus. In contrast to the transition zone based on the bulk modulus, the constraint based on compliance in (11) generates the greatest number of intermediate shapes. Similar trends are observed in the transition zones for connecting microstructures 2 and 3 or 3 and 4 in Fig. 13. This is expected as the compliance is an overall measure.

L-shaped beam with simultaneously designed material and structure topologies
In this example, both material and structural topologies are optimized simultaneously for an L-shape beam with the configuration shown in Fig. 14, and three regions are specified, where each region can take a different microstructure. The beam is fixed at the top edge, and a point load is applied at the tip as shown. The objective function is to minimize the overall compliance under linear elasticity. The optimization problem is subjected to a global volume constraint of 18% with 40% distributed at the macroscale design domain and 45% at the microscale. It is formulated in (17) Using the formulation proposed in Sivapuram et al. (2016), design solutions are found for the initial designs given in Fig. 15a and Fig. 16. As shown in Fig. 16, the optimal microstructure for region 1 does not connect well with the one for region 2 with the connectivity index β = 0.46, and the connection between microstructures for regions 2 and 3 is reasonable with β = 0.7.
A set of intermediate shapes is generated by applying the shape metamorphosis to each pair of the microstructures. We employ the bulk and shear modulus-based constraints. The histories of bulk and shear moduli are plotted in Fig. 17 to examine the mechanical functions for the intermediate shapes. It demonstrates that all of the three constraints can successfully morph the shapes preserving mechanical functionalities. The intermediate shapes generated from  the bulk modulus and the compliance constraints do not have significant discrepancies with microstructures 1 and 2, while the bulk moduli of some intermediate shapes generated from the shear modulus constraint are around a half of the value for microstructure 1. The applied constraint helps to preserve the shear behavior of an intermediate microstructure, but it scarifies the bulk modulus. For microstructures 2 and 3, similar observations have been obtained.
The transition zones built between microstructures 1 and 2 are shown in Fig. 18 under a target connectivity index of 0.85. All of the three investigated constraints can produce effective transition zones to connect the two microstructures. The two constraints based on the bulk modulus and the compliance perform similarly as indicated in the constructed transition zones in terms of the number of microstructures, while the transition zone generated from the constraint based on shear modulus is slightly different. The selected microstructures have the desired mechanical features, where microstructures grade from tension or compression dominant structures to those with strong shear resistance. For the interface of the microstructures 2 and 3, all of the three constraints produce transition zones with smoothly graded microstructures as shown in Fig. 19. To illustrate, the L-shape beam made of these optimized microstructures is depicted in Fig. 20 together with the transition zones built from graded microstructures generated by using the bulk modulus constraint. These numerical examples therefore demonstrate that the proposed Fig. 20 An illustration of the L-shape beam with the optimized microstructures approach can build a transition zone to connect two microstructures with smoothly grading geometries and mechanical functionality.

Conclusions
A novel shape metamorphosis formulation is presented to generate geometrically graded microstructures to connect the microstructures designed by multiscale topology optimization. It is based on maximizing the similarity between the two microstructures and on satisfying a constraint based on mechanical properties to ensure the generated microstructures maintain their load-carrying capacities. This also addresses the shrinkage and breakage of geometries previously observed in other shape metamorphosis methods. This approach naturally integrates with level set-based topology optimization. In addition, a selection scheme is developed to build a transition zone with a desired connectivity level.
A series of numerical studies demonstrate the feasibility, efficiency, and effectiveness of the proposed approach. In the first set of the numerical studies, the proposed method is studied to morph between the microstructures with the distinct extremal properties. The second set of the numerical studies shows the applicability of the proposed method in the context of multiscale topology optimization. These examples demonstrate that the proposed method can be useful in designing functionally graded materials with architected microstructures. The proposed method provides a general method for connecting multiple microstructures and the choice of the constraint depends on the underlying mechanical behaviors of the microstructures. Nevertheless, further numerical investigations will continue to study the robustness of the algorithm.

Replication of results
The presented results are produced by subroutines using our in-house code Multiscale and Multiphysics Design Optimization (M2DO). The code and data for producing the presented results will be made available by request.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.