1 Introduction
Infill graded microstructural (IGM) configurations have shown their vast engineering potentials in achieving exceptional material/structural properties Lakes_Nature1993 ; DongHW_JMPS2017 ; Sigmund_APL1996 ; Kushwaha_PRL1993 ; LiuC_PRAppl2015 ; Aage_Nature2017 . For IGM design, early treatments encompass the use of asymptotic homogenisation (AH) approaches (Bensoussan1978, )
. The basic idea is to treat a spatially periodic configuration as a continuum, whose elasticity tensor at every point is calculated through solving a single set of partial differential equations (PDEs) defined over one constituting cell. The homogenised formulation was then adopted for spatiallyvarying infill configurations, where every macroscopic pixel is associated with an individual cell
(Bendsoe_CMAME1988, ; Bendsoe_StructOpt1989, ; Sigmund_IJSS1994, ). The method was further extended to other cases, such as concurrent multiscale design Rodrigues_SMO2002 ; LiuL_CompStruct2008 , multifunctional design (NiuB_SMO2009, ; DengJD_SMO2013, ; YanJ_CompMech2016, ) and problems concerning load uncertainties DengJD_SMO2017 . However, several limiting issues arise against the treatments based on conventional AH formulation. For instance, the smoothness across neighbouring cells may not be guaranteed. Moreover, the asymptotic analytical results may become inaccurate, since the presumption of microstructural periodicity no longer holds.For improvement, various types of methods for IGM representation and stress analysis have been proposed. An intuitive way of thinking is to pile up buildingblock structures, whose topology description functions (TDFs) originate from a single set of continuous functions. This can be achieved either explicitly or implicitly. The explicit way involves cell configurations composed of parameterised members (ZhouSW_JMaterSci2008, ; Radman_JMaterSci2013, ; Radman_CompMaterSci2014, ; ChengL_CMAME2019, ; LiQH_CompStruct2019, ). Thus an IGM configuration is generated by letting the corresponding parameters vary gradually (but discretely) in space. The implicit treatment is usually achieved through the introduction of a levelset function, whose multiple (but discrete) levelset contours are used to characterise the boundaries of a family of prototype cells WangYQ_CMAME2017 ; Gao_AES2018 ; ZhangY_CMS2018 ; ZhangY_SMO2019 . By either mean, graded change in microstructure should be achieved, but several challenging issues still persist. First, slight nonsmoothness across the boundaries of neighbouring cells is still an unresolved issue, since only discrete values of the controlling functions are concerned. Second, the variation of the microstructural infill patterns is only permitted to take place along directions parallel to cell edges. This greatly squeezes the resulting IGM variety.
Recently, a number of novel approaches underpinning fast IGM design are proposed, where smoothness connections within IGM configurations can be fully guaranteed. Interestingly, the restriction on oriented rectangle/cuboidshape cells is also removed in these recently proposed approaches. Besides simulations directly conducted on fine scales Alexandersen_CMAME2015 ; LiuC_JAM2017
, the novelly introduced treatments can be roughly classified into three categories: 1) conformalmappingbased approaches; 2) the dehomogenisation approaches; 3) the modified asymptotic homogenisation topology optimisation (AHTO plus) approaches.
The conformalmappingbased approaches Vogiatzis_CMAME2018 mainly concern decorating microstructural configurations on a given surface in threedimensional space. In this scenario, conforming mapping is used to ensure that the local orthogonality in microstructure is preserved as a planar microstructure is deployed on the prescribed surface. As a cost, the Ricci flow equation, a useful differential equation defined on manifolds, needs to be numerically solved to a steady state.
The dehomogenisation methods Groen_IJNME2018 ; Groen_CMAME2019 ; Allaire_CompMathAppl2019 ; WuJ_IEEE2019 ; GeoffroyDonders_JCP2020 focus on cells of laminate configurations, which theoretically guarantee the solution optimality for most cases in compliance optimisation. They first determine an optimised distribution of the laminate volume fraction and orientation (on a coarser grid). Then the actual microstructure is resolved by solving another optimisation problem projecting the laminate configurations in alignment with the desired directions of the local principal stresses Pantz_SIAMJCO2008 .
The AHTO plus approaches ZhuYC_JMPS2019 describe an IGM configuration with function composition. The generated IGM configuration can always be mapped to a spatially periodic configuration (in a fictitious space) by means of a continuous mapping function. Then the TDF of the IGM configuration is obtained through the function composition of the mapping function with the TDF describing the cell configuration (of one period in the fictitious space). Under this setting, asymptotic analysis is carried out for deriving a homogenised formulation governing the mechanical behaviour of the resulting IGM configuration. Compared with the abovementioned two methods, no postprocessing procedure, such as solving the Ricci flow equation or conducting projecting operations, is needed in the AHTO plus framework. Nevertheless, other challenging issues arise over the following three aspects. Firstly, since the constituting cells (of parallelogram/parallelepiped shape) originate from a single generating cell, thus they must be selfsimilar to each other. This dramatically reduces the IGM variety in concern. For instance, spatial change in microstructural volume fraction is still disabled. Secondly, the number of design variables, which parameterise the controlling function, is relatively small under the AHTO plus framework. Thus the identification of an appropriate set of design variables becomes a critical issue, which has yet been properly addressed. Finally, the derived homogenisation formulation involves solving various PDE systems at individual macroscopic points, generating a huge computational burden. Hence finding an efficient way of solving these PDEs XueDC_arxiv2019 is also a key factor affecting the applicability of the method for IGM design.
The present article is aimed for addressing the first aforementioned issue, that is, to essentially improve the IGM variety under the AHTO plus framework. Compared to the original treatment in ZhuYC_JMPS2019 where local microstructures have to stem from a single generating cell, here we consider a continuous menu of generating cell configurations defined in a fictitious space, which can be systematically identified by the continuous level sets of a single function, say, , termed as a “structural generating function” in this article. Then a macroscopic continuous function defined in the actual space, say , is introduced to indicate which level set it takes around point . Here the function is termed as a “levelset indicator” function. Incorporating with the original function composition treatment, a TDF associated with an IGM configuration is thus introduced, and the “graded microstructural” behaviour is carried out by the combinative effects of the following two sets of operations: manipulation of a single cell and transition among cells. The action of cell manipulation involves 1) cell rescaling; 2) cell rotation; 3) cell distortion in an isovolumetric manner. Mathematically, these operations over a single cell can be carried out by manipulating the local Jacobian matrix, which consists of the spatial derivatives of the mapping function. Such an idea for IGM representation originates from the use of levelset contours for modelling dislocation loops in crystalline materials XiangY_JMPS2009 ; ZhuYC_JMPS2015 ; Chapman_SIAP2016 . It bears certain similarity with the treatment by WangYQ_CMAME2017 , but essential differences exist. For instance, level sets of discrete values are used in WangYQ_CMAME2017 , but continuous level sets are considered here. Moreover, local rotation and distortion of the generating cells are not taken into account by WangYQ_CMAME2017 , resulting in pixelwise microstructures only. This is, however, no longer a limiting issue if the present scheme is implemented. In fact, we will show that the TDF form given by Eq. (4) effectively offers a platform for unifying most existing treatments for IGM generation, and detailed discussion is carried out in Sec. 3.4. The present algorithm also sees its generalised applications in two important aspects. Firstly, it offers a simple treatment for continuously patching arbitrarily given cell configurations, which are termed as “seminal cell” configurations in the present article. Secondly, it breeds a general algorithm for generating an IGM surface.
The IGM configurations obtained here inherit all the attractive features from that generated based on the original AHTO plus framework, such as microstructural smoothness and simple forms in TDFs. Moreover, the corresponding asymptotic formulation for stress analysis is also available, and its accuracy is further shown against finescale simulation results. Under the present framework, one can monitor an IGM configuration by varying the mapping function , the levelset indicator function and the structural generating function . Mathematically, this can be fulfilled by a topology optimisation formulation with the associated sensitivity analysis results given. All these should pave way for adapting the AHTO plus framework for fast digital design of IGM configurations.
The rest of the article is arranged as follows. In Sec. 2, the original AHTO plus method is briefly revisited with the challenging issues identified. In Sec. 3, our algorithm for IGM generation from a continuous menu of cell configurations is introduced, and two of its generalised applications are discussed: 1) IGM generation continuously patching different seminal cells; 2) IGM generation on a given surface. The section concludes with discussion over the unification of the present methods for IGM representation under the present framework. In Sec. 4, the asymptotic formulation for analysing IGM mechanical behaviour is firstly introduced. It is followed by the corresponding topology optimisation formulation, along with the associated sensitivity analysis results. In Sec. 5, twodimensional examples will be firstly presented with several key issues of the algorithm addressed. Threedimensional cases will also be shown in brief in the end of Sec. 5. The article concludes with further discussion in Sec. 6.
2 IGM description under the original AHTO plus framework and challenging issues
2.1 IGM description under the original AHTO plus framework
A graded microstructural configuration filling a certain domain, as shown in the bottomleft panel of Fig. 1, normally involves at least two length scales:
a microscopic length scale on which the structural details are fully resolved, and a macroscopic length scale on which the trend in microstructural variation is observed. Here we characterise the “microscale” by a length parameter , and the “macroscale” by the domain size . Practically, we have
(1) 
The mathematical representation of a structural configuration is usually realised through the introduction of a topology description function, denoted by here. For a porous region with the solid parts occupying a domain of , its associated TDF can be defined by
(2) 
As for an IGM configuration, evaluating its associated TDF piece by piece on the fine scale is sometimes too timeconsuming. For improvement, a treatment using function composition has been introduced Pantz_SIAMJCO2008 ; LiuC_JAM2017 , and further developped ZhuYC_JMPS2019 . The key is to make use of the smoothlyvarying property of an IGM configuration.
As shown in Fig. 1, a macroscopically smooth mapping function is firstly introduced, so as to map an IGM configuration into a spatially periodic one in a fictitious space, whose constituting cells are uniform and of a same length scale . Then one may pick any cell (of rectangle/cuboid shape in the fictitious space) as a generating cell, and nondimensionalise it into a region of unit dimensions, say with denoting the dimensionality number (as shown in the bottomright panel of Fig. 1). Now a TDF , can be assigned to represent the structural configuration within the (nondimensional) fictitious unit cell. Finally, the TDF characterising the actual porous region should be obtained by a composition of the (macroscopically smooth) mapping function and the TDF associated with the (fictitious) porous region in :
(3) 
The TDF generated through function composition based on Eq. (3) has shown many attractive features. Firstly, by Eq. (3) is a continous function, because is smooth and is periodic over its input. Thus the resulting IGM configuration always takes a smoothlyvarying profile. Secondly, Eq. (3) has a simple form, and this facilitates its implementation in a computer aided design (CAD) platform. Thirdly, it overcomes the limitation that the actual building blocks have to be of rectangle/cuboid shape. Actually, a variety of IMG configurations can be described by using Eq. (3), as demonstrated in ZhuYC_JMPS2019 . Fourthly, Eq. (3) poses low demand on digital resources needed for memorising it. In theory, the mapping function can be defined on a coarser grid, and the TDF which is of high resolution spans only a single cell. Finally, a homogenisation formulation used for analysing the mechanical behaviour of an IGM configuration generated by Eq. (3) can be derived through mathematically rigorous asymptotic analysis ZhuYC_JMPS2019 , as tends to 0. Thus the stress analysis of an IGM configuration can still be performed at a coarsegrained level, and this renders significantly bring down the requirement on computational resources, when compared to simulations directly conducted on a fine scale.
2.2 Challenging issues
Albeit many attractive properties as mentioned above, the TDF generated by function composition through Eq. (3) still sees severe limitation. For instance, it struggles to describe an IGM configuration, where the microstructural volume fraction varies in space. This is because the mapping function simply rescales, rotates and distorts a single generating cell configuration. Thus the volume fraction stays unchanged after being transformed by . For a detailed illustration over this point, one may refer to Appendix B in ZhuYC_JMPS2019 . Thus the present study is motivated: how to expand the coverage of the AHTO plus framework for describing more general IGM configurations.
3 Representation of general IGM configurations
In this section, we will first discuss on how to extend the function composition treatment described in Eq. (3) for a much more general use. Then two of its generalised applications will be discussed in depth. We first introduce an algorithm for generating IGM configurations smoothly patching various seminal cells. Then the formulation is further extended for representing graded microstructures decorated on a given manifold. Finally, we will show that the present formulation effectively offers a theoretical framework that unifies various ways of IGM representation in literature.
3.1 General formulation
For the purpose of greatly improve the variety of describable IGM configurations using the function composition treatment, the TDF given by Eq. (3) is now modified as
(4) 
where , termed as the structural generating function in this article, is periodic over its input; is another smooth function in . A similar treatment has also been used in GuoX_CMAME2013 to account for possible boundary perturbations in topology optimization considering shape uncertainties.
The TDF defined by Eq. (4) can be rationalised in a viewpoint of levelset functions. In a traditional levelsetbased topology optimisation framework WangYu_CMAME2003 ; Allaire_JCP2004
, a levelset function is defined within the whole design domain, where only its zerolevelset contour has a geometric meaning of representing the structural boundaries. This idea is now generalised through Eq. (
4), where all levelset values are geometrically meaningful. Practically, one may introduce a function , which is periodic over the fictitious nondimensional domain of unit dimension measured by . As shown in Fig. 2, its corresponding levelset contour of height , given by , specifies a certain structural boundary. Since is continuous, the change in boundary profiles as varies should be smooth, too. Hence a continuous menu of generating cell patterns are obtained in the fictitious (nondimensional) space measured by . Then for a point in the actual space, one needs to specify and . Here we let . Thus Eq. (4) indicates that the actual microstructural configuration around is generated by manipulating the generating cell identified by the levelset contour of of height .The IGM configurations generated from Eq. (4) naturally inherit all the attractive features from the original treatment with function composition by Eq.(3). For instance, the resulting IGM configuration varies smoothly in space, because and are both continuous over their inputs and is periodic over its input.
The transition between generating cells is thus enabled, and so does the volume change in microstructures. To formulate the microstructural volume fraction, we first identify the correlation between the contour of of height and the volume fraction associated with the corresponding generating cells. Mathematically, this is given by
(5) 
where is the Heaviside function. It is noted that the correlation formulated by Eq. (5) is already set as the continuous menu of generating cells are present, which proceeds the representation of the final IGM configurations. Besides, cell manipulation, involving rescaling, rotation and isovolumetric distortion, carried out through , does not alter the corresponding local volume fraction. Thus the microstructural volume fraction around a spatial point solely depends on the value of the levelset indicator function , but is independent of . Therefore, the overall volume fraction of an IGM configuration characterised by Eq. (4) is formulated by
(6) 
where is given by Eq. (5); “” denotes the volume occupied by the corresponding domain.
In theory, the significantly expanded IGM describability by Eq. (4) is demonstrated. In the following, we will further explore its applicability in a more practically meaningful context.
3.2 IGM generation from seminal cell patterns
A number of microstructural configurations have been proposed for achieving exceptional functions (e.g. Sigmund_APL1996, ; Frenzel_Science2017, ). One extended application of Eq. (4) is to smoothly transit between cells outputting different highend performances. For instance, two seminal cell structures are present: a smiling face and a weeping face, as shown in Fig. 3.
Then the question is how to design an IGM configuration where the facial expression continuously turns from smiling to weeping. Note that this is by no means a trivial task. For instance, the two seminal structures are topologically different. Hence it is of difficulty to use a single set of parameters to consistently express both configurations. Besides, the heights of “cheeks” of the two faces are different. Using spatially discrete parameters to patch them, e.g. WangYQ_CMAME2017 ; ChengL_CMAME2019 , may result in persistent mismatch between neighbouring cells.
Suppose the two facial configurations shown in Fig. 3 are identified by the zerolevelset contours of two functions, and
, respectively. Then smooth connection between the two configurations should be obtained through interpolation between them. Hence we first construct a structural generating function by
(7) 
where is the interpolation parameter. Compared to the structural generating function constructed in the previous subsection, defined by Eq. (7) contains an extra parameter . This time, the structural boundaries of the generating cell patterns are still identified by the implicit relationship given by . It is checked that when , becomes , corresponding to the structural boundary of the seminal cell given by Fig. 3(a). So on for the case when , corresponding to the weeping facial configuration identified by .
By first letting in Eq. (7) and incorporating it into Eq. (4), we derive the overall TDF associated with an IGM configuration by
(8) 
Although coming from Eq. (4), the TDF expression defined by Eq. (8) has a clear mathematical meaning. It is generated through direct interpolation between the TDFs identifying the two seminal cell configurations, which is controlled by the levelset indicator function .
With reference to Eq. (8), Fig. 4 is generated, where the mapping function is set with , and is a linear function of , the spatial variable indicating the horizontal direction.
Note that tiny structural evolution takes place within every single cell. Thus the cells constituting the final IGM configuration are normally not exactly the same as any of the generating cells. This can be visualised in the third face to the left in Fig. 4, where tear drops only appear beneath its left eye. This again addresses the key different feature of the present method, as compared with that in WangYQ_CMAME2017 , where the constituting cells have to coincide with one of the prototype cells and tiny mismatch persists across cell boundaries.
When the effect from the mapping function is taken into account, a wall of IGM configuration with mixed expressions, for instance, can be generated, as shown in Fig. 5.
As a mapping treatment is now involved, several faces get rotated or distorted in the output IGM configuration.
The treatment can be extended to smoothly connect a number of seminal configurations. Suppose we have seminal cell configurations, represented by the zero level sets of , , , , , respectively. Now we introduce a number set , and our goal is to generate a TDF such that it becomes when the interpolation parameter equals . This can be achieved in a similar sense as Lagrangian interpolation, and a general TDF expression is thus defined by
(9) 
where the coefficients are formulated by
(10) 
It can be checked that when , the solid region is identified by , which actually corresponds to the th seminal cell configuration. Therefore, the proposed method can be used for generating IGM configurations from any given set of seminal cell patterns.
In fact, Eq. (9) may be extended to a more general form for describing a continuous menu of parameterised cell patterns, say identified by the generating function , where is a parameter. The corresponding TDF can be given by
(11) 
By this mean, the pool of generating cells can be enriched by “mixing” seminal cells at a certain ratio given by .
3.3 Generation of surfaces with smoothlyvarying microstructures
The method can be naturally generalised for deploying microstructures on a given manifold. A surface in three dimensions can be parameterised by two free variables, say, and , and it can thus be represented by
(12) 
where
are the vector containing the coordinates of the points lying on the manifold;
is a vectorvalued function from to . Here a superscript “” is attached to a variable indicating its affiliation with a manifold. Now our target is to decorate it with smoothlyvarying microstructures. For implementing the proposed method, we need to determine the TDFs of the generating cells and the mapping function .This time, the (continuous) menu of generating cells should be planar structures in twodimensional space, and suppose they are given by a periodic structural generating function , where the levelset of value identifies a (twodimensional) cell configuration. Then we consider how to construct the mapping function so as to decorate the menu of cells on the manifold given by Eq. (12). It takes the following steps.
First, since is actually a map of
, an extra degree of freedom (in conjunction with
and ), say , should be introduced. A natural way of thinking is to let it be in alignment with the surface normal , which is mathematically given by(13) 
where , and , denote the derivative of each entry of with respect to its th entry (e.g. ), and the symbol “” denotes the cross product of two vectors. Then a map between and all points around the manifold is established:
(14) 
where the subscript indicates the two inputs of the corresponding function. Then we choose , and . A mapping function from the fictitious space of to the actual space of is defined:
(15) 
and the desired mapping function appearing in Eq. (4) should be given by . Normally, an actual manifold is assigned with a thickness, say , and this can be controlled by constraining the range of by letting . Finally, one may use the treatment of function composition to define the TDF for a given manifold decorated with microstructures (generated from , ) by
(16) 
where is an arbitrary positive number; denotes the th entry of . As from Eq. (16), the TDF values for all points falling outside the thickness of the manifold is negative, and no solid materials are distributed there according to Eq. (2).
One illustrative microstructural surface is shown in Fig. 6.
The generating cell has a circular hole in its middle, which is given by with . The manifold is chosen as a parabolic surface given by . Following the aforementioned steps, the map between the fictitious and the actual spaces is found to be
(17) 
Incorporating this into Eq. (16) and letting , , , we get the microstructural surface in Fig. 6.
It is noted that for a same manifold, different parameterisations lead to disparate microstructural surfaces. For instance, both
(18a)  
and  
(18b) 
correspond to a same hyperboloid. But they give rise to different manifold decorations, as shown in Fig 7.
3.4 Comparison among existing methods for IGM representation
The present treatment based on Eq. (4) effectively offers a unified theoretical framework, which encompasses most existing approaches for IGM representation. Actually, Eq. (4) (and its extended form by Eq. (8)) enables two sets of operations: selection of generating cell configurations by the value of , and the manipulation of the selected cell (involving cell rescaling, cell rotation and cell distortion in an isovolumetric manner) controlled by the mapping function , or alternatively, by the corresponding Jacobian matrix given by
(19) 
In theory, controls the operation exerted on a single cell, and controls the transition between different cells. In this viewpoint, most existing treatments for IGM generation can be considered as being carried out through varying and , and a summary is made in Table 1.
Smooth  
Cell shape  ness  
Classical  Oriented  Piecewise  
HBMs  cuboid  No  constant  
Oriented  Slightly  Discrete  
MPCs  cuboid  not  step function  
Oriented  Slightly  Discrete  
MLSMs  cuboid  not  step function  
De  Rotational  
homogenisation  cuboid  Yes  Continuous  
Original  
AHTO plus  parallelepiped  Yes  Contant 
denotes an orthogonal matrix.
Here the reported IGM generation approaches in literature are in general classified into five groups:

classical homogenisationbased methods (HBMs) Bendsoe_CMAME1988 ; Bendsoe_StructOpt1989 ; Rodrigues_SMO2002 ; Coelho_SMO2008 ; NiuB_SMO2009 where each pixel is associated with a microstructural cell;

methods with parameterised cells (MPCs) Radman_JMaterSci2013 ; WangC_SMO2018 ; ChengL_CMAME2019 ; LiQH_CompStruct2019 where the constituting cells are formed by parameterised building blocks;

multiple levelset methods (MLSMs) WangYQ_CMAME2017 ; ZhangY_SMO2019 where the constituting cells are identified by multiple levelset contours of discrete values from a same generating function;

dehomogenisation methods Groen_IJNME2018 ; Allaire_CompMathAppl2019 ;

the original AHTO plus method ZhuYC_JMPS2019 ; XueDC_arxiv2019 .
For the first three groups of methods, the constituting cells must take a rectangle/cuboid shape. In a viewpoint of Eq. (4), each is only a function , for , , , or alternatively, the associated Jacobian matrix must be diagonal, i.e., , with a diagonal matrix. In contrast, with the projection algorithm Pantz_SIAMJCO2008 used in dehomogenisation methods, cell rotation is also permitted. In the context of Eq. (4), with an orthogonal matrix. And for the original AHTO plus method, the choice of is more arbitrary provided that . But in theory, a pool with should suffice for optimising the IGM mechanical behaviour, and this is well demonstrated by Groen et al. Groen_IJNME2018 . It should also be noted that, the edges of the constituting cells, in categories 4) and 5), are slightly curved in the actual space so as to accommodate smooth connection across cell boundaries.
The IGM smoothness is also controlled by the continuity in function appearing in Eq. (4). For 1) classical homogenisationbased approaches, each pixel is individually assigned with a microstructure. Thus is effectively piecewisely defined over individual pixels. Therefore, nonsmoothness emerges in the resulting IGM configurations. For category 2) MPCs, the microstructure is generated through gradually varying the parameter values of its building blocks, but stays the same in every single cell. Hence should be a step function, where the gaps between neighbouring steps are tiny. Consequently, the resulting IGM configurations appear graduallyvarying, but slightly nonsmooth connection still emerges across cell boundaries. Similar case arises for the category 3) MLSMs, where the microstructure is generated through gradually controlling the outputs of a single levelset function. In contrast, the dehomogenisation methods enable smooth transition across cell boundaries in the actual space, which should correspond to a continuous function . As for the original AHTO plus method, no transition between cells is permitted, giving rise to a constantvalue function of . But this limitation is removed in the present article by Eq. (4).
Therefore, most existing methods for IGM generation should find their corresponding formulation through Eq. (4) by selecting appropriate and .
4 Asymptotic formulation towards optimising IGM compliance
In this section, we will first consider how to formulate the mechanical behaviour of IGM configurations obtained based on Eq. (4). Then an optimisation problem for IGM compliance is established, with the corresponding sensitivity analysis being discussed at the end of the section.
4.1 IGM mechanical behaviour
Under the original AHTO plus framework ZhuYC_JMPS2019 , the multiscale mechanical response of a loaded IGM configuration can be predicted through asymptotic analysis. This is carried out by approximating the (multiscale) displacement field by a series expanded in terms of the small parameter defined by Eq. (1). Each term in the asymptotic expansion is a function containing two sets of spatial variables, say and . The variable captures the homogenised behaviour as if the IGM configuration is envisaged as a continuum. The (nondimensional) variable captures the large displacement gradient accommodating the local microstructures. Through asymptotic analysis, the originally multiplescale governing equations are solved on different scales individually.
In the following, the key equations constituting the original AHTO plus formulation ZhuYC_JMPS2019 are outlined. A weak form of the equilibrium equation is set up for the homogenised displacement field , given by
(20a)  
where denotes the effective elasticity tensor of the homogenised continuum; is a virtual displacement field belonging to the functional space of ; denotes the boundary on where displacement boundary condition is imposed; is the boundary on where a surface traction of density is applied. 
Meanwhile, a set of characteristic functions
, , , , , , are also needed, and they are obtained by solving the following set of PDEs:(20b) 
subjected to periodic boundary conditions over , where is the Jacobian matrix defined by Eq. (19); is a piecewisely constant function: it equals the corresponding elasticity tensor at a solid point, and vanishes elsewhere in . It is noted that the spatial variables of are , while simply serves as a parameter. But this means one has to solve Eq. (20b) individually at every single macroscopic point in the actual space.
When the solutions to Eq. (20b) become available, one may finally evaluate the effective elasticity tensor needed for solving Eq. (20a):
(20c) 
where is recalled to be the solid part in .
Such a scaleseparation treatment roots in the fact that an IGM configuration from the original AHTO plus framework becomes periodic in the fictitious space measured by , or equivalently by . For the present case, however, the generated IGM configurations are not exactly periodic in the nondimensional fictitious space measured by . Here we still adopt the same asymptotic formulae as given by Eqs. (20a)  (20c). The accuracy of the treatment here will be demonstrated through a comparison with the corresponding finescale simulation results in Sec. 5.
4.2 Topology optimisation
In this subsection, based on the results above, a mathematical formulation for optimising IGM compliance can be set up as follows
(21)  
Subject to  
where denotes the functional space of continuous functions; is a certain functional space formed by periodic functions defined over the nondimensional unit cell ; is the homogenised structural compliance; is the upper limit of the IGM solid volume. The last inequality in Eq. (21), as from Eq. (6), formulates the volume fraction constraint.
Eq. (21) establishes a scheme for minimising IGM compliance. Compared to the original AHTO plus framework, it exhibits several advantageous features. Firstly, the available IGM configurations are far more abundant. This is because one can not only conduct IGM design by rotating and distorting a generating cell through the mapping function , but also select among generating cells through the levelset indicator function . Secondly, under the original AHTO plus framework, one can not increase the volume of solid materials to strengthen the local microstructure, because the volume fraction of all constituting cells must be identical to each other. But this limitation is fully removed now. It will be shown in Sec. 5 that enabling transition among different transition cells should bring down the compliance value of an optimised structure by almost a half.
4.3 Gradientbased sensitivity analysis
For speeding up the optimisation process based on Eq. (21), its associated sensitivity analysis has been performed. Suppose that the vector quantifies all the design variables. Through using the idea of the adjoint method LiuST_SMO2002 , the derivative of the IGM compliance with respect to the th design variable can be formulated by
(22) 
Thus one simply needs to compute the derivative of the elasticity tensor with respect to all the design variables. Here the design variables are classified into two groups: the macroscopic design variables parameterising the mapping function and the levelset indicator function ; the microscopic design variables parameterising the structural generating function .
For sensitivity analysis with respect to typical macroscale variables, say, , we follow XueDC_arxiv2019 to write down
(23) 
where
(24) 
and
(25) 
respectively.
As for the sensitivity with respect to the microscale variables, say, , we also follow XueDC_arxiv2019 to write down
(26) 
where , , , , , is the Kronecker delta.
5 Numerical results and analysis
In this section, numerical results based on the proposed scheme are presented. Twodimensional examples will first be examined in depth, so as to reveal the key features of the present method. Then threedimensional numerical examples will also be briefly studied at the end of the section.
5.1 General problem setting
For all the examples examined in this section, linearly elastic material response is assumed with the isotropic parameters chosen as follows: the (nondimensional) Young’s modulus and Poisson’s ratio . For twodimensional examples, the specimen’s (nondimensional) thickness in the third dimension is fixed to be 1, and scenarios of plane stress are considered. The upper limit of IGM volume fraction is set to be 0.3 unless otherwise specified.
Here the continuous menu of generating cell configurations, once chosen, are set fixed during the optimisation processes. The reason for this action is that IGM generation based on Eq. (4) should output abundant configurations. This is in particular reasonable, when IGM configurations are produced through continuously patching seminal cell structures. Thus the evolution of the design variables parameterising the structural generating function is temporarily suspended here. As for the parameterisation of the mapping function , we follow XueDC_arxiv2019 to assume
(27) 
where symmetric conditions and are imposed.
As for the levelset indicator function , we let
(28) 
where for , , , , is required.
Upon this setup, the design variables are collected by
(29) 
Hence we have 24 design variables for twodimensional problems, and 67 design variables for threedimensional cases. Here the lengthscale characterising parameter is set to be 0.05, and the method of moving asymptotes (MMA) Svanberg_IJNME1987 are employed as the optimiser in this work.
Note that the homogenisation formulation for stress analysis involves individually solving a cell problem defined by Eq. (20b) at every single macroscopic point . This significantly brings up the computational cost. For overcoming this issue, a zoning scheme XueDC_arxiv2019 is employed, where the design domain is partitioned into several subdomains. In each subdomain, a point is selected and the effective elasticity tensor at that point obtained through Eq. (20b) is used to represent that of the whole subdomain. Here all design domains are divided into rectangles/cuboids, and the center of each rectangle/cuboild is chosen as the representative point. It has been shown XueDC_arxiv2019 that, when a subdomain contains no more than 25 finite elements (for twodimensional cases), the deviation from the corresponding finescale results is less than 0.04. We will resume using this criterion for conducting simulations for this work. As for finite element analysis, fournode quadrilateral elements are used in all twodimensional examples while eightnode hexahedral brick elements are adopted in all threedimensional examples.
5.2 Twodimensional cases
5.2.1 Laminar IGM configurations
We start with an example of a short beam fixed on its right side, as shown in Fig. 8.
A uniformly distributed compression of magnitude 0.1 is applied over its top side. The design domain is of size , and it is divided into identical squareshape subdomains. A mesh is used for solving the homogenised problem.
The continuous menu of generating cells are selected as a set of “X”shape configurations, which are controlled by the thickness of the solid region. Mathematically, all generating cells can be represented by the multiple levelset contours of a single function:
(30) 
and the structural boundaries of several generating cell configurations are shown in the left panel of Fig. 9.
Then according to Eq. (4), the TDF associated with the corresponding IGM configuration is thus given by
(31) 
Note that the range of by Eq. (30) defines the bounds for , that is, . In practice, if an output of Eq. (28) exceeds , we require it equal , and so does an output falling below . Under the present setting, the overall IGM volume fraction, according to Eq. (6), is formulated by
(32) 
The optimised result is then shown in the right panel of Fig. 9, and the IGM volume fraction is roughly 29.98%. The structural compliance based on the homogenisation formulation is computed to be 126.60, which delivers a 2.0% deviation from the corresponding finescale simulation result (124.24). Fig. 10 shows the evolution of both the compliance value and the corresponding volume fraction during optimisation. Compared to the initial structure composed of a spatially periodic microstructure, the structural compliance value drops more than 90%. It is noted that the start of evolution of the levelset indicator function is 30 steps later than that of the mapping function . This is because for minimising structural compliance, rotation and distortion of a single cell seems more effective than simply transiting among generating cells. This will be discussed with greater details later.
Compared to the spatially periodic configuration, the optimised one arranges the constituting fibres such that the external load on the top side of the design domain smoothly transited to the fixed boundary. Besides, the material distribution is denser at the topright corner, where more materials are needed for accommodating the relatively high local bending moment.
Eq. (4) effectively generates an IGM configuration combinatively by two means. One is through distorting and rotating a single cell by the mapping function , and the other is through transition among generating cells by the levelset indicator . Their individual roles in compliance optimisation are examined with the results summarised in Fig. 11 and Table 2.
Periodic  Cell manipulation  Cell transition  Both  
structure  enabled only  enabled only  enabled  
Compliance  1601.91  205.15  742.61  126.60 
Volume  
fraction (%)  30.00  30.00  29.99  29.93 
Optimised  
design  Fig. 11(a)  Fig. 11(b)  Fig. 11(c)  Fig. 11(d) 
It is shown that optimisation by simply permitting rotation and distortion of a single cell reduces the compliance value to 1/8 of that of the periodic structure shown in Fig. 11(a). This is in contrast to a drop of 1/2 if only transition between generating cells is enabled. This means arranging laminar orientations in align with the local directions of principal stress provides a more effective mean than simply redistributing materials by means of cell transition. However, compared to the optimised configuration in Fig. 11(b), further enabling transition among the generating cells results in an extra 50% drop in the optimised compliance. This demonstrates the necessities of involving transition among different cell configurations during IGM generation.
5.2.2 IGM generated through interpolating between seminal cells
In this case, the continuous menu of generating cells are obtained by interpolating between two seminal cells, which are chosen as a full solid square with a diminishing circular hole at its centre and a unit square with a hyperelliptic hole in its middle. Mathematically, a fullsolid configuration with a diminishing hole can be identified by the zero levelset contour of
(33) 
The cell with a hyperelliptic hole at its centre is indicated by the zero levelset contour of
(34) 
Thus based on Eq. (7), one introduces a structural generating function by interpolating and controlled by a parameter :
(35) 
such that the implicit relation corresponds to the structural boundary of one generating cell configuration in the fictitious space measured by . The structural boundaries of several generating cells are shown in the left panel of Fig. 12.
It can be checked that when is near 0, casts a dominant weight over . In this scenario, identifies a single point, and this corresponds to a square of full solid. As keeps growing, the hole shape looks more like a mixture of a circle and a hyperellipse. When , becomes , corresponding to a cell with a hyperelliptic hole. Note that the transition in the hole shape is carried out smoothly, since given by Eq. (35) is a continuous function (in ).
The volume fraction of the generating cells, obtained based on Eq. (35), spans from 100% () to 3.81% (). In contrast to the previous case, no explicit expressions linking a generating cell with its corresponding volume fraction are available. Nevertheless, one may precalculate a data set correlating the two quantities, and fit for a simple function approximating the volume fraction of generating cells against the parameter by
(36) 
Finally, by introducing the mapping function by , and the levelset indicator function , we express the TDF of the IGM configurations for topology optimisation:
(37)  
which is effectively Eq. (8). The volume constraint, based on (6) and (36), is thus formulated by
(38) 
The corresponding optimised IGM design is shown in the right panel of Fig. 12 with the volume fraction being 29.94%. The homogenised result for structural compliance is 124.40, which delivers a deviation of 2.1% from the finescale computation (121.73). Fig. 13 shows the evolution in the values of both the structural compliance and the corresponding volume fraction during optimisation.
Again, the individual effects of cell manipulation (through ) and cell transition (through ) are examined for this case. The results are summarised in Fig. 14 and Table 3, and a conclusion similar as that in Sec. 5.2.1 can be drawn.
Periodic  Cell manipulation  Cell transition  Both  
structure  enabled only  enabled only  enabled  
Compliance  717.21  192.85  485.87  124.40 
Volume  
fraction (%)  30.01  30.01  30.00  29.94 
Optimised  
design  Fig. 14(a)  Fig. 14(b)  Fig. 14(c)  Fig. 14(d) 
Finally, we compare in Fig. 15 the optimised designs based on the two sets of generating cells in Sec. 5.2.1 and here, so as to summarise for some features of the present approach.
Several aspects of similarities are read. Firstly, microstructural conduits transiting the load to the fixed boundaries are seen in both optimised configurations. Secondly, the trends in volume change are the same. Thirdly, both optimised configurations output similar compliance values. But the design given by Fig. 12 tends to reveal how a tradeoff is made during optimisation. It has been widely agreed that Allaire_book2002 , for compliance optimisation, laminate microstructures should deliver a nearoptimal performance, provided that the laminate orientation properly align with the desired principal stresses. Compared to circular holes, configurations with a hyperelliptic hole are morphologically closer to laminate structures. Fig. 12 indicates laminate microstructures are more desirable away from the fixed boundary and the loaded edge, while denser material distribution close to the fixed and loaded boundaries seems to be more effective for compliance optimisation.
5.3 Three dimensional examples
In this subsection, topology optimisation results of threedimensional cases are briefly presented, and further investigations are highly desired for its improvement. Here we consider a short beam of size with its right side fixed, as shown in Fig. 16.
On its upper surface, a uniformly distributed load is applied. The homogenised design domain is partitioned into identical cubic subdomains. For FE analysis, a mesh of is adopted.
The first example concerns a continuous menu of generating cells, whose representatives are shown in Fig. 17.
Every cell consists of three identical pillars, with each pillar crossing the other two in their middles. The crosssectional shape continuously changes from a unit circle to an almost diminishing hyperellipse, as the volume fraction drops. Mathematically, the solid parts of these generating cells can be described using a parameter :
(39) 
Starting from an initially periodic structure shown in Fig. 18(a), the optimised IGM configuration is given in Fig. 18(b). In a side view, the threedimensional optimised result bears a strong similarity with the corresponding twodimensional configuration shown in Fig. 12, that is, the microstructural cells arrange themselves for more effectively transiting the load to the fixed boundary.
Note that for better visualising the results, cell rotation is only permitted within the planes perpendicular to the direction of the design domain.
Another set of generating cells are also concerned with some of their representatives shown in Fig. 19.
This time, every generating cell is obtained by digging holes from a solid cube. The void region can be envisaged as the union of three cylinders, one crossing the other two. The crosssection of each cylinder smoothly transits from an almost diminishing circle to a hyperellipse of size 1, as the volume fraction drops. Mathematically, the solid regions within the generating cells can be described by
(40) 
with .
Starting from an initially periodic configuration in Fig. 20(a), the optimised configuration is shown in Fig. 20(b), where the microstructural orientation obtained is roughly the same as that in Fig. 18(b).
6 Conclusion and discussion
The present article proposes an algorithm for generating infill graded microstructural configurations from a continuous menu of simple generating cells. With the present scheme, the IGM profile is controlled through the introduction of two functions: a levelset indicator function identifying which cell from the menu is of concern around a macroscopic point ; and a mapping function indicating how the cell is rotated and distorted locally. With these two macroscopically continuous functions, IGM configurations can be generated based on TDF expressed by Eq. (4). Eq. (4) effectively offers a unified expression where most existing methods for IGM generation find their traces, with detailed information summarised in Table 1. The present scheme also sees its generalised applications in two aspects: 1) IGM generation through interpolation among a set of desired seminal cell configurations; 2) IGM configurations decorated on a given manifold. The IGM mechanical behaviour can then be predicted using asymptotic analysis, and the accuracy of the homogenised formulation is further demonstrated through comparisons with the corresponding finescale simulation results. The sensitivity analysis is also studied, and a (homogenised) topology optimisation formulation for IGM design is also proposed. Numerical examples for twodimensional cases are systematically analysed, and threedimensional results are also briefly shown.
The present approach brings about considerable expansions in IGM describability under a topology optimisation framework underpinned by asymptotic analysis, which is termed as an AHTO plus framework here. Compared to its original version ZhuYC_JMPS2019 , the IGM configurations obtained from the present generation scheme inherit all the attractive IGM features, such as smooth connection and simple form in TDF expression. Moreover, IGM configurations with other practically useful features, such as smooth change of volume fraction in space, are also enabled by the present algorithm. As discussed in the beginning of this article, the original AHTO plus method ZhuYC_JMPS2019 needs essential improvements mainly over three issues: 1) to substantially raise the computational efficiency when implementing the homogenised formulation for stress analysis; 2) to overcome the limitation where all constituting cells are restricted to be selfsimilar to each other; 3) to find more appropriate basis functions governing the IGM generation through Eq. (4). The present article aims for addressing the second issue mentioned above, while the first issue is partly resolved by XueDC_arxiv2019 and the third issue still needs further investigations. This also points to one promising direction for future studies.
Acknowledgement
The financial supports from National Key Research and Development Plan (2016YFB0201601) from the Ministry of Science and Technology of the People’s Republic of China, the National Natural Science Foundation of China (11772076, 11732004, 11821202), Program for Changjiang Scholars, Innovative Research Team in University (PCSIRT) and 111 project (B14013) are gratefully acknowledged.
References
 (1) Y. C. Zhu, S. S. Li, Z. L. Du, C. Liu, X. Guo, and W. S. Zhang. A novel asymptoticanalysisbased homogenisation approach towards fast design of infill graded microstructures. J. Mech. Phys. Solids, 124:612–633, 2019.
 (2) R. Lakes. Materials with structural hierarchy. Nature, 361(6412):511–515, 1993.
 (3) H. W. Dong, S. D. Zhao, Y. S. Wang, and C. Z. Zhang. Topology optimization of anisotropic broadband doublenegative elastic metamaterials. J. Mech. Phys. Solids, 105:54–80, 2017.
 (4) O. Sigmund and S. Torquato. Composites with extremal thermal expansion coefficients. Appl. Phys. Lett., 69(21):3203–3205, 1996.
 (5) M. S. Kushwaha, P. Halevi, L. Dobrzynski, and B. DjafariRouhani. Acoustic band structure of periodic elastic composites. Phys. Rev. Lett., 71(13):2022–2025, 1993.
 (6) C. Liu, Z. L. Du, Z. Sun, H. J. Gao, and X. Guo. Frequencypreserved acoustic diode model with high forwardpowertransmission rate. Phys. Rev. Appl., 3(6):064014, 2015.
 (7) N. Aage, E. Andreassen, B. S. Lazarov, and O. Sigmund. Gigavoxel computational morphogenesis for structural design. Nature, 550(7674):84–86, 2017.
 (8) A. Bensoussan, J. L. Lions, and G. Papanicolaou. Asymptotic analysis for periodic structures. Studies in mathematics and its applications. North Holland, Amsterdam, 1978.
 (9) M. P. Bendsoe and N. Kikuchi. Generating optimal topologies in structural design using a homogenization method. Comput. Methods Appl. Mech. Engrg., 71(2):197–224, 1988.
 (10) M. P. Bendsoe. Optimal shape design as a material distribution problem. Struct. Optim., 1(4):193–202, 1989.
 (11) O. Sigmund. Materials with prescribed constitutive parameters: An inverse homogenization problem. Int. J. Solids Struct., 31(17):2313–2329, 1994.
 (12) H. Rodrigues, J. M. Guedes, and M. P. Bendsoe. Hierarchical optimization of material and structure. Struct. Multidiscip. Optim., 24(1):1–10, 2002.
 (13) L. Liu, J. Yan, and G. D. Cheng. Optimum structure with homogeneous optimum trusslike material. Comput. Struct., 86(13):1417–1425, 2008.
 (14) B. Niu, J. Yan, and G. D. Cheng. Optimum structure with homogeneous optimum cellular material for maximum fundamental frequency. Struct. Multidiscip. Optim., 39(2):115–132, 2009.
 (15) J. D. Deng, J. Yan, and G. D. Cheng. Multiobjective concurrent topology optimization of thermoelastic structures composed of homogeneous porous material. Struct. Multidiscip. Optim., 47(4):583–597, 2013.
 (16) J. Yan, X. Guo, and G. D. Cheng. Multiscale concurrent material and structural design under mechanical and thermal loads. Comput. Mech., 57(3):437–446, 2016.
 (17) J. D. Deng and W. Chen. Concurrent topology optimization of multiscale structures with multiple porous materials under random field loading uncertainty. Struct. Multidiscip. Optim., 56(1):1–19, 2017.
 (18) S. W. Zhou and Q. Li. Design of graded twophase microstructures for tailored elasticity gradients. J. Mater. Sci., 43(15):5157, 2008.
 (19) A. Radman, X. Huang, and Y. M. Xie. Topology optimization of functionally graded cellular materials. J. Mater. Sci., 48(4):1503–1510, 2013.
 (20) A. Radman, X. Huang, and Y. M. Xie. Maximizing stiffness of functionally graded materials with prescribed variation of thermal conductivity. Comp. Mater. Sci., 82:457–463, 2014.
 (21) L. Cheng, J. X. Bai, and A. C. To. Functionally graded lattice structure topology optimization for the design of additive manufactured components with stress constraints. Comput. Methods Appl. Mech. Engrg., 344:334–359, 2019.
 (22) Q. H. Li, R. Xu, J. Liu, S. T. Liu, and S. Zhang. Topology optimization design of multiscale structures with alterable microstructural lengthwidth ratios. Comp. Struct., 230:111454, 2019.
 (23) Y. Q. Wang, F. F. Chen, and M. Y. Wang. Concurrent design with connectable graded microstructures. Comput. Methods Appl. Mech. Engrg., 317:84–101, 2017.
 (24) J. Gao, H. Li, L. Gao, and M. Xiao. Topological shape optimization of 3d microstructured materials using energybased homogenization method. Adv. Engrg. Softw., 116:89–102, 2018.
 (25) Y. Zhang, M. Xiao, H. Li, L. Gao, and S. Chu. Multiscale concurrent topology optimization for cellular structures with multiple microstructures based on ordered simp interpolation. Comput. Mater. Sci., 155(4):74–91, 2018.
 (26) Y. Zhang, H. Li, M. Xiao, L. Gao, S. Chu, and J. Zhang. Concurrent topology optimization for cellular structures with nonuniform microstructures based on the kriging metamodel. Struct. Multidiscip. Optim., 59(4):1273–1299, 2019.
 (27) J. Alexandersen and B. S. Lazarov. Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner. Comput. Methods Appl. Mech. Engrg., 290:156–182, 2015.
 (28) C. Liu, Z. L. Du, W. S. Zhang, Y. C. Zhu, and X. Guo. Additive manufacturingoriented design of graded lattice structures through explicit topology optimization. ASME J. Appl. Mech., 84(8):081008–1—081008–12, 2017.
 (29) P. Vogiatzis, M. Ma, S. K. Chen, and X. F. D. Gu. Computational design and additive manufacturing of periodic conformal metasurfaces by synthesizing topology optimization with conformal mapping. Comput. Methods Appl. Mech. Engrg., 328:477–497, 2018.
 (30) J. P. Groen and O. Sigmund. Homogenizationbased topology optimization for highresolution manufacturable microstructures. Internat. J. Numer. Methods Engrg., 113(8):1148–1163, 2018.
 (31) J. P. Groen, J. Wu, and O. Sigmund. Homogenizationbased stiffness optimization and projection of 2d coated structures with orthotropic infill. Comput. Methods Appl. Mech. Engrg., 349:722–742, 2019.
 (32) G. Allaire, P. GeoffroyDonders, and O. Pantz. Topology optimization of modulated and oriented periodic microstructures by the homogenization method. Computers and Mathematics with Applications, 78(7):2197–2229, 2019.
 (33) J. Wu, W. M. Wang, and X. F. Gao. Design and optimization of conforming lattice structures. IEEE Transactions on Visualization and Computer Graphics, pages 1–1, 2019.
 (34) P. GeoffroyDonders, G. Allaire, and O. Pantz. 3d topology optimization of modulated and oriented periodic microstructures by the homogenization method. J. Comp. Phys., 401:108994, 2020.
 (35) O. Pantz and K. Trabelsi. A posttreatment of the homogenization method for shape optimization. SIAM J. Control Optim., 47(3):1380–1398, 2008.
 (36) D. C. Xue, Y. C. Zhu, S. S. Li, C. Liu, W. S. Zhang, and X. Guo. On speeding up an asymptoticanalysisbased homogenisation scheme for designing gradient porous structured materials using a zoning strategy. arXiv:1909.11182, 2019.
 (37) Y. Xiang. Continuum approximation of the peachkoehler force on dislocations in a slip plane. J. Mech. Phys. Solids, 57:728–743, 2009.
 (38) Y. C. Zhu and Y. Xiang. A continuum model for dislocation dynamics in three dimensions using the dislocation density potential functions and its application to micropillars. J. Mech. Phys. Solids, 84:230–253, 2015.
 (39) S. J. Chapman, Y. Xiang, and Y. C. Zhu. Homogenization of a row of dislocation dipoles from discrete dislocation dynamics. SIAM J. Appl. Math., 76:750–775, 2016.
 (40) X. Guo, W. S. Zhang, and L. Zhang. Robust structural topology optimization considering boundary uncertainties. Comput. Methods Appl. Mech. Engrg., 253:356–368, 2013.
 (41) M. Y. Wang, X. M. Wang, and D. M. Guo. A level set method for structural topology optimization. Comput. Methods Appl. Mech. Engrg., 192(1):227–246, 2003.
 (42) G. Allaire, F. Jouve, and A. M. Toader. Structural optimization using sensitivity analysis and a levelset method. J. Comp. Phys., 194(1):363–393, 2004.
 (43) Tobias Frenzel, Muamer Kadic, and Martin Wegener. Threedimensional mechanical metamaterials with a twist. Science, 358(6366):1072, 2017.
 (44) P. G. Coelho, P. R. Fernandes, J. M. Guedes, and H. C. Rodrigues. A hierarchical model for concurrent material and topology optimisation of threedimensional structures. Struct. Multidiscip. Optim., 35(2):107–115, 2008.
 (45) C. Wang, J. H. Zhu, W. H. Zhang, S. Y. Li, and J. Kong. Concurrent topology optimization design of structures and nonuniform parameterized lattice microstructures. Struct Multidiscip O, 58(1):35–50, 2018.
 (46) S. T. Liu, G. D. Cheng, Y. Gu, and X. G. Zheng. Mapping method for sensitivity analysis of composite material property. Struct. Multidiscip. Optim., 24(3):212–217, 2002.
 (47) K. Svanberg. The method of moving asymptotes—a new method for structural optimization. Internat. J. Numer. Methods Engrg., 24:359–373, 1987.
 (48) G. Allaire. Shape Optimization by the Homogenization Method, volume 146. Springer Science and Business Media, 2002.
Comments
There are no comments yet.