De-homogenization of optimal multi-scale 3D topologies
Jeroen Groen, Florian Stutz, Niels Aage, J. Andreas Bærentzen, Ole Sigmund
DDe-homogenization of optimal multi-scale 3D topologies
Jeroen P. Groen a,1, ∗ , Florian C. Stutz b , Niels Aage a , Jakob A. Bærentzen b , Ole Sigmund a a Department of Mechanical Engineering, Solid Mechanics, Technical University of Denmark b Department of Applied Mathematics and Computer Science, Technical University of Denmark
Abstract
This paper presents a highly efficient method to obtain high-resolution, near-optimal 3D topologies opti-mized for minimum compliance on a standard PC. Using an implicit geometry description we derive a single-scale interpretation of optimal multi-scale designs on a very fine mesh (de-homogenization). By performinghomogenization-based topology optimization, optimal multi-scale designs are obtained on a relatively coarsemesh resulting in a low computational cost. As microstructure parameterization we use orthogonal rank-3microstructures, which are known to be optimal for a single loading case. Furthermore, a method to getexplicit control of the minimum feature size and complexity of the final shapes will be discussed. Numeri-cal examples show excellent performance of these fine-scale designs resulting in objective values similar tothe homogenization-based designs. Comparisons with well-established density-based topology optimizationmethods show a reduction in computational cost of 3 orders of magnitude, paving the way for giga-scaledesigns on a standard PC.
Keywords: optimal microstructures, giga-scale topology optimization, numerical efficiency, length-scaleenforcement, de-homogenization
1. Introduction
Topology optimization is an advanced design tool with the power to provide engineers with novel in-sights about optimal design. Nowadays, availability of high performance computing resources allows for theapplication of topology optimization on realistic design problems using different types of physics. Examplesinclude, optimization of heat sinks using natural convection (Alexandersen et al., 2016), optimizing fluid flowsystems by modeling turbulence using Reynolds-averaged Navier-Stokes (Dilgen et al., 2018), and acoustichorn optimization (Schmidt et al., 2016). However, the most studied (and arguably the most simple) typeof optimization problem is compliance minimization ( i.e. stiffness maximization) of linear elastic structuressubject to one or more load-cases and an upper bound on the material useage. Recently, Aage et al. (2017)extended the state-of-the-art of solution methods for these type of problems by optimizing an airplane wingusing more than 1 billion design variables. To do so, 8000 cores were employed on a high performance com-puter system up to 5 days. Hence, a significant reduction in computational cost is still required if large-scaletopology optimization is to be adapted interactively in engineering practice.Theoretically, it is known that the optimal shape for a compliance minimization problem contains periodicdetails on several length-scales (Kohn and Strang, 1986; Allaire and Aubry, 1999). Instead of modeling thesemicroscopic details on an extremely fine mesh, the problem can be decomposed into a multi-scale problem.The theory of homogenization can be used to calculate the effective macroscopic properties of these complexbut periodic microstructures (Bensoussan et al., 1978). A class of microstructures that is able to reach thebounds on maximum strain energy are the so-called rank-N laminates (Lurie and Cherkaev, 1984; Norris, ∗ Correspondence to: J. P. Groen, Department of Mechanical Engineering, Solid Mechanics, Technical University of Denmark,Nils Koppels All´e, Building 404, 2800 Kgs. Lyngby, Denmark E-mail: [email protected] a r X i v : . [ c s . C E ] O c t
2. Homogenization-based topology optimization
It is well-known that the optimal design for compliance minimization problems subject to a single loadingcase can be described by macroscopically varying orthogonal rank-3 laminates (Francfort and Murat, 1986;Avellaneda, 1987; Gibiansky and Cherkaev, 1987). These multi-scale microstructures consist of two differentmaterials, a stiff isotropic material (+) and a weak/compliant isotropic material (-) mimicking void, using2oung’s moduli E + and E − , respectively, and both with identical Poisson’s ratio ν . To visualize such alaminate, we consider the orthogonal rank-2 laminate shown in Figure 1, which is optimal for a planarproblem subject to a single load case. The first layering of this microstructure is a periodic rank-1 laminate,shown in Figure 1 (a). The orientation of the rank-1 laminate is described by the layer normal n and layertangent t . Furthermore, the relative layer width of the stiff material is described by parameter µ ∈ [0 , − µ ). The elasticity tensor E R of the rank-1 laminatecan be analytically calculated using the theory of homogenization, by assuming a periodic microstructureand perfect bonding between the two material phases, see e.g. (Allaire, 2002; Bendsøe and Sigmund, 2004). θ x / ε x / ε n t μ (1- μ (+) (-) (a) Rank-1 laminate in local frame. θ n μ (1- μ ) t x / ε (+) (Rank-1) x / ε (b) Rank-2 laminate in local frame. (c) Periodic rank-2 laminate.Figure 1: Visualization of how a rank-2 microstructure is constructed from a stiff isotropic material (+) and a weak/compliantisotropic material (-). Please note the different length-scales. The rank-1 microstructure is shown in relation to the global frame of reference ( x , x ); however, pleasenote that the depicted length-scale is x /ε with ε →
0. This means that the microstructure can be assumeduniform at the length-scale x /ε . The rank-2 microstructure is then constructed on this length-scale, ( i.e. x /ε ) by combining the stiff material phase and the rank-1 microstructure as is shown in Figure 1 (b). Bysetting n = t and t = n orthogonality of the layers is enforced. In a similar fashion as before theelasticity tensor E R of the rank-2 laminate can be analytically derived as a function of orientation angle θ and relative layer widths µ and µ .The procedure for constructing a rank-3 laminate in 3D is exactly the same; however, now each layer i has two tangents t i, and t i, to describe the plane spanned by the stiff material. The elasticity tensor fora rank-3 laminate E R can be described by, E R = E + − (1 − µ )(1 − µ )(1 − µ ) (cid:16) ( E + − E − ) − − (1 + ν )(1 − ν ) µ Λ ( n , t , , t , ) + (1 − µ )( µ Λ ( n , t , , t , ) + (1 − µ ) µ Λ ( n , t , , t , )) E + (cid:17) − . (1)The elasticity tensor is thus described by the layer normals n i and tangents t i, , t i, as well as the relativelayer widths µ i , for i = 1 , ,
3. The influence of the orientation of each layer is captured by the fourth-ordertensor Λ , Λ ( n i , t i, , t i, ) = 1(1 − ν ) n i ⊗ n i ⊗ n i ⊗ n i + 12(1 − ν ) (cid:16) t i, ⊗ n i ⊗ t i, ⊗ n i + n i ⊗ t i, ⊗ t i, ⊗ n i + t i, ⊗ n i ⊗ n i ⊗ t i, + n i ⊗ t i, ⊗ n i ⊗ t i, + t i, ⊗ n i ⊗ t i, ⊗ n i + n i ⊗ t i, ⊗ t i, ⊗ n i + t i, ⊗ n i ⊗ n i ⊗ t i, + n i ⊗ t i, ⊗ n i ⊗ t i, (cid:17) , (2)where ⊗ indicates the dyadic product. The normal and tangent vectors of the three layers are linked since weuse an orthogonal rank-3 microstructure. It is well-known that three Euler angles are required to represent3n orthogonal frame in 3D. Hence, we use angles θ , θ and θ to define the different layer normals andtangents, n = t , = t , = cos( θ )cos( θ )sin( θ ) + sin( θ )sin( θ )cos( θ )sin( θ )sin( θ ) − cos( θ )sin( θ )cos( θ )cos( θ ) , n = t , = t , = cos( θ )sin( θ )sin( θ ) − cos( θ )sin( θ )sin( θ )sin( θ )sin( θ ) + cos( θ )cos( θ )sin( θ )cos( θ ) , n = t , = t , = cos( θ )cos( θ )cos( θ )sin( θ ) − sin( θ ) . (3)This means that the elasticity tensor of a rank-3 laminate can be described by six parameters µ , µ , µ , θ , θ and θ . And it is possible to find the derivatives of E R w.r.t to these variables to perform gradient-basedoptimization. These expressions for the gradients are long but not necessarily difficult to derive; furthermore,the material volume fraction ρ of a rank-3 microstructure is defined as, ρ = µ + µ + µ − µ µ − µ µ − µ µ + µ µ µ . (4) The homogenization-based topology optimization problem will be solved on a finite element mesh dis-cretized using tri-linear finite elements, which each hold a uniform microstructure. As is shown by D´ıazand Sigmund (1995) a checkerboard-like pattern analyzed using linear finite elements can be stiffer than anoptimal microstructure with the same average density. To avoid these artificially stiff patterns, a classicaldensity filter is used to obtain filtered relative layer widths ˜ µ i from µ i (Bourdin, 2001; Bruns and Tortorelli,2001). As filter radius we use R = 1 . h c , with h c the length of an element.Additionally, we want to avoid very thin and very thick layers. Instead we want the layers to be eithervoid, completely solid or in the interval [ η, (1 − η )], with η = 0 .
05 used in this study. To do so, we use thesame interpolation scheme as proposed in (Groen and Sigmund, 2018) that links the filtered relative layerwidths ˜ µ i to the physical relative layer widths ¯˜ µ i ,¯˜ µ i = ˜ µ i (cid:0) − ¯ H ( β, (1 − η ) , ˜ µ i ) (cid:1) ¯ H ( β, η, ˜ µ i ) + (cid:32) β − β + ˜ µ i β (cid:33) ¯ H ( β, (1 − η ) , ˜ µ i ) . (5)Where ¯ H is the smoothed Heaviside function (Wang et al., 2010),¯ H ( β, η, ˜ µ i ) = tanh( βη ) + tanh( β (˜ µ i − η ))tanh( βη ) + tanh( β (1 − η )) , (6)with parameter β controlling the sharpness of the projection and η the threshold parameter. The interpo-lation curves for different values of β and η is shown in Figure 2. The order of lines in the legend shows thecontinuation approach that is taken, with 50 iterations per step. This means that the material interpolationscheme begins close to a linear function, gradually η is increased to enforce a length-scale on the relativelayer widths. Finally, β is increased to a high value to ensure that ¯˜ µ i is either 0, 1 or in the region [0 . , . It is known that a microstructure with an elasticity tensor having orthotropic symmetry conditions isoptimally aligned using the principal stress directions when a single load case problem is considered (Ped-ersen, 1990; Cowin, 1994; Norris, 2005). It is therefore appealing to align the microstructure normals n , n , n with the eigenvectors v I , v I , v III corresponding to the principal stresses. Unfortunately, solvingthe cubic equation for the 3D eigenvalue problem leads to an arbitrary order of the principal stresses as4 igure 2: Interpolation scheme plotted for the intervals where the behavior is non-linear, for different values of η and β , thatfollow the order of the continuation approach. σ I ≥ σ II ≥ σ III . This means that eigenvectors can swap 90 degrees when there is multiplicity in eigenvalues,in turn leading to instability in the optimization when the layer normals interchange. To circumvent theseinstabilities we update the orientation vectors based on the gradients w.r.t. θ , θ and θ .For the de-homogenization approach to work it is very important that n , n , n are smooth and con-tinuous throughout the design domain Ω. This is unfortunately neither possible using the principal stressdirections, nor using gradient-based optimization of the Euler angles. As possible solution one can regularizethe orientation field using the approach by Geoffroy-Donders et al. (2018); however, in our experience thismay lead to misaligned members w.r.t. the load path. To have a de-homogenized design that performswell, we propose a computational approach in which the smooth and continuous vector fields ˜ n , ˜ n , ˜ n arereconstructed from n , n , n . This approach, that will be described in detail in the next Section, requiresthat n , n , n are smoothly varying through Ω with a rotational symmetry of π/
2. Hence interchangingvectors are allowed; however, large changes in frame orientation should be avoided.We perform a regularization step on the discretized mesh and consider all faces n f connecting twoelements. A penalization function P fi ∈ [0 ,
1] is introduced that penalizes the difference between twovectors n i ( x f, ) and n i ( x f, ) that are connected by face f , P fi = 4 (cid:0) n i ( x f, ) · n i ( x f, ) (cid:1) − (cid:0) n i ( x f, ) · n i ( x f, ) (cid:1) . (7)This function has a minimum value if the vectors have the same orientation or an angle difference of kπ/ k , and a maximum value for a relative angle difference of π/ kπ/
2. To demonstrate this,consider face f and corresponding normal vectors n i ( x f, ) and n i ( x f, ) shown in Figure 3(a). The valueof penalization function P fi plotted against the inner product of the two corresponding vectors is shown inFigure 3(b).By looping over the three normal vector directions, and all faces n f connecting two elements, we canobtain a single regularization objective F θ , F θ = (cid:16) n f (cid:88) f =1 3 (cid:88) i =1 P qfi (cid:17) /q . (8)Numerical experiments have shown that a norm aggregation of q = 1 yields the best results. Finally, itshould be noted that regularization objective F θ can be used to augment the optimization objective or canbe imposed as a constraint. The sensitivities w.r.t. Euler angles θ , θ and θ can be derived analyticallyto allow for gradient-based optimization. The goal of the homogenization-based topology optimization problem is to minimize objective functional F , which is a combination of the compliance J ( i.e. the external work), and regularization objective F θ .5 i ( x f,1 ) n i ( x f,2 ) (a) Orientation vectors n i corresponding to face f . (b) P fi as a function of the inner product.Figure 3: Visualization of the penalization of difference in orientation vectors in two elements connected by face f . The domain is discretized using n x × n y × n z tri-linear finite elements, and hence the design variables canbe discretized in design vectors µ , µ , µ , θ , θ and θ . The optimization problem is solved in a nestedform, which means that for each design iteration we solve the state equation after which the design vectorsare updated. Hence, the discretized optimization problem can be written as,min µ , µ , µ , θ , θ , θ : F ( µ , µ , µ , θ , θ , θ , U ) = γ c J ( µ , µ , µ , θ , θ , θ , U ) + γ θ F θ ( θ , θ , θ ) , s.t. : K ( µ , µ , µ , θ , θ , θ ) U = F , : v T ρ ( µ , µ , µ ) − V maxf V ≤ , : ≤ µ , µ , µ ≤ , : -4 π ≤ θ , θ , θ ≤ π , (9)here v is the vector containing the element volumes and V maxf is the maximum allowed fraction of thematerial in Ω, with V the volume of Ω. K is the stiffness matrix and vector F describes the loads actingon the domain. We solve for the displacement vector U using a conjugate gradient method in combinationwith a geometrical multigrid pre-conditioner (Amir et al., 2013). For the design update the MATLABimplementation of the Method of Moving Asymptotes (MMA) introduced by Svanberg (1987) is used.As a starting guess for the layer widths we use µ = µ = µ , such that the volume constraint is exactlysatisfied. The starting guess for the orientation is based on a pre-analysis using isotropic microstructures,the corresponding principal stress directions are used to determine θ , θ and θ .Finally it should be mentioned that the scaling parameters γ c and γ θ have a large influence on theoptimization procedure. γ c is based on compliance of the first analysis step J (1) , while γ θ is based on theregularization objective for the starting guess F (1) θ , such that, γ c = 1 J (1) , γ θ = 12 F (1) θ . (10) In this paper we will use four different examples, comprising the Michell cantilever, Michell’s torsionsphere, an electrical mast example and the L-shaped beam all shown in Figure 4.For the Michell cantilever the load is applied in a distributed sense over a patch of L/ × L/ /L . Furthermore, this patch of elements is set to solid with a depth of L/
24 into thedomain. For the torsion sphere both the load and the Dirichlet boundary condition are applied on a squarewith dimensions L/ × L/
12. The load is applied as a line load along the boundary of this square using amagnitude of 3/L. Finally, there are solid elements at both boundary conditions. The electrical mast example6 L F L Michell cantilever.Michell's torsion sphere. Electrical mast. ? L F L-shaped beam. L F LL L F LL 3LL/4L/4L/8
Figure 4: Dimensions and boundary conditions for the four numerical examples used in this work. is inspired by Geoffroy-Donders et al. (2018). Like them, we only model a fourth of the full structure, hencethe red shaded boundaries represent symmetry conditions. The load is applied in a distributed sense over apatch of L/ × L/ /L . The L-shaped design domain, consists of passive elementsto show a torsion bending coupling. The load is applied in a distributed sense over a square patch of solidmaterials of L/ × L/ /L . Finally, it has to be mentioned that all examples areobtained using a Young’s modulus for the stiff material E + = 1, and E − = 10 − L = 1 and a maximum allowed volume fraction V maxf = 0 . J c are shown in Table 1. We demonstrate the differencebetween gradient-based alignment of the microstructure and the use of a stress-based alignment procedure.Furthermore, we show the effect of the starting guess for the microstructure orientation, which can be eitherbased on stress directions in a pre-analysis or by setting the Euler angles to zero. Finally, we show the effectof regularization on the angles and the regularization scheme that avoids thin features in the relative widthsof the layers. The method used in this work is denoted in boldface.As can be seen, the best results are obtained using gradient-based optimization of the orientation anglesin combination with the principal directions as starting guess. In general the effect of adding the projectionscheme (see Figure 2) to avoid thin relative widths does not increase the compliance, in fact it improvesthe performance. Furthermore, the regularization of the angles only reduces the performance by 1 − able 1: Compliance values J c for the electrical mast example, optimized using different discretizations, alignment methods,starting guesses and regularization schemes. The method used in this work is denoted in boldface. Mesh size alignment method Starting orientation Regularization method J c × ×
72 gradient-based principal directions Both widths and angles 98 . × ×
72 gradient-based principal directions Only on widths 97.1224 × ×
72 gradient-based principal directions Only on angles 99.1724 × ×
72 gradient-based principal directions No regularization 97.1724 × ×
72 gradient-based θ = θ = θ = Both widths and angles 100.1024 × ×
72 gradient-based θ = θ = θ = Only on widths 97.3924 × ×
72 stress-based principal directions Only on widths 97.69 × ×
144 gradient-based principal directions Both widths and angles 98 . × ×
144 gradient-based principal directions Only on widths 97.4448 × ×
144 gradient-based principal directions Only on angles 99.6348 × ×
144 gradient-based principal directions No regularization 98.4348 × ×
144 gradient-based θ = θ = θ = Both widths and angles 102.1348 × ×
144 gradient-based θ = θ = θ = Only on angles 98.4348 × ×
144 stress-based principal directions Only on widths 97.57
3. De-homogenization method
A design with rank-3 microstructures can be approximated on single scale using an approach similar tothe one presented in (Groen and Sigmund, 2018; Groen et al., 2019). However, some extra steps need to betaken, since the vector fields that describe the layer normals should be smooth and continuous. First of all,it is not directly possible to interchange the layer normals n i , n j and the physical relative layer widths ¯˜ µ i and ¯˜ µ j for i (cid:54) = j , due to the hierarchy contained in ¯˜ µ i Hence, before the three layers are sorted to be smoothand continuous throughout domain Ω it is important that this hierarchy is taken into account. Afterwards,we have three orthogonal vector fields that are aligned up to π/
2. However, we should note that the normalvectors n , n , n describing the microstructure orientation are rotationally symmetric up to π , i.e. using − n instead of n results in the same microstructures. Hence, we have a 3-dimensional 6-direction fieldfrom which 3 orthogonal, smooth and continuous 1-direction fields ˜ n , ˜ n , ˜ n have to be extracted. As discussed in Tr¨aff et al. (2019) for the 2D case, it is possible to approximate a multi-scale rank-3laminate with small loss in performance on a single scale. Here, we extend the idea to orthogonal rank-3microstructures in 3D. To do so we make use of the relative layer contributions p i . These layer contributionsare linked to the physical relative layer width ¯˜ µ i using, p = ¯˜ µ ρ ,p = (1 − ¯˜ µ )¯˜ µ ρ ,p = (1 − ¯˜ µ )(1 − ¯˜ µ )¯˜ µ ρ . (11)In a subsequent step we can obtain the single scale layer widths w i = αp i , with α > ρ = α ( p + p + p ) − α ( p p + p p + p p ) + α p p p . (12)8 .2. Obtaining smooth and continuous vector fields Consider the 6-direction vector field consisting of layer normals ± n , ± n , ± n and layer widths w , w , w .From this we want to extract three smooth and continuous 1-direction fields ˜ n , ˜ n , ˜ n with sorted widths˜ w , ˜ w , and ˜ w . Where the normal vectors ˜ n i are smooth and continuous throughout the domain.Similar to the method described by Geoffroy-Donders et al. (2018) we require that the vector fields usedfor the de-homogenization are free from singularities. The effect of singularities is that if you trace a curvearound them and follow a given vector as you move along this curve, then the vector will be rotated whenyou come full circle. A more detailed explanation about vector fields, direction fields, singularities andpossible treatments can be found in (Klberer et al., 2007; Nieser et al., 2011; Vaxman et al., 2016).Fortunately, we observe that singularities tend to occur in (nearly) void regions outside of the mechanicalstructure, and this points to an effective way of computing the 1-direction fields such that the singularities donot affect our results. We observe that in the absence of singularities, we can simply propagate a consistentchoice of vectors for our 1-direction fields from an initial element. If we propagate in such a way that wevisit elements which might have singularities ( i.e. nearly void elements) only after we have visited all of theelements pertinent to the mechanical structure, then we could stop once all significantly non-void elementshave been visited and it will not be possible to draw a loop containing a singularity in this region. Clearly,this does not protect us from singularities inside the mechanical structure (regardless of density) but theseseem to occur only for specific boundary conditions.Informed by the observations above, we have designed a front propagation approach that visits allelements in density sorted order. Initially, we take a starting element with all layers widths w i ∈ [0 . , . n = n , ˜ n = n , ˜ n = n , hence the widths follow such that ˜ w = w , ˜ w = w , and ˜ w = w .Then we add its neighbors to a priority queue, Q . The priority is given by | . − ρ | where smaller valuescorrespond to higher priority. Subsequently, when we take an element out of the queue, we fix its 1-directionfields in a way discussed below, and add its non-visited neighbors to Q . Finally, we mark this element asvisited in a vector V . This procedure leads to a traversal of all elements in order of density closest to 0 . e out of Q , we find the right handed frame ˜ F e that describes the layer widths,˜ F e = (cid:2) ˜ n ( x e ) ˜ n ( x e ) ˜ n ( x e ) (cid:3) . (13)There are j = 24 possible frame orientations F je that have to be tested to find the best ˜ F e . To do this,we identify the number of neighbor elements n n that already have been visited, and identify the possiblerotation matrix R je,i with respect to each of the already defined frame orientations of neighbor elements i ∈ n n , R je,i = ˜ F Ti F je . (14)The corresponding orientation angle ψ e,i,j that defines the frame orientation can be calculated as, | ψ e,i,j | = arccos (cid:32) trace( R je,i ) − (cid:33) , (15)Hence, ψ e,i,j = 0 would mean that the frame in e coincides with the frame in i for a given possibility j . Thebest orientation follows as, ˜ F e = F ke , for k = arg min j =1 ,..., (cid:88) i | ψ e,i,j | . (16)Once, we have the sorted vectors for element e , ˜ n , ˜ n , ˜ n that describe frame ˜ F e we store the correspondingwidths ˜ w i . Subsequently we remove element e from the queue and mark it as visited; furthermore, we addthe unvisited neighbors of element e to the queue and sort again based on density. Subsequently, we takethe element with the highest priority out of the queue and perform the same procedure again. This processis repeated until we have three smooth 1-direction fields.9 .3. De-homogenization of multi-scale designs Similar to the work presented in (Groen and Sigmund, 2018) we need to calculate a mapping field φ i foreach of the three layers. Using this mapping function we can create an implicit geometry description ˜ ρ i ( x )of the i - th layer, ˜ ρ i ( x ) = H (cid:32)(cid:0)
12 + 12
S { P i φ i ( x ) } (cid:1) − ˜ w i ( x ) (cid:33) . (17)Here H is the Heaviside function and S ∈ [ − ,
1] corresponds to the sawtooth function in MATLAB.Furthermore, P i is a periodicity scaling, which as will be discussed later, depends on mapping function φ i . The three implicit geometry functions for each layer can be combined to create an implicit geometrydescription of the de-homogenized structure,˜ ρ ( x ) = min (cid:40) , (cid:88) i =1 ˜ ρ i ( x ) (cid:41) . (18)Since small widths are avoided using the continuation scheme presented in Figure 2 we only need an accuratedescription of φ i in ˜Ω i , x ∈ ˜Ω i if ˜ w i ( x ) > .
01 and ρ ( x ) < . . (19)To solve for φ i we solve the following least-squares problem,min φ i ( x ) : I ( φ i ( x )) = 12 (cid:90) Ω α i ( x ) (cid:13)(cid:13) ∇ φ i ( x ) − ˜ n i ( x ) (cid:13)(cid:13) dΩ , s.t. : α i ( x ) ∇ φ i ( x ) · ˜ t i, ( x ) = 0 , s.t. : α i ( x ) ∇ φ i ( x ) · ˜ t i, ( x ) = 0 . (20)The domain is split into three parts, which dictate the weights on the objective α i and the weights on theconstraints α i , α i ( x ) = .
01 if ˜ w i ( x ) < . , . ρ ( x ) > . , x ∈ ˜Ω i . , α i ( x ) = w i ( x ) < . , ρ ( x ) > . , x ∈ ˜Ω i . (21)The term α i is introduced to relax the requirements for φ i in regions that are either solid or void, where thelow values still ensure some regularization to the lattice spacing. Furthermore, the term α i is used to turn offexact angular enforcement in these regions. Numerically, we solve the above mentioned problem using a finiteelement approach, where the constraints are enforced in an augmented setting using a penalty parameter γ φ . Furthermore, the complete sequence of topology optimization, creating mapping fields, and creating animplicit geometry description can be solved in a multi-scale manner. It is known that homogenization-basedtopology optimization can be performed on a relatively coarse mesh T c , while it can still contain a lot ofdetails. The implicit geometry function ˜ ρ is a continuous description of the de-homogenized shape as long asthe mapping functions φ i and widths ˜ w i have a continuous description. For practical purposes ˜ ρ is evaluatedusing a discrete number of points; however, on a fine mesh T f , such that h f ≤ h c / ε , which can be interpretedas the unit-cell size. To do so, we define the periodicity scaling parameter P i based on the average latticespacing in the domain of interest ˜Ω i , P i = 2 π(cid:15) (cid:82) ˜Ω i d ˜Ω i (cid:82) ˜Ω i ||∇ φ i ( x ) || d ˜Ω i . (22)10 .4. Verification and clean up of de-homogenized designs The de-homogenized designs can be interpolated on a very fine voxel grid, since an interpolation basiscan be made for mapping functions φ i and widths ˜ w i . Similar to the 2D case (see Groen and Sigmund(2018)) we can identify: 1) de-homogenized layers that consist of thin widths, 2) regions that do not carryload and 3) disconnected patches of material.De-homogenized layers consisting of thin widths can be identified by evaluating the local spacing λ i ,corresponding to layer i , λ i ( x ) = 2 πP i (cid:107)∇ φ i ( x ) (cid:107) . (23)This local spacing can be used to get a description of the actual feature size on the solid f i for layer i . f i ( x ) = ˜ w i ( x ) λ i ( x ) . (24)Similar to Groen and Sigmund (2018), we add a minimum feature size f min to the solid at places where thisfeature size was violated. To do so we obtain a modified local width ˜ w ∗ i in cases where the feature size onthe solid is violated, such that ˜ w ∗ i ( x ) = (cid:40) ˜ w i ( x ) , if f i ( x ) ≥ f min , f min λ i ( x ) , if f i ( x ) < f min . (25)This new width ˜ w ∗ i is then used in Equation 18. However, this means that the volume of the de-homogenizedshape is slightly increased. This can be alleviated by proportionally removing material in the rest of thedomain; however, this option is not considered in this study.An example of the Michell cantilever evaluated on a mesh of 960 × ×
480 elements is shown inFigure 5(a) where it can be seen that there are disconnected patches of material. These patches of materialcan be identified using the in-built connected component labeling algorithm in MATLAB ( i.e. bwconncomp ),which can identify the separate components of solid material. By only keeping the component with the largestnumber of solid voxels, the disconnected patches are removed as can be seen in Figure 5(b). Unfortunately,this scheme does not take care of solid regions that do not carry any load. To remove these solid elements,we make use of the simple iterative update scheme as proposed in (Groen and Sigmund, 2018). First weperform a fine-scale finite element analysis using a slightly modified version of the publicly available topologyoptimization code using the PETSc framework Aage et al. (2014). Second, the solid elements that have astrain energy density E lower than 10 − . of the mean strain energy density ¯ E are set to void. Third, to makesure that the length-scale f min is still satisfied after each iteration, an open-close filter operation (Sigmund,2007) is applied. These steps are then repeated until no changes are made. The final design for the Michellcantilever can be seen in Figures 5(c).The combination of removal of solid elements with a low strain energy density, followed by the open-closefilter operation, generally convergences within 5-10 iterations. Besides a clean and well-connected design,the performance of the design on the fine-scale J f is immediately known. Furthermore, it has to be notedthat for the fine-scale analysis void is modeled using a Young’s modulus E − = 10 − E + .
4. Numerical examples
To demonstrate that our approach requires moderate computational resources, we perform both thehomogenization-based topology optimization step and the de-homogenization step on a modern workstationPC using a single core MATLAB code. To be more specific, the PC uses Ubuntu 16.04.6 as operating system,contains an Intel Xeon Platinum 8160 processor, with 64GB RAM memory. Although the homogenization-based designs can be evaluated on an infinitely fine voxel grid, we have been restricted by memory size toexamples in the range of 200 million voxels. These large-scale designs have been verified on the DTU Sophiacluster using 100 nodes each containing 2 AMD EPYC 7351 processors and 128GB RAM memory. Hence,a total of 3200 cores was used. 11 a) With floating material. (b) De-homogenized design. (c) After post-processing.Figure 5: A Michell cantilever de-homogenized on a fine mesh of 960 × ×
480 using ε = 30 h f and no minimum length-scaleenforcement. The first step in the procedure is the homogenization-based topology optimization. All examples shownin Figure 4 are optimized on differently discretized coarse meshes T c . The corresponding compliance valueson the coarse optimization meshes J c , the volume fraction V cf , number of iterations n iter and the total timeused for the optimization T c are shown in Table 2. Table 2: Compliance J c , volume fraction V cf , the number of iterations n iter , and run time T c for different optimizationexamples optimized using homogenization-based topology optimization. Reference compliance values J ref obtained usingdensity-based topology optimization are shown as well. Example T c J c J ref V cf n iter T c [ hh : mm : ss ]Michell cantilever 48 × ×
24 227 .
89 269 .
763 0 .
100 400 01:23:51Michell cantilever 96 × ×
48 226 .
68 265 .
52 0 .
100 400 09:45:35Michell’s torsion sphere 48 × ×
48 12 .
33 12 .
51 0 .
100 400 04:28:35Michell’s torsion sphere 72 × ×
72 14 .
00 14 .
14 0 .
100 400 15:14:11Electrical mast 24 × ×
72 98 .
23 107 .
53 0 .
100 400 01:14:20Electrical mast 48 × ×
144 98 .
09 104 .
28 0 .
100 400 07:16:38L-beam 48 × ×
24 590 .
56 668 .
37 0 .
100 400 02:02:02L-beam 96 × ×
48 567 .
62 608 .
56 0 .
100 400 13:32:11It can be seen that even the example with the finest discretization (using approximately 1.2 milliondegrees of freedom) does not require more than 16 CPU hours on the workstation, making the computationalcost manageable. Furthermore, it should be mentioned that there is great potential for a speed up if the codeis run in parallel or in a lower-level programming language (e.g. C++) compared to MATLAB. Furthermore,it is interesting to note that in most parts of the design domain, only one or two layers have a finite width.This observation is in agreement with the work of Gibiansky and Cherkaev (1987), where it is shown thatnot all loading situations result in laminates of third rank. Finally, reference compliance values J ref fordesigns obtained using the well-known Solid Isotropic Microstructure with Penalty (SIMP) method obtainedon the same mesh sizes are shown in Table 2 as well. These reference values have been obtained using acontinuation scheme on the penalization factor ( p ). We start with p = 1 and slowly increase this to p = 4such that the design converges within 450 iterations. Furthermore, it should be noted that the density filteris used in these examples with a filter radius of R = 2 h c , where the filter is turned off in the final iterationsto allow for black and white designs. Furthermore, it has to be noted that the void material is modeledusing a Young’s modulus E − = 10 − E + . 12s expected the homogenization-based designs perform much better than the single-scale density-baseddesigns. The reason is obvious, for the multi-scale designs optimal rank-3 microstructures are used, whichcontain much more information on this coarse mesh than the isotropic material using the SIMP method.Only for Michell’s torsion sphere the values are close, this is expected since it is well-known that the optimalsolution is a closed sphere of solid material (Sigmund et al., 2016). As discussed before, the de-homogenization consist of two parts. The first part is the sorting of the vectorfields. The second part is calculating mapping functions φ i , applying an average unit-cell spacing ε and alength-scale f min . The de-homogenization is performed in a multi-scale fashion, where the vector fields aresorted on coarse optimization mesh T c , the mapping functions are calculated on T i , while the design isprojected on fine mesh T f . The total computational cost to de-homogenize each of the above-mentionedexamples to a fine mesh of approximately 200 million elements is shown in Table 3. For each example anaverage unit-cell spacing ε = 40 h f is chosen, with f min = 2 h f ; however, it is noted that changing thesevalues does not affect the computational cost. Table 3: Different mesh sizes T c , and T i , T f , and computational cost for the de-homogenization step T φ . Example T c T i T f T φ [ hh : mm : ss ]Michell cantilever 48 × ×
24 192 × ×
96 960 × ×
480 00:15:36Michell cantilever 96 × ×
48 192 × ×
96 960 × ×
480 00:45:39Michell’s torsion sphere 48 × ×
48 144 × ×
144 576 × ×
576 00:36:28Michell’s torsion sphere 72 × ×
72 144 × ×
144 576 × ×
576 01:19:44Electrical mast 24 × ×
72 96 × ×
288 384 × × × ×
144 96 × ×
288 384 × × × ×
24 192 ×
96 768 × ×
384 00:23:32L-beam 96 × ×
48 192 ×
96 768 × ×
384 01:08:06As can be seen all examples can be de-homogenized into very fine designs in less than one and a halfhours, rendering the total computational cost for obtaining very fine designs in the range of approximately2 −
17 hours, utilizing only a single core on the workstation PC.
The next step is to demonstrate the performance of the presented approach. First, we show the effectof changing minimum length-scale f min and average unit-cell spacing ε on the volume fraction of the de-homogenized design V φf . To do so, we use the Electrical mast example projected on a fine mesh of 384 × × T c are shown in Table 4.As can be seen, adding a minimum length-scale adds a significant amount of material, hence the volumeconstraint can be violated up to 0 . V max . Nevertheless, this effect can be minimized if reasonable valuesfor ε and f min are chosen. When f min is large compared to ε , a large amount of material will be addedto satisfy the minimum length-scale, which results in a large violation of the volume constraint. Generally,a good rule of thumb is that εf min >
10, such that the violation of the volume constraint stays below 10%of V max . Furthermore, it is obvious that ε has to be small compared to the domain size, otherwise thede-homogenized design cannot represent the homogenization-based design. Ideally, we would thus like tode-homogenize the multi-scale design to a finer mesh, see e.g. Groen and Sigmund (2018), since this wouldallow finer details and smaller or no violations of the volume constraint, while a minimum length-scale canstill be guaranteed, ideally related to η used in the homogenization-based topology optimization procedure.Furthermore, all the de-homogenized designs perform close to the homogenization-based solutions. Thiscan be verified by analyzing the compliance of the de-homogenized designs on the fine mesh J φ , whichis also shown in Table 4. As expected the designs which have a higher volume have a lower compliance.Nevertheless, it can be seen that all values are close to the homogenization-based compliance J c . It can be13 able 4: The volume fraction V φf , compliance J φ , stiffness per weight measure S φ and total computation time T tot of thede-homogenized electrical mast on a fine mesh of 384 × × V postf andcompliance J post and stiffness per weight S post after the post-processing step are shown. Results are shown for different T c , f min and ε . T c ε f min V φf J φ S φ T tot [ hh : mm : ss ] V postf J post S post × ×
72 24 h f h f × ×
72 24 h f h f × ×
72 24 h f h f × ×
72 24 h f h f × ×
72 32 h f h f × ×
72 32 h f h f × ×
72 32 h f h f × ×
72 32 h f h f × ×
72 40 h f h f × ×
72 40 h f h f × ×
72 40 h f h f × ×
72 40 h f h f × ×
144 24 h f h f × ×
144 24 h f h f × ×
144 24 h f h f × ×
144 24 h f h f × ×
144 32 h f h f × ×
144 32 h f h f × ×
144 32 h f h f × ×
144 32 h f h f × ×
144 40 h f h f × ×
144 40 h f h f × ×
144 40 h f h f × ×
144 40 h f h f ε results in a better performance, which is expected, since the de-homogenized designsbetter resemble the multi-scale designs. Hence, this again shows the need for an even finer resolution ofthe de-homogenized designs. Furthermore, we introduce an additional measure S φ , which is the compliancemultiplied by the volume and can be seen as a measure of stiffness per volume. The lower this value, thebetter the material usage and hence performance. The use of S φ allows us to quantitatively compare thedifferent de-homogenized designs with each other to show that an almost constant performance is achievedfor the different values of ε and f min . However, with a slight advantage for the designs without length-scaleenforcement. The values for the stiffness per volume of the homogenization-based structures are S c = 9 . S c = 9 .
809 for the coarse and slightly finer mesh respectively.It is also interesting to see that the mesh on which the homogenization-based designs are achieved canbe very coarse. This demonstrates how much information is captured by using a rank-3 parameterization.Furthermore, this shows that the total time T tot to obtain the de-homogenized designs without evaluatingthe fine-scale compliance can be as low as one and a half hour on a single PC. The times are only shown onceper optimization mesh, since the choice of ε and f min does not affect the calculation of mapping functions φ i .Section views of the de-homogenized designs for different values of ε and f min are shown in Figures 6(a),(b)and (c).From these Figures it can also be seen that the average unit-cell spacing ε has to be small enough toallow for the de-homogenized design to represent the homogenization-based design. Furthermore, it canbe seen that the de-homogenized design obtained using T c = 48 × ×
144 contains more microstructuraldetails. Finally, the post-processing method described in the previous section can be used to get rid of some14 a) ε = 24 h f , f min = 0 h f . (b) ε = 40 h f , f min = 2 h f . (c) ε = 24 h f , f min = 3 h f . (d) ε = 24 h f , f min = 3 h f ,after post-processing.Figure 6: Cut-outs of the de-homogenized designs of the electrical mast example de-homogenized on a fine mesh of 1152 × ×
384 voxels for different values of ε and f min . Example (a) is optimized on T c = 24 × ×
72, while the other examplesare obtained using T c = 48 × × × ×
384 voxels for ε = 24 h f and f min = 3 h f , after post-processing. Furthermore, a reference design using density-based topology optimization is shown. V postf , compliance J post and measure for stiffness per weight S post can be seen in Table 4 as well. Here itshould be noted that in general the post-processing step not only removes material but also improves thestiffness per weight measure. Unfortunately however, this requires several fine-scale analyses that come ata large computational cost, i.e. using 3200 cores on the DTU Sophia cluster.The performance of the de-homogenized designs of the Michell cantilever, Michell’s torsion sphere andthe L-shaped beam are shown in Tables 5-7. The values for the stiffness per volume of the homogenization-based structures are S c = 22 .
79 and S c = 22 .
67 for the coarse and slightly finer mesh respectively for theMichell cantilever, S c = 1 .
233 and S c = 1 .
400 for Michell’s torsion sphere and S c = 59 .
06 and S c = 56 . Table 5: The volume fraction V φf , compliance J φ , stiffness per weight measure S φ and total computation time T tot of thede-homogenized Michell cantilever on a fine mesh of 960 × ×
480 elements. Furthermore, the volume fraction V postf andcompliance J post and stiffness per weight S post after the post-processing step are shown. Results are shown for different T c , f min and ε . T c ε f min V φf J φ S φ T tot [ hh : mm : ss ] V postf J post S post × ×
24 24 h f h f × ×
24 24 h f h f × ×
24 24 h f h f × ×
24 24 h f h f × ×
24 32 h f h f × ×
24 32 h f h f × ×
24 32 h f h f × ×
24 32 h f h f × ×
24 40 h f h f × ×
24 40 h f h f × ×
24 40 h f h f × ×
24 40 h f h f × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 40 h f h f × ×
48 40 h f h f × ×
48 40 h f h f × ×
48 40 h f h f .
1. This means that for asmall unit-cell spacing, the minimum width can be violated and a lot of material is added. Hence, for thisexample a slightly larger value of ε is actually better. Furthermore, it can be seen that there is no noticeabledifference between the two mesh sizes on which the homogenization-based designs are obtained. Finally, itcan be seen that for one case the post-processing scheme resulted in a significantly worse compliance, whichis caused by the open-close filter step in which also load carrying material has been removed.Moving on to Michell’s torsion sphere example, we observe a large difference in compliance between16 a) De-homogenized design. (b) After post-processing. (c) After post-processing.Figure 8: A Michell cantilever optimized on T c = 96 × ×
48 elements and de-homogenized on a fine mesh of 960 × × ε = 40 h f and f min = 2 h f .Table 6: The volume fraction V φf , compliance J φ , stiffness per weight measure S φ and total computation time T tot ofMichell’s torsion sphere de-homogenized on a fine mesh of 576 × ×
576 elements. Furthermore, the volume fraction V postf and compliance J post and stiffness per weight S post after the post-processing step are shown. Results are shown for different T c , f min and ε . T c ε f min V φf J φ S φ T tot [ hh : mm : ss ] V postf J post S post × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 40 h f h f × ×
48 40 h f h f × ×
48 40 h f h f × ×
48 40 h f h f × ×
72 24 h f h f × ×
72 24 h f h f × ×
72 24 h f h f × ×
72 24 h f h f × ×
72 32 h f h f × ×
72 32 h f h f × ×
72 32 h f h f × ×
72 32 h f h f × ×
72 40 h f h f × ×
72 40 h f h f × ×
72 40 h f h f × ×
72 40 h f h f a) Full view. (b) Section view.Figure 9: Michell’s torsion sphere optimized on T c = 48 × ×
48 and de-homogenized on a fine mesh of 576 × ×
576 using ε = 32 h f and f min = 2 h f .(a) Side view. (b) Section view from the top.Figure 10: The L-shaped beam example optimized on T c = 96 × ×
48 and de-homogenized on a fine mesh of 768 × × ε = 24 h f and f min = 2 h f . the mesh size. Nevertheless, the goal of this example was to show that the optimized solution visuallyrepresents the analytical solution of Michell’s torsion sphere which is a close-walled sphere (Sigmund et al.,2016). Interestingly, this is the case; however, the de-homogenized solution contains several closed spheresof different radii as can be seen in Figure 9. The reason for the multiple spheres is the relatively high volumefraction, the rank-3 material model and the fact that a relatively coarse mesh is used for the homogenization-based topology optimization. Finally, it should be mentioned that using the post-processing procedure alarge number of passive solid elements at the boundary condition have been removed, making the comparisonwith the homogenization-based design unfair.The orientation field for the L-shaped beam actually contains a singularity at the corner with the passivevoid domain. However, since this singularity is at the boundary of the domain, the sorting algorithm managedto return smooth and continuous vector fields, although we have to note that this cannot be guaranteed forother cases including singularities. In this example the microstructures and layer orientations are rapidly18 able 7: The volume fraction V φf , compliance J φ , stiffness per weight measure S φ and total computation time T tot de-homogenized L-shaped beam on a fine mesh of 768 × ×
384 elements. Furthermore, the volume fraction V postf andcompliance J post and stiffness per weight S post after the post-processing step are shown. Results are shown for different T c , f min and ε . T c ε f min V φf J φ S φ T tot [ hh : mm : ss ] V postf J post S post × ×
24 24 h f h f × ×
24 24 h f h f × ×
24 24 h f h f × ×
24 24 h f h f × ×
24 32 h f h f × ×
24 32 h f h f × ×
24 32 h f h f × ×
24 32 h f h f × ×
24 40 h f h f × ×
24 40 h f h f × ×
24 40 h f h f × ×
24 40 h f h f × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 24 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 32 h f h f × ×
48 40 h f h f × ×
48 40 h f h f × ×
48 40 h f h f × ×
48 40 h f h f T c = 96 × ×
48. As can be seen in Figure 10(a) and(b), the optimized solution is a combination of a cantilever optimized for bending stiffness and a hollow boxoptimized for torsion.
Besides verifying that the de-homogenized designs perform similar to the homogenization-based designs,it is interesting to compare the results with a well-established density-based topology optimization method.To do so, we make use of the publicly available topology optimization code using the PETSc framework (Aageet al., 2014). The code has been slightly modified to allow for different boundary conditions, passive elementsand a continuation scheme on the penalization factor ( p ). We start with p = 1 and slowly increase this to p = 4, and stop the optimization procedure after 450 iterations when the change in design is negligible.Furthermore, it should be noted that the density filter is used in these examples with a filter radius of R = 3 h f , where the filter is turned off in the final iterations to allow for black and white designs.All four examples shown in Figure 4 are optimized on the same resolution as the de-homogenized designsdiscussed above. The corresponding fine-scale compliance J f , measure of stiffness per weight S f , volumefraction V ff , the number of iterations n iter , and run time T f are shown in Table 8. Furthermore, the designfor the cantilever, and the electrical mast example are shown in Figures 7(c) and (d) respectively. Finallyit is mentioned that the same high-performance computing cluster is used as for the fine scale validations.Hence, for each optimization example 100 nodes each with 32 cores are used.19 able 8: Compliance J f , volume fraction V ff , the number of iterations n iter , and run time T f for different optimizationexamples optimized using a density-based optimization procedure. Example Mesh size J f S f V ff n iter T f Michell cantilever 960 × ×
480 236 .
262 23 .
626 0 .
100 450 08:13:28Michell’s torsion sphere 576 × ×
576 20 .
916 2 . .
100 450 01:03:06Electrical mast 384 × × .
961 9 . .
100 450 02:16:28L-beam 768 × ×
384 585 .
114 58 .
511 0 .
100 450 04:38:50It can be seen that the compliance values of the density-based designs are close in general even 5 −
5. Concluding remarks
We have presented a highly efficient approach to obtain manufacturable and ultra-high resolution designsfrom homogenization-based topology optimization. By doing homogenization-based topology optimizationusing optimal microstructures, an optimal design can be represented on a relatively coarse mesh. A good ruleof thumb would be that at least 48 elements in each direction of the mesh are required. The subsequentlyde-homogenized designs perform within 5 −
10% compared to designs obtained using well-established density-based topology optimization. However, instead of using high performance computing facilities with morethan 3000 cores the presented designs have been obtained using a single core MATLAB process on a modernworkstation PC. Hence, the presented procedure is a first step on the way to achieve giga-scale interactivetopology optimization.Besides obtaining near-optimal designs, we have presented a method to control the shape and minimumfeature size of these de-homogenized designs ensuring manufacturability. From the results it can also be seenthat even better performing designs can be obtained when de-homogenizing on a finer mesh. This paves theway for future studies into different methods to represent the de-homogenized designs, other than the voxelgrid that has been used now. Note however that the presented procedure only works for examples wherethe homogenization-based designs are free of singularities. Hence, a natural extension would be to betterunderstand the occurence of singularities and extend the de-homogenization procedure such that it takesthese singularities into account. We are confident that this can be done to make topology optimization anintegrated part of the design process for large-scale problems.
Acknowledgments
The authors acknowledge the financial support from the Villum Foundation (InnoTop VILLUM investi-gator project). Furthermore, the authors would like to express their gratitude to the members of the TopOptgroup at DTU for valuable discussions during the preparation of this work. Finally, the authors wish tothank Krister Svanberg for the Matlab MMA code. 20 eferences
Aage, N., Andreassen, E., Lazarov, B.S., 2014. Topology optimization using PETSc: An easy-to-use, fully parallel, opensource topology optimization framework. Structural and Multidisciplinary Optimization 51, 565–572. doi: .Aage, N., Andreassen, E., Lazarov, B.S., Sigmund, O., 2017. Giga-voxel computational morphogenesis for structural design.Nature 550, 84–86. doi: .Alexandersen, J., Sigmund, O., Aage, N., 2016. Large scale three-dimensional topology optimisation of heat sinks cooled bynatural convection. International Journal of Heat and Mass Transfer 100, 876 – 891. URL: , doi: .Allaire, G., 2002. Shape Optimization by the Homogenization Method. Springer-Verlag, New York.Allaire, G., Aubry, S., 1999. On optimal microstructures for a plane shape optimization problem. Structural optimization 17,86–94. doi: .Allaire, G., Geoffroy-Donders, P., Pantz, O., 2019. Topology optimization of modulated and oriented periodic microstruc-tures by the homogenization method. Computers & Mathematics with Applications 78, 2197 – 2229. URL: , doi: https://doi.org/10.1016/j.camwa.2018.08.007 . sim-ulation for Additive Manufacturing.Amir, O., Aage, N., Lazarov, B.S., 2013. On multigrid-cg for efficient topology optimization. Structural and MultidisciplinaryOptimization 49, 815–829. doi: .Avellaneda, M., 1987. Optimal bounds and microgeometries for elastic two-phase composites. SIAM Journal on AppliedMathematics 47, 1216–1228. doi: .Bendsøe, M., Kikuchi, N., 1988. Generating optimal topologies in structural design using a homogenization method. ComputerMethods in Applied Mechanics and Engineering 71, 197–224. doi: .Bendsøe, M., Sigmund, O., 2004. Topology optimization-theory, methods, and applications. Springer, Berlin.Bendsøe, M.P., 1989. Optimal shape design as a material distribution problem. Structural Optimization 1, 193–202. doi: .Bensoussan, A., Lions, J.L., Papanicolaou, G., 1978. Asymptotic Analysis for Periodic Structures (Studies in mathematics andits applications). Elsevier Science Ltd.Bourdin, B., 2001. Filters in topology optimization. International Journal for Numerical Methods in Engineering 50, 2143–2158.doi: .Bruns, T., Tortorelli, D., 2001. Topology optimization of non-linear elastic structures and compliant mechanisms. ComputerMethods in Applied Mechanics and Engineering 190, 3443–3459. doi: .Coelho, P., Cardoso, J., Fernandes, P., Rodrigues, H., 2011. Parallel computing techniques applied to the simultaneous designof structure and material. Advances in Engineering Software 42, 219–227. doi: .Cowin, S.C., 1994. Optimization of the strain energy density in linear anisotropic elasticity. Journal of Elasticity 34, 45–68.doi: .D´ıaz, A., Sigmund, O., 1995. Checkerboard patterns in layout optimization. Structural Optimization 10, 40–45. doi: .Dilgen, C.B., Dilgen, S.B., Fuhrman, D.R., Sigmund, O., Lazarov, B.S., 2018. Topology optimization of turbulent flows.Computer Methods in Applied Mechanics and Engineering 331, 363 – 393. URL: , doi: .Francfort, G.A., Murat, F., 1986. Homogenization and optimal bounds in linear elasticity. Archive for Rational Mechanics andAnalysis 94, 307–334. doi: .Geoffroy-Donders, P., Allaire, G., Pantz, O., 2018. 3-d topology optimization of modulated and oriented periodic microstructuresby the homogenization method. Working paper or preprint.Gibiansky, L.V., Cherkaev, A.V., 1987. Microstructures of composites of extremal rigidity and exact estimates of the associatedenergy density. Report 1115. Ioffe Physico-Technical Institute, Acad of Sc, USSR. Leningrad, USSR. English translation in:Cherkaev, A.V., Kohn, R.V. (Eds), 1997. Topics in the mathematical modelling of composite materials. Birkh¨auser Boston,MA.Groen, J.P., Sigmund, O., 2018. Homogenization-based topology optimization for high-resolution manufacturable microstruc-tures. International Journal for Numerical Methods in Engineering 113, 1148–1163. doi: .Groen, J.P., Wu, J., Sigmund, O., 2019. Homogenization-based stiffness optimization and projection of 2d coated structureswith orthotropic infill. Computer Methods in Applied Mechanics and Engineering 349, 722 – 742. URL: , doi: .Kohn, R.V., Strang, G., 1986. Optimal design and relaxation of variational problems, i. Communications on Pure and AppliedMathematics 39, 113–137. doi: .Klberer, F., Nieser, M., Polthier, K., 2007. Quadcover - surface parameterization using branched coverings. Computer GraphicsForum 26, 375–384. doi: .Liu, L., Yan, J., Cheng, G., 2008. Optimum structure with homogeneous optimum truss-like material. Computers & Struc-tures 86, 1417 – 1425. URL: , doi: . structural Optimization.Lurie, K.A., Cherkaev, A.V., 1984. G-closure of some particular sets of admissible material characteristics for the problem ofbending of thin elastic plates. Journal of Optimization Theory and Applications 42, 305–316. doi: .Milton, G.W., 1986. Modelling the properties of composites by laminates, in: Ericksen, J.L., Kinderlehrer, D., Kohn, R., Lions,J.L. (Eds.), Homogenization and Effective Moduli of Materials and Media, Springer New York, New York, NY. pp. 150–174. ieser, M., Reitebuch, U., Polthier, K., 2011. CubeCover- parameterization of 3d volumes. Computer Graphics Forum 30,1397–1406. doi: .Norris, A., 1985. A differential scheme for the effective moduli of composites. Mechanics of Materials 4, 1–16. doi: .Norris, A.N., 2005. Optimal orientation of anisotropic solids. The Quarterly Journal of Mechanics and Applied Mathematics59, 29–53. doi: .Pantz, O., Trabelsi, K., 2008. A post-treatment of the homogenization method for shape optimization. SIAM Journal onControl and Optimization 47, 1380–1398. doi: .Pedersen, P., 1990. Bounds on elastic energy in solids of orthotropic materials. Structural optimization 2, 55–63. doi: .Rodrigues, H., Guedes, J., Bendsøe, M., 2002. Hierarchical optimization of material and structure. Structural and Multidisci-plinary Optimization 24, 1–10. doi: .Schmidt, S., Wadbro, E., Berggren, M., 2016. Large-scale three-dimensional acoustic horn optimization. SIAM Jour-nal on Scientific Computing 38, B917–B940. URL: https://doi.org/10.1137/15M1021131 , doi: , arXiv:https://doi.org/10.1137/15M1021131 .Sigmund, O., 1994. Materials with prescribed constitutive parameters: An inverse homogenization problem. InternationalJournal of Solids and Structures 31, 2313–2329. doi: .Sigmund, O., 2007. Morphology-based black and white filters for topology optimization. Structural and MultidisciplinaryOptimization 33, 401–424. doi: .Sigmund, O., Aage, N., Andreassen, E., 2016. On the (non-)optimality of michell structures. Structural and MultidisciplinaryOptimization 54, 361–373. doi: .Sivapuram, R., Dunning, P.D., Kim, H.A., 2016. Simultaneous material and structural optimization by multiscale topologyoptimization. Structural and Multidisciplinary Optimization 54, 1267–1281. doi: .Svanberg, K., 1987. The method of moving asymptotes–a new method for structural optimization. International Journal forNumerical Methods in Engineering 24, 359–373. doi: .Tr¨aff, E., Sigmund, O., Groen, J.P., 2019. Simple single-scale microstructures based on optimal rank-3 laminates. Struc-tural and Multidisciplinary Optimization 59, 1021–1031. URL: https://doi.org/10.1007/s00158-018-2180-3 , doi: .Vaxman, A., Campen, M., Diamanti, O., Panozzo, D., Bommes, D., Hildebrandt, K., Ben-Chen, M., 2016. Directional fieldsynthesis, design, and processing, in: Computer Graphics Forum, Wiley Online Library. pp. 545–572. doi: .Wang, F., Lazarov, B.S., Sigmund, O., 2010. On projection methods, convergence and robust formulations in topologyoptimization. Structural and Multidisciplinary Optimization 43, 767–784. doi: .Wu, J., Wang, W., Gao, X., 2019. Design and optimization of conforming lattice structures. IEEE Transactions on Visualizationand Computer Graphics , 1–1doi: ..