We present a new method for the simultaneous topology optimization and material selection of structures made by the union of discrete geometric components, where each component is made of one of multiple available materials. Our approach is based on the geometry projection method, whereby an analytical description of the geometric components is smoothly mapped onto a density field on a fixed analysis grid. In addition to the parameters that dictate the dimensions, position, and orientation of the component, a size variable per available material is ascribed to each component. A size variable value of unity indicates that the component is made of the corresponding material. Moreover, all size variables can be zero, signifying the component is entirely removed from the design. We penalize intermediate values of the size variables via an aggregate constraint in the optimization. We also introduce a mutual material exclusion constraint that ensures that at most one material has a unity size variable in each geometric component. In addition to these constraints, we propose a novel aggregation scheme to perform the union of geometric components with dissimilar materials. These ingredients facilitate treatment of the multi-material case. Our formulation can be readily extended to any number of materials. We demonstrate our method with several numerical examples.

## Introduction

Topology optimization has been used extensively in order to generate novel designs that improve structural performance while decreasing the cost. One important option to improve the structure is to employ multiple materials. By having an appropriate distribution of materials, multimaterial designs can outperform designs made of each of the materials separately. Moreover, besides decreasing cost, multimaterial designs can take advantage of substantial differences in properties to perform multiple functions.

Methods for multimaterial topology optimization were first introduced in density-based approaches. In [1], a method is developed to determine the optimal distribution of multiple phases to obtain composite materials with extreme thermal expansion behavior. This work employs an extension of the power-law interpolation used in solid isotropic material penalization (SIMP) [2,3] to three-phase material designs (two solid materials and void). The interpolation uses two design variables, one that indicates where to put material or void and another that determines which material to choose. This and similar formulations are applied to various problems, such as design of multiphase composites with extremal bulk modulus [4], design of multiphase piezoelectric actuators [5,6], combined optimization of material and voltage distribution [7], and optimal reinforcement of concrete structures [8]. As indicated in Ref. [9], this formulation violates the Hashin–Shtrikman bounds and renders different designs if the phases are interchanged. In Ref. [9], a formulation that combines the power-law penalization with an interpolation of every material property within its Hashin–Shtrikman bounds is proposed to circumvent the foregoing limitations to design multimaterial actuators. In all of these applications, the maximum number of phases involved is three (including two solid phases and a void phase) as the generalization of their interpolation schemes to more than three phases becomes quite involved.

A method for topology optimization of multi-material compliant mechanisms using an alternative material interpolation scheme that employs only one density variable is proposed in Ref. [10]. As opposed to the single-material topology optimization, where the optimal density indicates either void (zero) or solid (unity), in this method, the density variable in the optimal design can take any real value. The proximity of this value to one of the means of a sum of Gaussian distributions with given mean locations and with modes corresponding to the property values for different materials indicates the choice of material. To avoid premature convergence to undesired local minima, this method starts with a large standard deviation for each material's property distribution and decreases it gradually in order to sharpen the peaks.

The discrete material optimization method (DMO) [11] simultaneously optimizes the stacking sequence, reinforcement orientation, and choice of material of composite shell structures. In this method, there is a set of materials from which each finite element can be made in order to minimize the objective function. A material interpolation scheme is formulated, whereby an increase in the weighting factor for one material decreases the weighting factors for all other materials. The weighting factor for a given material indicates the relative influence of that material on the properties of the corresponding finite element. If the weighting factor for a given material is unity, it indicates that finite element is made solely of that material (e.g., pure phases). One challenge with this formulation is that the weighting factors do not add up to unity except for pure phases. The proposed solution is to normalize each weighting factor by the summation of all weighting factors. However, this alters the penalization effect and the monotonic convergence. The DMO method is used in Ref. [12] to optimize the buckling behavior of multimaterial composite shell structures, and in Ref. [13], to compare topology optimization of multimaterial structures with a mass constraint using both this method and the one presented in Ref. [1]. The work in Ref. [14] formulates simpler multimaterial interpolations by imposing one linear constraint per finite element that ensures that it is made of only one material (or no material).

Other density-based methods employ a heuristic rule for the multimaterial interpolation. The work in Ref. [15] uses a homogenization approach with a unit cell made of two materials and a square hole. The optimization allows for continuous mixtures of the three phases. Then, a postprocessing heuristic rule is applied to every finite element in the optimal design to determine if it is made of one of the two materials or void. In Ref. [16], material is iteratively removed from a fully solid initial design by applying a rule to each finite element that changes its material, or that makes it solid or void based on the elemental change in compliance incurred by these changes.

There also exist several level set methods for topology optimization of multimaterial structures. In the color level-set method [17], the boundaries of multiple regions made of different materials are given by the zero level-sets of multiple functions. Instead of using a level set function per material, this method employs *m* level-set functions whose combinations can represent up to $2m$ materials. This strategy circumvents the need for an equality constraint that would be required to avoid overlapping materials if one level set per material was used. Since the combination of level sets renders one and only one material at any point, this method eliminates the need for material interpolation. This method is applied to topology optimization of multimaterial compliant mechanisms in Ref. [18] and stress-based topology optimization of continuum structures involving multiple materials in Ref. [19]. In Ref. [20], the color level-set method is combined with a variational approach for optimization of structures made of functionally graded materials.

In the level set method of Ref. [21], only one level-set, piece-wise constant function is used as an “index” that indicates the choice of material. A constraint is applied to ensure that the level set function converges to the index values. This method is applied to design of piezoelectric actuators [22]. In Ref. [23], *n* level-set functions are used to represent *n* solid materials and void (i.e., *n* + 1 phases). This method employs a material interpolation similar to the one proposed in Ref. [1], except it uses the Heaviside of the level-set functions instead of density variables.

The design of multimaterial structures has also been studied with phase-field topology optimization methods. The generalized SIMP interpolation scheme proposed in Ref. [11] in combination with the mutual material exclusion constraint of Ref. [14] is employed in the phase field methods of Refs. [24–26]. In Ref. [27], the volume fraction of each phase in phase field model directly represents the contribution of the corresponding material to the material properties. A constraint on these volume fractions ensures that they add up to unity, while a penalization term added to the objective function ensures that the volume fractions are zero or unity upon convergence. Phase field methods for topology optimization of multimaterial structures require a very large number of iterations, typically in the order of thousands.

All the aforementioned methods produce organic designs that cannot be readily manufactured with stock material. The geometry projection method presented in Refs. [28–30] generates designs that can be made of stock material, such as bars and plates, by smoothly mapping an analytical description of the geometric elements onto a density field over a fixed finite element grid. A size variable is ascribed to each geometric component that is penalized in the spirit of SIMP, allowing the optimizer to entirely remove that component from the design. In the moving morphable components method of Refs. [31] and [32], a parametric description of bars is mapped onto a level set representation (termed the “topological description function”) for the analysis. The foregoing methods for design with discrete geometric elements use a single material. Recently, the work in Ref. [33] proposed a geometry projection method for the topology optimization of multimaterial, three-dimensional (3D) lattice structures. This method adapts the material interpolation scheme of Ref. [1] to the geometry projection framework, by assigning additional variables to each geometric component that determine the choice of material. Although possible, the extension of this interpolation scheme to more than two materials is not straightforward and it is more prone to getting locked into undesired local minima, as reported in Ref. [11]. Also recently, the work in Ref. [34] extends the moving morphable components method to accommodate multiple materials by seeding the initial design with geometric components made of different materials (i.e., the material of a component does not change during the optimization). The material interpolation scheme in regions where two or more geometric components intersect chooses the stiffer material of any intersecting component.

In this paper, we extend the geometry projection method so that it can be readily applied to structures whose geometric components can be made of one of any number of available materials. We achieve this by means of several key ingredients. The first ingredient corresponds to a new aggregation function used to perform the union of geometric components, since the function used in the previous geometry projection schemes cannot accommodate components made of dissimilar materials. A second important ingredient pertains to the interpolation of properties from the multiple material candidates. Here, we adapt the DMO formulation to accommodate discrete geometric components made of different materials. Unlike DMO, however, we impose two separate aggregate constraints in the optimization to (a) penalize the size variables of geometric components (discreteness constraint), and (b) to ensure that each component is made of at most one of the available materials (mutual material exclusion constraint). We demonstrate the effectiveness of our method via numerical examples.

The rest of the paper is organized as follows. Section 2 describes the projection of geometric components onto the analysis mesh for analysis, including the aggregation function for the union of geometric components and the interpolation of material properties. Sections 3 and 4 detail the constraints to penalize intermediate values of size variables and to ensure that geometric components are made of at most one material, respectively. The modification of the geometry projection to enforce a symmetric design is presented in Sec. 5. Section 6 describes the optimization problem. We present numerical examples to demonstrate our method in Sec. 7 and draw conclusions of our work in Sec. 8.

## Geometry Projection

To perform the analysis for any given design made of the union of geometric components, we use the geometry projection method [28,29], wherein a parametric description of the components is smoothly mapped onto a density field over a fixed grid. In this work, we model bars as the solid enclosed by the surface defined by the set of points equidistant to a straight line segment (the medial axis of the bar). In two-dimensional (2D), these bars correspond to semicircular-ended rectangles of width *w* and out-of-plane thickness *t* [29], and in three-dimensional, they correspond to cylinders with semi-spherical ends, both with diameter *w*. Here, we follow closely [29] to define the projected density for a single component; the reader is referred to that work for more details.

**p**in the design domain as the area fraction (for 2D components) or volume fraction (for 3D components) of the portion of the sample window $Bpr$ that intersects the solid structure

*ω*(cf. Fig. 1), i.e.,

*q*can be approximated for 2D components as the area ratio of the circular segment of height $r\u2212\varphi q$ (cf. Fig. 1)

where $\varphi q(p)$ is the signed distance from point $p$ to the boundary $\u2202\omega q$ of geometric component *q*. That is, the projected density at a point $p$ with respect to a component *q* is entirely dictated by its signed distance to the component (for a fixed sampling window radius *r*).

*q*is parameterized by the positions of the endpoints of its medial axis, $xqo$ and $xqf$. Thus, the signed distance is given by

*d*is the unsigned distance from point $p$ to the medial axis of bar

_{q}*q*. As indicated by the arguments in the equation above, the minimum distance

*d*is a function of the bar parameters (cf. Fig. 1) and can be computed as

_{q}In these expressions, $\Vert \xb7\Vert $ denotes the Euclidean 2-norm, ⊗ denotes the dyadic (or tensor) product, and $Pa\u22a5$ is the perpendicular projector on $a$. We fix *w* as we desire to design structures with fixed-width members, but note that this could readily be a design parameter.

Since several components can intersect, in the previous work we aggregated the corresponding densities in the intersection by using a *p*-norm, which smoothly approximates the maximum projected density of any of the intersecting components. While this approach is effective for single-material structures, it does not accommodate components made of multiple materials. In this work, we consider a different aggregation strategy, which we explain in the following.

*q*to the projected density at point $p$ to be 1.0 if $p$ is in the interior of the bar, and 0.0 if it is outside the bar. We can achieve this by using the Heaviside function of the signed distance to express the effective density as

*q*.

*N*is the number of bars in the design. With the Heaviside weighting factor, a component contributes to the effective density when $\varphi q\u22640$, i.e., if $p$ is inside the component. The sum of all weighting factors in the denominator ensures $\rho eff\u2208[0,1]$. Since the exact Heaviside function is not differentiable, we replace it with a smooth approximation $H\u0303\epsilon $ [35] so that design sensitivities are well defined and we can use gradient-based optimizers

_{b}We introduce the exponent *p* > 1 as a parameter to sharpen the Heaviside (i.e., to reduce the range over which it attains an intermediate value). In the case of single-material structures, the material properties are modified by some function of the foregoing effective density, as in our previous work. With the smooth Heaviside, it is clear that material properties vary continuously within a narrow band around the boundary of the bars.

*N*materials, we express the effective elasticity tensor at point $p$ as

_{m}*i*, $\u2102min$ is the elasticity tensor of a weak isotropic material that prevents the analysis from being ill-posed, and $\rho effi$ is an

*effective density per material*given by

*q*and material

*i*. Here, we adapt the DMO interpolation scheme to our method, and define these weights as

The main difference between this expression and the original DMO approach for density-based topology optimization is in the variables *α*. In the density-based approach, these variables correspond to element-wise densities per material. In our approach, we ascribe a size variable $\alpha iq\u2208{0,1}$ to geometric component *q* that indicates if it is made of material *i* when $\alpha iq=1$. Correspondingly, the vector of design variables for bar *q* is now given by $zq=[xq0T\u2009xqfT\u2009\alpha qT]T$ (note that $wiq$ only depends on $\alpha q$ and not on $xq0T$ or $xqfT$, but we use the foregoing notation for simplicity).

These size variables are discrete, which precludes the use of nonlinear programming methods; therefore, we relax them so that they can take any value between 0.0 and 1.0. To ensure that the optimizer converges to a design with pure phases, we penalize intermediate values in the same spirit as the SIMP method. The second difference between our approach and the original DMO method lies in that the penalization is not imposed through the weights of Eq. (16), but through an optimization constraint, as described in Sec. 3. Furthermore, as we desire each component to be made of one and only one material, a mechanism to ensure that $\alpha iq$ is 1.0 for at most one material is required. We detail a constraint to enforce this requirement in Sec. 4.

*i*(i.e., with $\alpha iq=1$ and $\alpha jq=0$ for $j\u2260i$) and a void bar (i.e., with $\alpha iq=0\u2009\u2200i$), the effective density for material

*i*is 0.5. Consequently, the effective elasticity tensor of Eq. (14) equals $0.5(\u2102i+\u2102min)$, which is clearly incorrect. This occurs because the denominator is counting

*all*the bars in the intersection and not only the solid bars. To remedy this, we modify the denominator so that it only counts the nonvoid bars (i.e., those for which $\u2211iNm\alpha iq>0$) as follows:

*B*equals 1 if all intersecting bars have zero size variables, and 0 if at least one size variable equals 1. We note that both situations occur only when all intersecting bars have pure phases. Finally, the maximum function is not differentiable, hence we replace it with a Kreisselmeier–Steinhauser (KS) approximation [36]

The KS function approximates better the maximum as the parameter *k* increases.

## Discreteness Constraint

*lower-bound*KS function

where $\alpha iq$ is the size variable for component *q* and material *i*. The lower-bound KS function approaches the maximum from below, therefore we circumvent the need to adaptively adjust the constraint limit (cf. [37]) and guarantee that $gd\u2208[0,1]$. As before, we denote *g _{d}* being dependent on $z$ for notational simplicity, however it only depends on the size variables. If all size variables are 0.5, the above constraint will attain its largest value, 1.0. If all size variables are either zero or unity, the constraint value will be zero.

*I*+ 1 is defined as

## Mutual Material Exclusion Constraint

*N*=

*N*.

_{b}*q*satisfy the discreteness constraint, the term in parenthesis equals 1.0 if the bar is solid (in which case it is made of one and only material) or 0.0 if the bar is void. We employ a similar continuation strategy as in Sec. 3 to avoid premature convergence to an undesired local minimum. To this end, we write

Note that $\epsilon m(0)$ must be smaller than 1.0, as otherwise it is possible for a bar to be made of more than one pure phase. In this situation, a decrease in one of the material size variables can increase the violation of the discreteness constraint, making it difficult for the optimizer to find feasible designs with bars made of at most one material.

## Symmetry

If the boundary conditions and design envelope are such that a symmetric design is expected, density-based and level set topology optimization methods have significant design freedom that allows them to readily produce symmetric designs. However, as discussed in Ref. [29], the more restrictive design representation enforced by discrete geometric components can lead to asymmetric designs as the optimizer satisfies exactly the resource constraint. Specifically, if the design can only be made of a finite number of discrete geometric components, it is entirely possible to find an asymmetric design that exactly satisfies the weight fraction constraint and it has lower compliance than the best symmetric design that can be obtained with the available set of components.

In cases where a symmetric design is expected or desired, we employ a simple approach to enforce symmetry that consists of reflecting the geometry projection for elements on the opposite side of the symmetry plane. This approach was introduced in Ref. [33]. We define geometric components only on one side of the symmetry plane and bound their location to ensure they remain on that side. To compute the projected density for a point on the reflected side, we reflect the point with respect to the symmetry plane, and then compute the projected density for the reflected point as usual using Eq. (3). The sensitivities are modified accordingly to account for this reflection.

where $s$ is a vector from the origin to the closest point on the line of symmetry, $\varphi $ is the angle between the symmetry line and the $e1$ axis, and $R$ is the reflection matrix. We use this approach to enforce symmetry in the example of Sec. 7.3.

## Optimization Problem

*ω*and Ω, respectively, with $\omega \u2282\Omega $. We consider linearly elastic problems without a body load. The compliance minimization problem is then stated as

*γ*is the physical density for material

_{i}*i*,

*γ*

_{ref}is a reference density that we choose to be 1 here, $|\Omega |$ denotes the volume of the design region, $v$ denotes the virtual displacement, $t$ denotes the design-independent traction, and

*w*is the weight fraction. $U\Omega :={u|u\u2208H1(\Omega ),u|\Gamma u=0}$ is the set of admissible displacements, and $u$ is the displacement obtained from the solution to Eq. (36), with $a$ and $l$ being the energy bilinear and load linear forms, respectively, given by

_{f}We stop the optimization when the discreteness and mutual material exclusion constraints are feasible, and when the relative change of the objective function between consecutive iterations falls below the specified value $\Delta f*$.

## Examples

To illustrate the effectiveness of our method, we present several examples. In all examples, we employ bilinear quadrilateral elements for the analysis. The entire code, including the finite element analysis, the sensitivities calculation, and the optimization, is implemented in matlab. To solve the optimization problem, we use the method of moving asymptotes [38,39] with the default parameters described in Ref. [39], i.e., $a0=1$ for the objective function, and *a _{l}* = 0,

*c*= 1000, and

_{l}*d*= 1 for every constraint

_{l}*l*in the optimization (we refer the reader to [39] for a description of these parameters). Unless noted, the stopping criterion on the relative change in compliance between consecutive iterations is $\Delta f*=10\u22124$. Also, unless specified, we employ a move limit of

*m*= 0.3. We use $\epsilon d0=1.0$ and $\epsilon d*=0.01$ for the discreteness constraint of Eq. (27), and $\epsilon m0=0.3$ and $\epsilon m0=0.01$ for the mutual material exclusion constraint of Eq. (31). We use a power of

*p*= 2 in the smooth Heaviside approximation of Eq. (13). The constant

*k*= 25 is used for the KS functions of Eqs. (23) and (25). All the materials considered are homogeneous, isotropic, and linearly elastic with Poisson's ratio $\nu =0.3$, but with different Young's moduli and material densities.

### Two-Bar Cantilever Beam.

The first example is a short cantilever beam made of two bar. The design envelope, boundary conditions, and initial design are shown in Fig. 2. The design envelope is meshed with a regular grid of 80 × 80 elements. There are two available materials with Young's moduli $E1=10$ and $E2=5$, and physical densities $\gamma 1=0.9$ and $\gamma 2=0.45$, respectively. The initial design, shown in the same figure, consists of two horizontal bars of width *w* = 0.25 and with $\alpha 1q=\alpha 2q=0.5$, *q* = 1, 2. For this problem, we use a looser stopping criterion on the relative change in compliance between consecutive iterations of $\Delta f*=10\u22123$.

We perform the optimization for several weight-fraction limits $wf*$, and the results are presented in Table 1. The choice of $wf*$ for each run corresponds to expected configurations. For example, $wf*=0.0520$ corresponds to a design made of a single horizontal bar made of material 2 and completely inside the design envelope. Some of the weight-fraction limits account for the fact that half of the horizontal bar may be outside of the design envelope.

Design ID | Optimal design | $wf*$ | w_{f} | C | $\alpha 11$ | $\alpha 21$ | $\alpha 12$ | $\alpha 22$ | Its. |
---|---|---|---|---|---|---|---|---|---|

1 | 0.0520 | 0.0506 | 1.0030 | 0.0000 | 1.0000 | 0.0010 | 0.0000 | 191 | |

2 | 0.1020 | 0.1019 | 0.0347 | 1.0000 | 0.0000 | 0.0009 | 0.9999 | 164 | |

3 | 0.1150 | 0.1150 | 0.0193 | 1.0000 | 0.0000 | 0.0008 | 0.9998 | 65 | |

4 | 0.1600 | 0.1598 | 0.0163 | 1.0000 | 0.0000 | 0.0000 | 1.0000 | 35 | |

5 | 0.1800 | 0.1800 | 0.0130 | 1.0000 | 0.0000 | 0.9996 | 0.0006 | 81 | |

6 | 0.2300 | 0.2258 | 0.0103 | 1.0000 | 0.0000 | 1.0000 | 0.0000 | 36 |

Design ID | Optimal design | $wf*$ | w_{f} | C | $\alpha 11$ | $\alpha 21$ | $\alpha 12$ | $\alpha 22$ | Its. |
---|---|---|---|---|---|---|---|---|---|

1 | 0.0520 | 0.0506 | 1.0030 | 0.0000 | 1.0000 | 0.0010 | 0.0000 | 191 | |

2 | 0.1020 | 0.1019 | 0.0347 | 1.0000 | 0.0000 | 0.0009 | 0.9999 | 164 | |

3 | 0.1150 | 0.1150 | 0.0193 | 1.0000 | 0.0000 | 0.0008 | 0.9998 | 65 | |

4 | 0.1600 | 0.1598 | 0.0163 | 1.0000 | 0.0000 | 0.0000 | 1.0000 | 35 | |

5 | 0.1800 | 0.1800 | 0.0130 | 1.0000 | 0.0000 | 0.9996 | 0.0006 | 81 | |

6 | 0.2300 | 0.2258 | 0.0103 | 1.0000 | 0.0000 | 1.0000 | 0.0000 | 36 |

Two interesting cases are worth noting. For the second run, the weight-fraction limit $wf*=0.1020$ corresponds to the weight fraction of a single horizontal bar made of material 1 and completely inside the design envelope. However, the optimization produces a better design by using a short diagonal bar made of material 2. For the third run, the weight-fraction limit $wf*=0.115$ corresponds to a V-shape design made of the lighter material. However, the optimizer finds a better design whereby the horizontal bar is partially outside the design envelope.

In Fig. 3, we plot the compliance versus the weight fraction for the optimal designs of all six runs to illustrate that the compliance decreases as the weight fraction increases. As expected, for the lowest weight-fraction limit, the optimal design corresponds to a single bar made of the lighter (and weaker) material. As the weight-fraction constraint increases, the optimizer first introduces a second bar made of the heavier (and stiffer) material and eventually obtains a two-bar design made of the stiffest material. Finally, Fig. 4 shows the effective density for each of the two materials for the initial and optimal designs of run 4 in Table 1.

### Messerschmitt–Bölkow–Blohm Beam.

The second example corresponds to the Messerschmitt–Bölkow–Blohm (MBB) beam widely studied in topology optimization. The design envelope, boundary conditions, and initial design are shown in Fig. 5. We note that, as we have discussed in the previous work (cf., [29]), geometry projection methods are more prone to converging to different local minima than free-form topology optimization methods due to the more restrictive design representation. The design envelope is meshed with a regular grid of 160 × 40 elements. We use the same two materials of the preceding example. The initial design consists of 21 bar of width *w* = 0.4 and with $\alpha 1q=\alpha 2q=0.5$, $q=1,\u2026,21$. We consider two configurations of bars. In one configuration, as in the previous example, the endpoints $xq0$ and $xqf$ for each bar are independent from other bars, so that bars are “floating” inside the design envelope. In the second configuration, bars share common endpoints, so that the design remains connected at all times.

We perform the optimization for several weight-fraction limits, $wf*=0.1,0.11,\u2026,0.19$. The optimal designs for the floating and connected configurations are shown in Figs. 6 and 7, respectively. Several of the designs obtained resemble known solutions for the MBB beam. As expected, the runs with floating bars produce better designs than their connected bars counterparts since they have more design freedom. Also, the smallest “solid” size variable $\alpha iq$ is 0.988249 for all floating bars runs and 0.983790 for all connected bars runs; and the largest “void” size variable $\alpha iq$ is 0.036276 for all floating bars runs and 0.038084 for all connected bar runs. This is an indication that the discreteness constraints are very effective in penalizing intermediate values of the size variables.

### Michell Cantilever.

We now present an example that shows that the proposed formulation can be readily extended to any number of materials. It corresponds to a cantilever frame considered in Michell's work [40], with design envelope and boundary conditions as shown in Fig. 8. The design envelope is meshed with 9700 elements. We enforce symmetry of the design with respect to the horizontal center line shown in Fig. 8. For this example, we use a tighter move limit of *m* = 0.2 to improve convergence. The initial design is made of 12 near-zero length bars of width *w* = 0.25 and with $\alpha iq=0.5,\u2009q=1,\u2026,12$.

We perform the optimization using four materials with Young's moduli $E1=6.5,\u2009E2=5.0,\u2009E3=4.5$, and $E4=3.5$ and physical densities $\gamma 1=0.55,\u2009\gamma 2=0.4,\u2009\gamma 1=0.35$, and $\gamma 2=0.25$. We use a weight-fraction limit of $wf*=0.028$ that allows us to obtain a design with four materials, as otherwise we would get designs made only of the stiffer material(s) if $wf*$ is larger, or made only of the weaker material(s) if $wf*$ is smaller. Figure 9 shows the optimal designs using one, two, three, and four materials. As expected, the designs improve as we increase the number of materials available. Also, the smallest solid size variable $\alpha iq$ is 0.999817 for all runs, and the largest void size variable $\alpha iq$ is 0.003536 for all runs, once again indicating the effectiveness of the penalization scheme.

### Three-Dimensional Cantilever Beam.

In the last example, we perform the optimization for a 3D cantilever beam. For this example, we migrated our matlab code to C++ using the deal.II library [41,42]. The design envelope, boundary conditions, and initial design are shown in Fig. 10. The design envelope is meshed with a regular grid of $64\xd732\xd732$ elements. We use the same two materials of the first example. The initial design consists of 42 bar of width *w* = 0.2 and with $\alpha 1q=\alpha 2q=0.5,\u2009q=1,\u2026,42$. We consider an initial design with connected bars as in the previous MBB beam example of Sec. 7.2. The optimal designs corresponding to weight-fraction limits of $wf*=$ 0.01, 0.02, 0.03, and 0.04 are shown in Table 2.

## Conclusions

We presented a method for the design via topology optimization of structures constructed as the union of geometric components, where each component is made of one of several available materials or removed from the design. Several examples that minimize the structural compliance subject to a constraint on the weight fraction demonstrate the proposed method. The available materials have different moduli but also different physical densities, hence a combination of materials is most advantageous for some weight fraction limits.

The examples demonstrate our method's effectiveness in producing structurally efficient multi-material designs. By penalizing intermediate size variables and enforcing the mutual material exclusion requirement as constraints in the optimization and not through the interpolation scheme, our technique makes it easier to incorporate any number of materials. Unlike density- and level set-based topology optimization methods for design with multiple materials, which produce material phases with geometries that are difficult to manufacture and assemble, the use of geometric components that are readily fabricated makes it easier to physically realize the multimaterial structures. Also, instead of imposing arbitrary volume fraction constraints on each of the available materials, we directly constrain the weight, which is a more natural design requirement.

We have demonstrated our method for two- and three-dimensional bars modeled with offset surfaces. Since the computation of the effective density *ρ*_{eff} of Eq. (21) is uncoupled from the calculation of the projected density *ρ _{q}* of Eqs. (3) and (4), the scheme to interpolate the material properties from various discrete geometric components should work for any geometry representation, and thus, we believe that our method can be applied to other component geometries; this will be demonstrated in future work. Future work will also focus on incorporating other important structural criteria.

## Acknowledgment

Support from the National Science Foundation, award CMMI-1634563 to conduct this work is gratefully acknowledged. We also thank Professor Krister Svanberg for kindly providing his MMA Matlab optimizer to perform the optimization.

## Funding Data

National Science Foundation (Grant No. CMMI-1634563).