{'fraction_non_alphanumeric': 0.03596761318280305, 'fraction_numerical': 0.015326719124187479, 'mean_word_length': 4.916340574821212, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal and continuous multilattice embedding

Voxel-based additive manufacturing enables multilattice topology optimized structures with continuously embedded microstructures.


INTRODUCTION
Many natural structures and materials such as wood, bone, and shells synthesize spatially varying mechanical properties and structural hierarchy to achieve compelling functionalities (1). For example, teeth have not only a hard and brittle outer layer (enamel) that is prone to cracking but also a tough inner layer (dentin) that mitigates crack propagation into the tooth's interior (2); the mantis shrimp's dactyl club exhibits functionally graded elasticity and hardness that leads to high impact resistance needed for striking its prey (3); bamboo's fiber reinforcement is functionally graded radially along the culm (stalk) cross section for high bending rigidity per unit mass (4); and the internal bone structure of the cuttlefish contains a porous, layered architecture that simultaneously resists high pressures experienced in the deep sea while remaining lightweight and enabling the cuttlefish to control its buoyancy (5,6).
Many of these natural features have been borrowed to enhance the performance of engineered structures and materials. For example, reinforced concrete is a multimaterial system with properties exceeding that of the constituent materials; sandwich panels use structural hierarchy to enhance their strength-to-weight ratio; and engineered foams provide energy absorption, thermal insulation, and buoyancy in various engineering applications (7). Topology optimization provides a rational way to further elicit functionalities from such systems [e.g., artificial materials with negative thermal expansion (8,9) and structures exhibiting prescribed deformations (10)(11)(12)].
Fundamental to what we know today as topology optimization is the homogenization-based method proposed by Bendsøe and Kikuchi in 1988 (13) in which an optimized material distribution is determined using homogenized microstructural-material properties, interpolated as a function of the microstructure's porosity and orientation.
Although additive manufacturing is pushing architected materials to new limits [e.g., unit cells with 300-nm minimum feature size by two-photon lithography (26,27) and hierarchical structures spanning length scales from tens of nanometers to tens of centimeters by projection microstereolithography (28,29)], such complex microstructures have yet to be cohesively coupled with macrostructures of equally complex geometry as obtained from topology optimization. For example, most manufactured topology-optimized parts composed of architected materials have box-like macrostructures with heterogeneity attributed only to unit cell density (30)(31)(32)(33)(34). On the other hand, free-form, topology-optimized macrostructures have been embedded with a homogeneous microstructural material (35). The realization of geometrically heterogeneous architected materials embedded and manufactured within free-form, three-dimensional (3D), topology-optimized macrostructures, remains a challenge. A limiting factor has been a means to effectively represent the 3D geometry of complex, spatially varying, multiscale parts in a way that is meaningful (to a 3D printer). Multimaterial and multiscale additive manufacturing is shifting from impractical surface representations to procedural (36) or voxel representations (37)(38)(39) that are less memory intensive. Here, we integrate this current trend in multiscale additive manufacturing with multiscale topology optimization to achieve structures with unprecedented complexity at both the macro-and microscales.
In an effort to achieve free-form, multiscale structures with spatially varying mechanical properties and microstructural geometries, we propose a multimicrostructure, density-based, topology optimization formulation that simultaneously determines an optimized macrostructural geometry and distribution of a set of porous microstructural materials with distinct, possibly anisotropic mechanical behavior. We focus on volume-constrained compliance minimization problems, but the ideas readily extend to other problems [e.g., stress-constrained topology optimization (40)]. We also overcome a number of challenges in manufacturing spatially varying, multiscale, topology-optimized parts with functionally graded microstructural-material interfaces via a multimaterial slicing and continuous multimicrostructure-embedding scheme at the level of the 3D printer slices. We demonstrate the capabilities of our proposed approach using a commercially available masked-stereolithography (m-SLA) 3D printer from Prusa Research.

Topology-optimized parts, manufactured with continuously embedded microstructures
We begin by using the proposed multimicrostructure topology optimization formulation (see Materials and Methods) and continuous multimicrostructure-embedding scheme to design and manufacture a cantilever beam, with minimized compliance, f, according to the domain and boundary conditions defined in Fig. 1A.
To demonstrate that the proposed multimicrostructure-embedding scheme is capable of generating smooth, well-connected transitions between microstructures composed of unit cells with similar or dissimilar geometries that may or may not have clear connectivity, we design and manufacture the cantilever beam considering three different combinations of two microstructural materials. In Fig. 1A, we consider materials composed of octahedron unit cells with different bar diameters. In Fig. 1B, we consider a simple cubic and a truncated octahedron unit cell. In Fig. 1C, we consider a face-x and a center-x unit cell. The unit cells and associated normalized, directional tensile and shear moduli, E 11 ′ / (E ˆ  i ) and G 12 ′ / (G ˆ  i ) , obtained using computational homogenization, are provided at the left of Fig. 1 for each of the three cases. In addition to normalizing by the tensile and shear moduli of the bulk material, E and G, we also normalize by the unit cell volume fraction, ˆ  i , to capture the conflicting volume and stiffness requirements in the volume-constrained compliance minimization problem. The directional shear moduli plots represent an envelope of shear stiffness for critical orientations of shear (see section S1). In all cases, the bulk materials associated with microstructural materials 1 and 2 are limited to a domain volume fraction of ¯  1 = 0.07 ˆ  1 and ¯  2 = 0.03 ˆ  2 , respectively. In each case, a smooth and continuous transition is demonstrated in the manufactured part. Note that the topology optimization formulation allows spatial freedom within the domain, leading to organic material interfaces. As a result, the microstructural materials may require connectivity at any arbitrary cross section (not necessarily at the unit cell boundaries), which poses a challenge in obtaining well-connected structures after manufacturing. A 1D illustration of this challenge is illustrated for the three different representative combinations of microstructures in Fig. 2. In Fig. 2A, a well-connected interface with abrupt transition is shown between the two octahedron unit cells that have the same geometry but different bar diameters. In Fig. 2B, a discontinuous interface is shown between the simple cubic and truncated octahedron unit cell for which connectivity at any arbitrary interface is, in general, not guaranteed. In Fig. 2C, a disconnected interface is shown between the face-x and center-x unit cell for which connectivity is only guaranteed at the unit cell boundaries. To achieve smooth and continuous transitions between the unit cells, we define a set of transitional unit cells that are obtained by interpolating the bar diameter in Fig. 2A, interpolating the unit cell geometry in Fig. 2B, and composing the two unit cells into a series of hybrid unit cells in Fig. 2C. The transitions in Fig. 2 are shown in 1D, but in the manufactured parts, they can occur in any arbitrary direction in 3D space.
We choose the transitional unit cells to achieve smooth and continuous geometry transitions and do not directly enforce any requirements on how the microstructural-material properties vary over the transition region (see plots of unit cell geometry versus normalized modulus in Fig. 2). The transition regions demonstrate monotonically decreasing stiffness over the transition in Fig. 2A, monotonically increasing stiffness followed by some oscillations in Fig. 2B, and a stiffening effect, much like standard connections used in engineering, in Fig. 2C. The microstructural-material property transitions in Fig. 2A can be approximated by an exponential function of the form E 11 (x) = E 11 (0)e −x , where the length scale of inhomogeneity in the tensile modulus (similar for the shear modulus) is characterized by 1/ = w/ ln (E 11 (0)/E 11 (w)), with w as a given transition region length (e.g., if w = 1 in Fig. 2A, then 1/ = 0.55). Other transitioning techniques can be used to achieve desired homogenized material property transitions [see, e.g., (41)]. The geometry and homogenized material property transitions are shown in more detail in movie S3.

Multimaterial slicing and continuous multimicrostructure embedding
The manufactured parts in Fig. 1 are a result of the multimaterial slicing and continuous multimicrostructure-embedding scheme, illustrated in Fig. 3 and a flowchart in fig. S3, which overcomes several existing challenges in manufacturing multimicrostructure topology-optimized parts. Homogenized material properties are used in the multimicrostructure topology optimization step, but we refer to the output simply as multimaterial density data since the microstructural-material geometries have yet to be embedded and we process the data as if it were multiple, solid, isotropic materials until the multimicrostructure-embedding step. Figure 3 (A to D) shows various ways of representing multimaterial topology optimization density data to highlight challenges in communicating the data to a 3D printer. Oftentimes, the topology optimization density field(s) coincide spatially with a mesh used for the finite element analysis [here, we use centroids of a hexahedral (hex) mesh]. By discarding densities less than a cutoff value (typically 0.5), the density data represented directly on the underlying hex mesh lead to poor resolution (Fig. 3A). Instead, it is typical to generate a smooth isosurface of the data; however, for multimaterial, care is needed to avoid disjointed material interfaces caused by low densities resulting from the density filter (see Eq. 2 in Materials and Methods) used to regularize the topology optimization problem (Fig. 3B). To promote well-connected material interfaces and capture material variations through the volume of the part, we project the density fields onto a tetrahedral (tet) mesh generated within an isosurface of the composite density field (Fig. 3C). Furthermore, we functionally grade the interfaces by applying a convolution filter (see Eq. 2 in Materials and Methods) to the multimaterial tet data (Fig. 3D), leading to controlled mixing at the interfaces (influence of the filter radius is investigated in fig. S5). Both the abrupt-transition and functionally graded multimaterial tet data can be directly processed for multimaterial 3D printing [e.g., inkjet (38,42), grayscale digital light processing (DLP) (43,44), and multi-vat DLP (45,46)]; however, the functionally graded, multimaterial tet data are needed for continuously embedding the microstructural materials into the part.
We choose an image-based manufacturing approach (m-SLA) that enables us to easily embed the microstructures directly in the pixelated 3D printer slices and avoid a prohibitively expensive surface representation (STL) of the multimicrostructure-embedded part. We start by defining a set of transitional unit cells that provide a continuous shape morphing between each of the microstructural materials in our design (Fig. 3E). These transitional microstructural materials map to the graded regions of the functionally graded tet mesh according to a color mapping. We slice the functionally graded tet mesh to obtain color-coded "macro-slices" for each print layer. Using a consistent layerheight, we also slice each unit cell and tile the slices over the build area to obtain a stack of slices representing a one-unit-cell-high tessellation of the unit cell over the entire build area (a small region of the first layer of the "micro-slices" for each transitional unit cell is shown in Fig. 3E). Next, we use the colored macro-slices to map into the appropriate micro-slices and replace each pixel in the macro-slice with the corresponding pixel of the micro-slice. Essentially, we perform a Boolean intersection of the macrostructure with tessellations of each of the unit cells, according to the color mapping; however, we do it at the slice level rather than in the 3D geometry to avoid an expensive STL of the multimicrostructure-embedded part (additional details on the multimaterial slicing and multimicrostructure-embedding scheme are provided in Materials and Methods and in a flowchart in fig. S3). By comparing the "embedded slices" of the multimaterial tet mesh and the functionally graded tet mesh in Fig. 3 (C and D), it is clear that functional grading is critical to ensure connectivity of the microstructural materials. Macro-slices, micro-slices, embedded slices, and the macro-to-micro mapping for the other two, two-microstructural-material cantilever beams are provided in fig. S4.

Effect of porous, anisotropic microstructural materials in topology optimization
Next, we investigate how porous, anisotropic microstructural materials affect the geometry, topology, microstructure distribution, and structural efficiency of the topology-optimized parts. We design the same cantilever beam, but considering free selection from different subsets of seven porous, anisotropic microstructural materials. For comparison, we first design a reference beam considering a single, solid, isotropic material with domain volume fraction limited to ¯  = 0.022 (Fig. 4A). Then, using the same volume limit on the bulk material, we design the beam considering microstructural materials composed of the unit cells indicated to the right of each design in Fig. 4 (B to G). The volume constraint is active in all cases, and all designs in Fig. 4   convergence (within a tolerance of 4% after postprocessing). Normalized directional tensile and shear moduli, E 11 ′ / (E ˆ  i ) and G 12 ′ / (G ˆ  i ) , are plotted for the solid, isotropic material and the porous, anisotropic microstructural materials in Fig. 4 (H and I, respectively). To more easily understand the microstructural-material placement, the multimaterial tet meshes (before applying functional grading) are plotted for each design in Fig. 4, and the reported objective function values correspond to these multimaterial tet meshes with abrupt interfaces.
As in Fig. 1, the microstructures in Fig. 4 (B to G) distribute themselves according to the mechanics of a cantilever beam. That is, microstructures with higher tensile stiffness in the x 3 direction tend toward the tension/compression regions farthest from the neutral axis of the beam, and more isotropic microstructural materials that are also stiffer in shear (in the x 2 − x 3 plane) tend toward the inclined members and regions of high shear. In addition, not all of the available microstructural materials are used. For example, in Fig. 4 (E and F), only the most efficient in tension in the x 3 direction and the most isotropic and efficient in shear in the x 2 − x 3 plane are selected. The beam geometry is also influenced by the available microstructures. For example, the tension/compression members farthest from the neutral axis become increasingly inclined from the supports to the load point as the microstructural material becomes more isotropic. Notice that the design in Fig. 4B avoids this inclination as much as possible and the inclination of these members gradually increases as the microstructural material becomes more isotropic from microstructure 3 in Fig. 4 (E and G) to microstructure 5 in  (47) and can provide other biomimetic functionalities, e.g., buoyancy and impact resistance (see additional discussion in section S3).

Scaling to larger build volumes
No printer is able to print structures at unlimited size; thus, scaling is an important practical consideration. To demonstrate scalability of the proposed multimicrostructure-embedding scheme, we design and manufacture two additional structures at a larger scale: a hyperbolic paraboloid canopy structure and an Eiffel Tower-inspired structure. Although we are limited by the printer's display area (2560 × 1440 pixels) and build height (15 cm), we slice the parts for a larger build volume and then divide the slices into smaller images that the printer can handle (i.e., the structures are printed in pieces and assembled). To slice for a larger build volume, we simply increase the number of pixels per slice (for increased length and width) and generate a larger number of slices (for increased height).
Since we process small subsets of the pixels at a time (in parallel), the increase in number of pixels does not lead to increased memory requirements but does lead to increased slice time.
The hyperbolic paraboloid canopy structure is subjected to a uniformly distributed, vertical load on the top surface of the canopy, which is defined on the domain x 1 , x 2 ∈ [ − 0.5,0.5] by the equation We demonstrate various ways of representing multimaterial topology optimization data for 3D printing and the proposed procedure for continuously embedding multiple microstructures into topology-optimized parts. In (A), the data are projected directly onto the underlying hex mesh, leading to stairstepping features at the boundaries and material interfaces. In (B), a separate isosurface is generated for each material individually, leading to disjointed interfaces. In (C) and (D), the data are projected onto a tet mesh generated inside an isosurface of the composite data. By applying a convolution filter to the multimaterial tet data, we avoid the abrupt material interfaces shown in (C) and obtain the functionally graded interfaces shown in (D).
In (E), we define a number of transitional unit cells that enable us to move smoothly between the two microstructural materials (8 of the 25 transitional unit cells are shown here). We slice each unit cell, tile the slices over the build area, and assign a color to each set of micro-slices. The colors are used to map the macro-slices to the embedded slices shown in (C) and (D). The embedded slices in (C) are disconnected, while those in (D) are well-connected because of the transitional unit cells mapped to the functionally graded regions.
face-x microstructural material (shown in blue in Fig. 5), i.e., it does not participate in the optimization. A short tube just above the structure's fixed support is another passive region occupied by a solid, isotropic material (shown in red in Fig. 5). No material can occupy the inner volume of the tube or the space above the canopy. In the remainder of the domain, we use the proposed multimicrostructure topology optimization formulation to design a structure composed of octet and truncated octahedron unit cells (shown in cyan and yellow in Fig. 5) to transfer the loads from the canopy to the fixed support. The total optimizable domain volume fraction is limited to ¯  = 0.0096 . The topology optimization problem is described in more detail in fig. S7.
The final design and multimicrostructure-embedded, manufactured canopy structure are shown in Fig. 5. Details of the transition regions between microstructural materials are provided in Fig. 6. The manufactured part fits inside a bounding box of dimensions 11.6 cm by 11.6 cm by 14.4 cm, and the embedded unit cells are scaled to have an edge length of 2 mm (at this scale, the bar diameters are

of 13
300 m for the octet unit cells and 400 m for the face-x and truncated octahedron unit cells). Because the part exceeds the printer's build volume, we generate 2560 × 2560 pixel slices and print half of the structure at a time. Additional images of the canopy structure, including support structures needed during printing and the printed parts before assembly, are provided in fig. S8. Movies S1 and S2 illustrate the full design and manufacturing process for the canopy structure. Inspired by Gustave Eiffel's open-lattice, multiscale Eiffel Tower, we also use the multimicrostructure topology optimization formulation to design an Eiffel Tower-inspired structure. The domain and boundary conditions provided in fig. S7C are roughly based on those of the actual Eiffel Tower. Three floors are defined between the base and the top of the tower. The width of each floor reduces to imitate the shape of the actual tower, which was chosen by Eiffel to efficiently resist wind loading (48,49). Here, we adopt the shape of the original tower and design the remaining form, considering only vertical (gravity) loads. The structure is fully fixed at the corner regions of its base. Uniformly distributed, vertical loads with total force equal to 1, 0.766, and 0.3 are applied at floors 1, 2, and 3, respectively, and a point load of magnitude 0.01 is applied at the top of the tower. The tower is designed considering face-x, octet, and truncated octahedron microstructures with total domain volume fraction limited to ¯  = 0.008 . The final design and multimicrostructure-embedded, manufactured Eiffel Tower-inspired structure are shown in Fig. 7 (A and B,  respectively). The manufactured part fits inside a bounding box of dimensions 8.2 cm by 8.2 cm by 26.0 cm, and the embedded unit cells have the same scale as that used for the canopy structure. Because the part exceeds the printer's build volume, we generate 2560 × 2560 pixel slices and print the bottom portion of the tower in two pieces and the top portion of the tower as a third piece, as shown in Fig. 7C. Additional images of the Eiffel Tower-inspired structure are provided in fig. S9, and a comparison of the computer and physical model dimensions for the Eiffel Tower-inspired structure (and all other models) is provided in table S1.

DISCUSSION
Aiming toward convergence of design and manufacturing, we integrate porous, anisotropic microstructural materials into a general multimaterial topology optimization formulation and establish a procedure to (i) translate the multimaterial topology optimization data to a 3D printable, well-connected part with functionally graded interfaces and (ii) continuously embed multiple microstructures into the functionally graded part by mapping slices of a set of continuously varying microstructures into slices of the macrostructure. Surface representations typically used to communicate with a 3D printer cannot adequately represent microstructures at the resolution of the 3D printer within a complex macrostructure geometry, as demonstrated by our approach. Furthermore, most existing structures composed of architected materials contain a single microstructure type and cannot attain spatially varying mechanical properties with wellconnected interfaces, as demonstrated by our approach. Thus, the ideas presented here enable a new class of multiscale optimized structures that enhance our ability to mimic nature.
The presented multiscale structures are designed to maximize stiffness, and we show that by including additional microstructural materials that increase the design space of directional stiffness at the microstructural material level, the global stiffness of the macrostructure

of 13
tends to increase. Although strength criteria are not considered in design here, when compared to discrete interfaces, we anticipate that the functionally graded interfaces between microstructural materials will redistribute stress concentrations and mitigate detrimental effects of interfaces on strength. Strength criteria can be considered using topology optimization with local stress constraints (40). Defects introduced by the manufacturing process (e.g., imperfect and nonuniform microstructure nodal connectivity and strut cross sections) are also not considered here, but are expected to influence the mechanical performance of the additively manufactured parts, especially as the minimum feature size of the microstructures approaches the resolution of the 3D printer (50,51).
Although the proposed multimaterial topology optimization formulation can handle a sufficiently general class of porous, anisotropic microstructural materials (i.e., those for which an elasticity tensor can be provided), the current multimicrostructure-embedding scheme requires a few simplifying assumptions: (i) Each microstructural-material must be periodic. (ii) The same translation operations must be used to create each periodic microstructural material from its associated unit cell. (iii) A set of transitional unit cells must be defined to ensure a smooth, well-connected transition between the different unit cell geometries. The transitional unit cells can be devised using the intuitive approaches proposed here or other methods that may have better control over material property transitions (41,52). In addition to the lattice-like unit cells considered here, unit cells consisting of plate elements or triply periodic minimal surfaces (53) can be directly integrated with our approach. In addition to cubic materials, periodic materials without cubic symmetry [see examples by Zok et al. (54)] can be handled. Nonperiodic materials [e.g., spinodal architectures (55-58)] will be the subject of future work. In characterizing microstructural materials using homogenized properties, we assume infinite periodicity and separation of length scales (59), neither of which can be verified in the manufactured parts. Two factors prevent infinite periodicity: the presence of multiple microstructures and truncation of the periodic microstructure tessellations at the structure boundaries. The first factor can be mitigated by avoiding abrupt transitions between the microstructural materials, which is facilitated by the proposed continuous multimicrostructureembedding scheme. The boundary remains a challenge. Our parts contain truncated unit cells at the boundaries that likely influence the mechanical behavior. These edge effects cannot be completely removed, but making the unit cells conform to the boundary (35,37) can alleviate them.
Inadequate separation of length scales is due to the fact that we cannot print the microstructure at an infinitely small scale or the macrostructure at an infinitely large scale. However, the approach pursued here is scalable; that is, the maximum macrostructure size and minimum microstructure feature size are dictated by the 3D printer and not the data representation. Thus, our approach provides a means to obtain a practical separation of length scales, e.g., by using large-area projection microstereolithography (29) or high-area rapid printing (60) or by assembling the part from a number of components manufactured at a practical scale. The latter approach, in combination with our multimicrostructure-embedding scheme, could make optimized, multiscale, architectural engineering-scale structures possible. Furthermore, the scalability and modularity of the proposed scheme facilitates extension to other existing and yet to be invented additive manufacturing technologies that will enable the method to be explored in ways not yet anticipated.

Problem setting and optimization formulation
The multiscale topology optimization formulation is stated as a multimaterial problem for volume-constrained compliance minimization of an elastostatic body that accommodates many candidate materials and many local or global volume constraints (61)(62)(63) min In Eq. 1, we define a density field, Z = { z ℓ1 , … , z ℓm } ℓ=1 N , where z ℓi is a density design variable for each of the m candidate materials at the centroid of each of the N elements used to discretize the design domain, . The elemental density field, Y = { y ℓ1 , … , y ℓm } ℓ=1 N , is obtained as y i = Pz i , where y i and z i are the ith column of Y and Z, respectively, and P is a regularization map (density filter) that enforces well-posedness of the problem and a minimum length scale (64,65). The coefficients of P are defined as where ‖x i − x j ‖ 2 is the Euclidean norm between the centroids of elements i and j, R is a filter radius, and q defines the order of the filter (e.g., linear filter when q = 1) (66). We specify j = 1, …, K volume constraints that control any subset of the candidate materials in any subregion of the domain. Hence, j and E j represent the set of material and element indices associated with constraint j, respectively. Furthermore, A represents the volume of element ℓ. V = {  ℓ1 , … ,  ℓm } ℓ=1 N is the material volume for each of the m candidate materials in each of the N elements, where v ℓi = m V ( y ℓi ) = y ℓi ˆ  i accounts for the material porosity via the unit cell volume, ˆ  i , of material i, and ¯  j is the volume fraction limit for constraint j. When the subscript, j, is omitted, it is understood that there is only one volume constraint.
The same discretization used for the optimization problem is also used to solve for the displacement field, U, via the discretized state equations of static elasticity, (K + I)U = F, which are derived from the principle of minimum potential energy with a Tikhonov regularization term,  _ 2 U T U , included in the expression of total potential energy. The Tikhonov regularization term prevents the stiffness matrix from becoming singular because of low-density regions of the domain (67). In the state equations, the stiffness matrix, K = ∑ ℓ=1 N k ℓ , is assembled from the element stiffness matrices, F i = ∫ ~  N t · N i ds is the vector of design-independent nodal loads, t is the traction applied on the portion, ~  N , of the domain boundary, and N is the vector of interpolation (shape) functions used to interpolate quantities between the mesh nodal points.
The formulation in Eq. 1 is not specific to any type of material. Instead, the material properties are embedded in the stiffness matrix, which is a function of the penalized element densities, W = { w ℓ1 , … , w ℓm } ℓ=1 N , where w ℓi = m W ( y ℓi ) = y ℓi p , with p > 1 [SIMP (14,15)] to penalize intermediate densities. Then, the element stiffness matrix of element is obtained via a material interpolation function, m M , that penalizes material mixing (63,68) where w ′ is the th row of W, ( k ℓi is the element stiffness matrix for material i in element ,  is the domain of element ℓ, B is the strain-displacement matrix of shape function derivatives, and D i H is the (homogenized) elasticity matrix characterizing material i, which is supplied as input for each of the candidate materials in Eq. 1. In addition, the parameter, 0 <  < 1, controls the amount of allowable mixing. In general, we seek solutions without material mixing (i.e.,  = 1); however, the problem in Eq. 1 is convex for  = 0 and p = 1, and a continuation scheme on these parameters can be used to bias the solution toward the convex one at the beginning of the optimization iterations (63).

Gradient-based solution scheme
To solve Eq. 1, we adopt a gradient-based approach in which we use derivatives of the objective and constraint functions with respect to the design variables to iteratively guide the design toward an optimal solution. We adopt the Zhang-Paulino-Ramos Jr. (ZPR) design variable update scheme (61), which uses Lagrangian duality to solve a series of convex approximate subproblems of Eq. 1 around the current design, Z 0 . The ZPR update scheme was derived specifically for the multimaterial formulation of interest in Eq. 1, in which each design variable is associated with a single volume constraint. Because of separability, the constraints can be updated independently, leading to an update scheme that efficiently handles a large number of volume constraints at a cost on par with that of the optimality criteria update scheme (69).
The original ZPR update scheme was derived for a monotonically decreasing objective function. Since the derivatives of the objective function in Eq. 1 may become positive in regions of material mixing (62,63), we integrate sensitivity separation into the ZPR update scheme (70,71) by decomposing the objective function gradient into positive and negative components to arrive at the following nonmonotonic, convex approximation of the objective function The approximation in Eq. 4 is the sum of a constant term that can be neglected for optimization, a monotonically decreasing convex function in exponential intermediate variables, ξ ℓi ( z ℓi ) = z ℓi −α , α > 0 , and a monotonically (linearly) increasing function. For decomposition of the objective gradient, we adopt the scheme proposed in (70), which uses approximate second-order information to account for the curvature of the objective function such that the negative and positive components, respectively, are where h ℓi 0 is a Broyden-Fletcher-Goldfarb-Shanno (72) approximation of the diagonal terms of the objective function's Hessian matrix.
For brevity, we omit the full derivation of the ZPR update scheme considering sensitivity separation with approximate second-order information but refer the reader to (70), where they arrive at the following update with z ℓi new as the design at the next iteration and z ℓi * as the candidate design for the next iteration that is accepted if it is within bounds defined by box constraints, ¯  and  _ , and move limit, M. The candidate design, z ℓi * , is obtained from a fixed-point iteration of the form where and Eq. 8 includes the heuristic ZPR filter introduced in (62, 63). The sensitivities needed in Eqs. 5 and 9 are found using the chain rule where ∂y i /∂z i = P T and the other components are as follows To compute Eq. 11, we also need the following derivative The iterative solution scheme is said to have converged when either the maximum number of iterations is reached or the change in the design is small

Topology optimization algorithmic parameters and computational resources
The topology optimization results were obtained using the following algorithmic parameters: ZPR move limit, M = 0.15; ZPR intermediate variable exponent,  = 1; and box constraints,  _ = 0 and ¯  = 1 . To bias the solution toward the convex one, we perform five continuation steps on the material interpolation parameters: p = [1,1.5,2,2.5,3] and  = [0,0.2,0.5,0.8,1]. Each continuation step is said to converge after reaching the maximum number of iterations, MaxIter = [100,100,100,100,200], or the convergence tolerance, tol = 0.01. In each problem, the initial guess is specified such that the volume fraction limit of each constraint is evenly distributed among the design variables associated with each available material. For example, in the seven-microstructural-material design in Fig. 4G, the initial value of each design variable is ¯ v / 7 (the effect of the initial guess is studied in section S4). For the cantilever problems, the filter exponent and radius are q = 3 and R = 0.064, respectively, and the problem is solved on half of the domain (symmetry enforced on the x 1 − x 3 plane) on a 3 × 0.625 × 1 hex mesh with 192 × 40 × 64 elements, for a total of 491,520 elements. For the canopy problem, the filter exponent and radius are q = 3 and R = 0.032, respectively, and the problem is solved on a 1 × 1 × 1.5 hex mesh with 80 × 80 × 120 elements; however, elements with centroid above the hyperbolic paraboloid surface and interior to the tube are removed, for a total of 510,411 elements. For the Eiffel Tower-inspired problem, the filter exponent and radius are q = 3 and R = 0.025, respectively, and the problem is solved on a quarter of the domain (symmetry enforced on the x 1 − x 3 and x 2 − x 3 planes) on a 0.5 × 0.5 × 2.75 hex mesh with 46 × 46 × 253 elements, for a total of 535,381 elements. All problems were run using MATLAB 2019a on an Intel(R) Xeon(R) central processing unit (CPU) ES-1660 v3 @ 3.0 GHz with 256-GB random-access memory (RAM) and NVIDIA Quadro M5000 graphics processing unit (GPU) or MATLAB 2018b on a Dual Intel(R) Xeon(R) Silver 4116 CPU@2.10 GHz with 256-GB RAM and NVIDIA Quadro P1000 GPU.

Homogenized properties of porous, anisotropic microstructural materials
The formulation in Eq. 1 can handle any material for which the full material tensor is available. For simplicity, we consider microstructural materials composed of a periodic tessellation of lattice-based unit cells defined on the unit cube, where the lattice elements are cylindrical bars. Effective macroscopic properties are obtained using computational homogenization (59,73), specifically using an educational MATLAB code (74). The geometry of the lattice unit cell is inscribed in a hex mesh that is used in the homogenization computations. Borrowing ideas from the educational polygonal mesh generator, PolyMesher (75), we use signed distance functions to compute the signed distance of each hex centroid from the boundary of the cylinders. Any hex element with a negative signed distance to one of the cylinders' boundaries is determined to be inside the unit cell structure and is assigned a value of one. All other hex elements are void and are assigned a value of zero. Such implementation facilitates extension to other types of unit cells (e.g., unit cells composed of noncylindrical bars or plates). The educational homogenization code (74) outputs the homogenized stiffness elasticity tensor of microstructural material i in matrix (Voigt) notation, D i H , and we can easily compute the volume fraction of microstructural material i's unit cell, ˆ  i , as the sum of the solid hex element volumes. These two properties, D i H and ˆ  i , are needed for each candidate microstructural material defined in Eq. 1. The directional tensile and shear moduli are extracted from D i H after performing a coordinate transformation [see section S1 and references (76, 77) for more details].
In all cases studied here, the bulk material has Young's modulus, E = 1, and Poisson's ratio,  = 0.3. In addition, the computational homogenization is performed using a hex mesh with at least 160 × 160 × 160 elements.

Multimaterial slicing and continuous multimicrostructure embedding
Although the multimaterial interpolation in Eq. 3 prevents material mixing, small mixing regions occur at the material interfaces as an artifact of the density filter used in topology optimization. Furthermore, in some elements that contain mixing, the total density of material may sum to greater than one (62). Thus, we first remove mixing from the converged topology-optimized result and limit the total element densities to one by assigning At this stage, we maintain intermediate densities at the structure boundaries and at the material interfaces that also result from the density filter. Using an appropriate isovalue (i.e., one that preserves the volume fraction specified in topology optimization to within a tolerance of 4%), we generate an isosurface of the composite density field, which is defined as y ℓ c = ∑ i=1 m y ℓi , ℓ = 1, … , N , and generate a tet mesh within, using the iso2mesh toolbox (78), which relies on TetGen (79). Making use of the disjointed isosurfaces generated for each material individually (see Fig. 3B), we use the inpolyhedron function from MATLAB's file exchange to determine which tet centroids fall within each of the material isosurfaces and assign a material accordingly. A relatively small number of tets near the material interfaces that do not fall within any of the disjointed isosurfaces (i.e., interface tets) are assigned a material according to the nearest density of the postprocessed topology optimization data. The functionally graded tet mesh is then obtained by applying the filter in Eq. 2 to tets within radius, R, of the interface tets ( fig. S5 shows the effect of the magnitude of R on the length scale of the graded region). For the cantilever problems, the filter power and radius used to generate the functionally graded tet mesh are q = 1 and R = 0.10, respectively, for the canopy problem, q = 1 and R = 0.04, and for Eiffel Tower-inspired problem, q = 1 and R = 0.08, where the filter radius, R, is relative to the dimensionless domain dimensions provided in Fig. 1 and fig. S7.
Next, we slice the multimaterial tet mesh by thinking of the printer's build volume as a 3D voxel matrix, where the voxel dimensions are controlled by the printer's pixel dimensions and the layerheight. From the functionally graded tet mesh, we generate a scattered interpolant, F(x 1 , x 2 , x 3 ), that takes in a voxel's 3D coordinate and outputs a nearest-neighbor interpolated material. The 3D printer has more than 11 billion voxels, but using the postprocessed topology optimization data on the underlying hex mesh, we reduce the number of voxels that need to be evaluated. First, all voxels that fall outside the hex mesh are known to be void. In addition, all voxels that fall within a hex element with density equal to one can be immediately assigned the corresponding material. Thus, we only need to determine a material for the subset of voxels that fall within grayscale hex elements, i.e., those near the structure boundaries and at the material interfaces. From this subset of voxels, we find those that fall inside of the composite isosurface (using inpolyhedron), evaluate them using the scattered interpolant, F(x 1 , x 2 , x 3 ), and assign a material accordingly. Memory requirements are controlled by evaluating small subsets of the voxel matrix in parallel.
Using the multimaterial voxel matrix, we can easily obtain pixelated slices for each layer, using color to represent the functionally graded materials. Next, we generate the mapping between color and microstructure and embed the microstructures into the pixelated slices as described in Results. Last, the images are converted to black and white and sent to the 3D printer in png format. A flowchart summarizing the overall process from design to manufacturing is provided in fig. S3. All computations needed to translate the multimaterial topology optimization data to microstructure-embedded 3D printer slices were done using MATLAB 2019a on an Intel(R) Xeon(R) CPU ES-1660 v3 @ 3.0 GHz with 256-GB RAM. m-SLA manufacturing All physical models were fabricated using the Original Prusa SL1 m-SLA 3D printer (Prusa Research, Czech Republic), which shines ultraviolet light onto the underside of a resin vat, masked by a 2560 × 1440 pixel liquid crystal display according to pixelated images (slices), to cure the part layer by layer. The pixel edge length is 47.25 m, and we print with a 50-m layer height. The build volume is 120.96 mm by 68.04 mm by 150 mm. All models are built using Prusa's Transparent Red Tough Resin with 6-s exposure time per layer. Slicing and multimicrostructure embedding are done with the in-house MATLAB code described previously, and the generated black-and-white png images for each layer are provided to the 3D printer. Support material was not required for any of the cantilever designs or the Eiffel Towerinspired structure. The canopy required a support structure, which was designed using Rhino software (see fig. S8). Because of build volume restrictions, the canopy and Eiffel Tower-inspired structure were printed in two parts and three parts, respectively, and bonded together using Krazy Glue (see photos in Fig. 7 and fig. S8).