1. Introduction
Depicting shapes by aggregating elements is a way for artists and designer to convey a specific look and feel, for instance creating furniture from wooden sticks ^{1}^{1}1http://www.marvelbuilding.com/uniquewoodenchairstackedstickscrossedstickchair.html, ropes ^{2}^{2}2https://www.dezeen.com/2014/12/22/azdventcalendarvermelhachaircampanabrothers/, or metal ornaments ^{3}^{3}3https://www.dezeen.com/2013/05/08/brazilianbaroqueexhibitionbythecampanabrothers/. This is often used to create a striking contrast between the appearance of the base elements and the purpose of the object formed by their aggregation. Famous examples are the candy house depicted in Hansel and Gretel, the iron throne (made of swords) of the Game of Thrones series, the paintings by Giuseppe Arcimboldo (portraits made of e.g. fruits, books, flowers), and the “accumulations” sculptures by artist Arman.
Aggregate elements are also ubiquitous in natural settings, e.g. a pile of rocks, a stash of fruits, a stone bridge.
Our goal is to allow the design of aggregate objects that can be manufactured and used in the real world. Challenges occur at two different scales: locally the elements have to be in contact with each others, but without significant overlaps that would destroy their appearance. Globally the elements have to form a rigid network so that the object is structurally sound.
Despite significant research on aggregates (see Section 2), no solution currently exists for our purpose. Available methods either optimize the manual assembly of specific 3D elements [Yoshida et al., 2015; Luo et al., 2015], or apply only to 2D domains such as planes or surfaces [Dumas et al., 2015; Martínez et al., 2015; Chen et al., 2016; Zehnder et al., 2016; Schumacher et al., 2016], or pack elements to closely approximate a target shape [Gal et al., 2007].
We present the first method to automatically generate 3D aggregates of elements arranged into globally rigid objects. We draw inspiration from modeling techniques that use rigidity as a design tool, by relying on topology optimization [Christiansen et al., 2015; Wu et al., 2016]. This lets the user specify the target shape through its desired structural properties. Our system thus requires as input only a set of 3D elements and a set of constraints: the output domain, a loading scenario (forces and attachments), and optionally spatial constraints (volumes to be avoided). The output of our system is a physically sound object that can be 3D printed. Our method is versatile and supports elements with shapes which are convex, concave, rigid or deformable.
Our algorithm optimizes aggregates by formulating a topology optimization problem with the position and orientation of individual elements as design variables. Each element contributes density within a material grid used for structural analysis. We maximize global rigidity through a dedicated gradient descent scheme that relates elements degrees of freedom to local material densities.
Our setting is made challenging by the 3D nature of the problem and the number of elements. The resolution of the simulation grid cannot be too high for reasonable speed and memory consumption. To achieve proper convergence while maintaining performance we rely on a continuation method that gradually reduces the individual influence of element samples over the density grid, and a connection step to eliminate small gaps that might be filtered away by the grid resolution.
The main contributions of this paper include:

The first system that automatically synthesizes structurally sound 3D aggregates from userspecified elements, structural objectives, and output domain constraints;

The formulation of a topology optimization problem to optimize the positions and orientations of elements;

The parameterization of element degrees of freedom for the support of convex, concave, rigid, as well as flexible elements of arbitrary shapes;

A gradientdescent based solver that links output material densities to input element degrees of freedom through a differentiable density field function of elements’ shapes;

A continuation method and a connection step to ensure convergence and quality while using coarse simulation grids for performance.
We verify our method through both numerical simulation and physical manufacturing of a variety of objects with different shapes and elements.
2. Related Work
Aggregate geometry has been extensively explored in a variety of fields in graphics, engineering, and manufacturing. We focus below on the works most relevant to ours.
Graphics
Aggregate geometry is ubiquitous and yet complex, and thus has been a research focus for modeling, animation, and rendering in computer graphics. Such repetitive geometry is often too complex for manual authoring and is better suited for procedural generation [Peytavie et al., 2009; Guérin et al., 2016], datadriven synthesis from exemplar elements [Sakurai and Miyata, 2014] and their distributions [Gal et al., 2007; Ma et al., 2011; Landes et al., 2013; Roveri et al., 2015]. While some of these techniques fill an existing shape by packing elements inside, they cannot fully synthesize a rigid shape having only desired structural properties as input.
Physical simulation often have to deal with the dynamics of aggregated shapes, which are computationally demanding and hard to control. Nailing specific aspects such as contact, friction, and gravity among a subset of elements can greatly help stability, control, and speed [Kaufman et al., 2008; Hsu and Keyser, 2010; Ritchie et al., 2015]. To reduce aliasing and increase rendering speed, [Cook et al., 2007] introduced a levelofdetail method to render aggregate geometry.
Design for manufacturing
Large objects can be manually assembled from specific base elements, e.g. sticks [Yoshida et al., 2015] or LEGO bricks [Luo et al., 2015]. These methods focus on computational support for manufacturing large, stable assemblies; they do not optimize the target shape itself and do not allow the user to choose the base elements.
Other approaches take into account both appearance and structure. Zehnder et al. [2016] focus on interactive editing, allowing the user to assemble curve elements into a fabricable connected network along a surface. The system automatically optimizes for deformations while indicating potential weak regions in the pattern. Several methods formulate a structure optimization problem combing userspecified loading scenarios and exemplar patterns [Dumas et al., 2015; Martínez et al., 2015; Chen et al., 2016]. These methods have been applied to 2D domains only, either planes or surfaces, but not to 3D volumes due to fundamental limitations in algorithms (e.g. [Dumas et al., 2015; Chen et al., 2016] are inherently 2D) or practical numerical feasibility (e.g. [Martínez et al., 2015] is computationally challenging for volumes).
The approach of Schumacher et al. [2016] covers a surface with holes having user specified shapes. The layout is optimized for both structural properties and aesthetics, and is coupled with the structural simulation through the rigidity of the surface triangles. Our method focuses on the complementary problem of optimizing a layout of elements that are in contact and form a globally rigid structure, in 3D, and using concave as well as deformable elements.
Note that most of the aforementioned techniques attempt to cover a domain with a pattern, while our technique, like [Martínez et al., 2015], attempts to completely synthesize a novel shape from structural and appearance constraints – in our case the base elements forming the structure itself.
Mechanical engineering
The problem of distributing solid elements or holes has been investigated in the field of mechanical engineering, as surveyed in [Lazarov et al., 2016].
There is a strong interest for jointly optimizing the position of embedded components and their support frame [Zhang et al., 2011, 2012; Kang and Wang, 2013; Xia et al., 2013; Wang et al., 2014; Zhang et al., 2015]. However, approaches that produce results comprised of aggregated elements only, i.e. without embedding them into an optimized support frame, have recently gained interest. These approaches relate the elements to an underlying simulation in a material grid, considering that each element generates a local density field [Overvelde, 2012; Guo et al., 2014; Norato et al., 2015; Zhang et al., 2016; Deng and Chen, 2016]. Avoiding overlaps between elements can be challenging. Explicit constraints can be used [Zhang et al., 2015; Gao et al., 2015; Kang et al., 2016], or elements can be combined by differentiable CSG [Chen et al., 2007, 2008; Liu et al., 2014] using rfunctions [Rvachev, 1982] or max operators relaxed as norms [Norato et al., 2015]. When applicable, the density field around each element can be modified to discourage overlaps [Overvelde, 2012; Guest, 2015].
Our approach is inspired by 2D methods using density fields around elements [Overvelde, 2012; Norato et al., 2015]. However, contrary to these techniques, we consider a 3D formulation and use elements that can have arbitrary shapes (elongated, concave, flexible), while these approaches typically focus on simpler, parametric shapes. We also target larger number of elements: the aforementioned techniques optimize for a few elements, while our smallest results have a hundred of them.
This makes convergence more challenging, and increases computational costs significantly. Besides the 3D formulation and parameterization of the element’s degrees of freedom, we introduce a continuation scheme on the element densities to allow for an improved convergence within reasonable computation times, as well as a connectivity improvement step.
3. Overview
Our algorithm takes as input a set of elements as colored 3D meshes, an output domain (size and shape), and a description of the desired object structural’s role: external loads to be supported as well as fixed attachment points. This is illustrated in Figures (b)b and (a)a
. Based on the exemplar elements and output domain, an initial configuration is created by distributing an estimated number of random instances from the exemplar elements.
Our optimizer iterates from this initial configuration to produce a rigid shape composed of the base elements, see Figure (c)c. Since the aggregated elements form a single connected structure optimized for rigidity, the object can be 3D printed and used as intended, see Figure (d)d.
The elements in the domain are defined by a set of parameters , which are the design variables used in the optimization. The structural objective we minimize is the compliance of the system, which can be understood as the inverse of the rigidy. It is computed from the current shape configuration and the user specified external forces . We thus seek to solve the following constrained optimization problem with respect to the set of parameters : equationparentequation
(1a)  
(1b)  subject to  
(1c) 
where the constraint forces elements to stay confined to a userdefined output domain, while and are box constraints on the element parameters , which are used to limit deformations on flexible elements and retain element’s centroids within the domain.
Due to the large number of design variables we seek to use a gradientbased method. To do so, we have to define the compliance and its gradient in terms of the element parameters .
The difficulty is that computing the compliance requires an underlying finite element simulation. Using the elements directly would impose an expensive union and tetrahedral remeshing step at every iteration. Instead, we build upon densitybased topology optimization methods [Bendsøe and Sigmund, 2004]. These approaches represent the optimized shape as a 3D density grid , where a density of is solid material and empty. A soft elastic material is assigned to the grid cells with density , and a rigid solid material is assigned to the grid cells with density
. In between, the material stiffness is interpolated according to
. Rather than directly manipulating the grid cell densities as design variables, we define the cell densities from the elements and their parameters . Intuitively, each element fills with solid material the density grid cells it overlaps, as illustrated in Figure 6. Given a differentiable density field function for the elements’ shapes, it is then possible to compute the compliance gradient through the chain rule.
However, using arbitrary shapes that optionally can deform makes finding a differentiable density field function challenging. We propose to approximate the geometry of the elements with multiple samples, each supporting a density kernel based on a differentiable smoothed Heaviside function. When the elements are allowed to deform, the samples are allowed a limited range of relative motions.
4. Problem formulation
In the following, we first discuss our choice of parameterization for elements in 3D, for rigid and flexible elements (Section 4.1). We next detail how we express domain constraints (Section 4.2) and how the densities and their gradients are computed from the point samples (Section 4.3). Finally we describe the formula for computing the compliance and its derivative (Section 4.4). For convenience, we give a table with all notations in the supplemental material.
4.1. Element Parameterization
We represent the elements in the domain by point samples [Ma et al., 2011]. The samples of an element can be selected manually, or they can be computed automatically. We chose the latter option, and used a centroidal Voronoi tessellation [Liu et al., 2009] to sample points regularly inside the input meshes (see Figure 6, left). Each sample is associated with a radius , that can vary from sample to sample, e.g. depending on how well the ball of center and radius locally approximates the shape of the input element.
We denote the set of all samples as , and their worldspace position by the matrix .
Rigid elements
When the user chooses to preserve the element’s shape during the optimization, they undergo only a rigid transformation parameterized by a translation and a rotation. For algorithm simplicity and computation efficiency, we do not allow elements to rescale, which can be achieved by providing input elements with different sizes.
Consider an element in the output domain. The element is defined by a set of parameters , where denotes the translational degrees of freedom, the rotational dofs, and the individual positions of the samples before transformation. The positions of the element samples in the output domain are given by the following relation:
(2) 
, where is the rotation matrix which is a function of the rotation parameters (exponential map in 3D, details in Section 5.3),
is a fixed linear transformation (fixed rotation + scaling),
is the translation matrix of the element centroid, and depends only on . Finally, are the (untransformed) sample positions. For rigid elements, is constant and determined once for each element type, while for elements allowed to deform is a matrix representing all the positional dofs for element .Deformable elements
Optionally, the user may choose to allow elements to deform during the optimization process. In addition to the rigid degree of freedom associated to an element , we use the sampled points inside each element (coordinates in the local element reference frame) as control points to drive a nonrigid deformation. At the end of the process we use ARAP [Sorkine and Alexa, 2007] and the implementation from libigl [Jacobson et al., 2016] to recover the final shape.
To prevent excessive deformation during the optimization, we limit how much the control points can stray from the initial configuration. One approach would be to add a set of constraints on the positions of the constituent samples in Equation 1. However, we have found that adding such nonlinear constraints makes the optimization problem too challenging to solve: the solver either tends to get stuck (constraints too tight), or produces excessive deformations (constraints too loose), without striking a good balance.
Instead we propose to rely on a parameterization of the point sample positions that eliminates the need for complex constraints. The parameterization we propose is described in Figure 7. We observe that most often, the geometry of deformable elements can be well captured by a skeleton. We thus embed an articulated skeleton within the element’s shape. The lengths of the skeleton’s bones remain constant, preventing longitudinal distortions. The angles can be restricted through simple box constraints (see in Equation 1) to prevent excessive rotations at joints.
Let us assume that we are given a skeleton (tree) connecting the samples within an element . (We describe an algorithm for the automatic construction of in the supplemental material.) We define the position of each sample relative to its parent, parameterized by a rotation :
(3) 
Note that is defined based on the input, undeformed element, but is fixed through the optimization. The root of the tree, , is also fixed/constant through the optimization (without loss of generality since each element has already a translational degree of freedom). In the undeformed configuration, for all samples and the reconstructed surface is the same as input element. Figure 41 shows a result where different limits are set for the acceptable deformations.
One limitation of the skeleton approach, however, is that we are restricted to a set of skeleton–like elements and cannot directly support deformable elements such as sheets or flexible rings.
4.2. Domain Constraints
Similarly to [Ma et al., 2011], we seek to restrict the sample positions to the interior of a usergiven boundary. In addition, each sample must be at a distance at least from the boundary surface, otherwise part of the element will fall outside the domain. Let be the signed distance from a point to the domain boundary, with negative values inside. We define the following functional:
(4) 
It follows that we have if and only if all the samples lie inside the domain at a distance at least .
We always activate the boundary constraint for the domain. Optionally, it can be used to forbid regions of the output domain, e.g. the surfaces along chair back and seat in Figure 5, for seating comfort.
4.3. Material Densities
To compute the compliance of the shape being synthesized, it is necessary to assign a density everywhere in space from the current configuration of the point samples.
Representation
Let be a point sample in the output domain, and let be the material density associated to the sample . We define
as a radial basis function, which depends only on the distance from the sample position
, and is parameterized by the sample radius . The RBF is chosen to have a compact support so that each element only has a local influence. In practice, we use a smoothed Heaviside step function:(5) 
where controls the smoothness of the approximation, and controls the radius of influence. In our implementation we consider that when . The term considers whether is inside the element so that small gaps will not be erroneously filled during FEM in Section 4.4:
(6) 
In contrast to works in mechanical engineering that use elements for structural optimization (see Section 2), in Equation 5 multiple samples correspond to a same element, and are driven by the same element’s parameters (a rigid transformation, or a skeletal parameterization if deformable). In addition, we propose to use and in a continuation scheme during optimization to improve convergence, as will be detailed in Section 5.2.
Gathering
The total material density at any point in is defined as the max of the densities induced by all the sample points :
(7) 
Defining the material density using a maximum instead of a sum discourages elements to overlap. Indeed, overlapping two elements under a max decreases the total density of the system, which increases the compliance. Thus, solutions where samples do not overlap have a lower compliance and are preferred. Conversely, minimizing the compliance tends to pull the elements closer together, as inter–element gaps result in fragile structures. This is precisely the behaviour we intend: pulling the elements together while discouraging large overlaps.
It should be noted that the function is technically not differentiable. A common workaround is to resort to a smoothed formulation, e.g. using a norm , with a high or . The drawback of this approach is that the actual density at any point in space depends on the number of samples , as the norm inevitably computes a form of weighted average over the domain. To retrieve a good approximation of the function, it is necessary to increase the exponent when there are a high number of samples, which leads to numerical inaccuracies.
In practice, we have found that simply ignoring this theoretical issue and retaining a hard formulation does not impede the overall gradient computation. Indeed, the only nondifferentiable points of the function are the points which are closest and equidistant from two different samples (assuming is the same for all samples). In other words, they are located on the edges of the Voronoi diagram formed by the current distribution of samples. Due to the limited numerical precision of floating point operations these singular points never occur in practice.
Discretization
In order to compute the compliance of the system via a finite element method, the material densities are discretized in a regular 3D grid, in our case composed of linear H8 cube elements. Let denote a cell in the regular grid , and let be its associated material density. The discretized cell density is given by the following relation:
(8) 
where is the domain covered by the grid cell .
Equation 8 simply means that is defined as the average density over the grid cell . In practice, the integral (8) can be computed either 1) analytically by an exact expression of the integral of , or 2) numerically by means of a Gaussian quadrature rule, i.e. evaluating at specified points inside the cell. While the expression of for a given sample is integrable analytically, the use of the function makes it difficult to derive a simple analytic expression of the resulting integral. For this reason, we opted for numerical integration of the expression given in Equation 8:
(9) 
where is the th evaluation point of the quadrature rule for the cell , and is the associated weight. In our experiments, we used a 2point quadrature rule (i.e. 8 points per cube).
4.4. Compliance and Sensitivities
Given the 3D density grid, and knowing the external forces applied to the system, we can now compute how the system deforms under load. The discrete displacement field is obtained with the finite element method, solving the equilibrium equation:
(10) 
where is the global stiffness matrix of the system, assembled from the individual matrices of every grid cell , and where is the stiffness matrix for the base solid material. We use a solver similar to Wu et al. [2016] to solve for Equation 10.
The compliance of the system is then computed as [Bendsøe and Sigmund, 2004]:
(11) 
where is the displacement vector associated to the node of the grid cell induced by the external force .
Using the adjoint method ([Bendsøe and Sigmund, 2004], §1.2.3), the partial derivatives of the compliance can be expressed as
(12) 
Chain Rule
The current pipeline for computing the compliance from the element parameters can be summarized as follows:
(13) 
where denotes the element parameters, denotes the sample positions in the output domain, are the grid cell densities, and is the scalar value of the compliance objective function. Equation 12 explains how to obtain the compliance gradient . Given , the gradients and can be computed efficiently via the chain rule. If we note the gradient as a rowvector, this can be expressed as:
(14) 
Note that the Jacobian matrices and are sparse matrices, so the products in Equation 14 can be implemented efficiently.
In the following, we explain how to compute the Jacobian matrices and . The sensitivity of the cell density with respect to the thsample thcoordinate is
(15) 
with
(16) 
When the in Equation 16 is reached by a single sample , then partial derivative exists, and is nonzero when . It is then equal to the partial derivative , whose expression is obtained by deriving Equation 5 (detailed in the supplemental material).
For the element positions, the partial derivative with respect to element parameters can be expressed as:
(17) 
If sample belongs to element , then we can write:
(18) 
The expression of Equation 18 can be simplified whether corresponds to the rotation parameter, the translation, or the sample positions. The different cases are presented in the supplemental.
Regarding deformable elements, in order to compute the gradient of the wordspace sample position with respect to the sample parameters (Equation 18), we need to be able to compute when is one of the that parameterize the sample position . The complete derivation of this partial derivative is also detailed in the supplemental material.
We now have all the definitions that form the basis of our formulation. We proceed with describing our numerical solver in Section 5.
5. Solver
The problem defined in Equation 1 is nonlinear and nonconvex, with constraints that are also nonlinear and nonconvex. In addition, computing the objective function involves an expensive FEM computation to compute the equilibrium state (Equation 10). A solver of choice for solving topology optimization problems is MMA [Svanberg, 1987], a gradientbased optimizer, which can handle such nonlinear constraints and does not require a linesearch at every update step (thus avoiding solving the expensive Equation 10). In practice, we use a custom implementation of MMA in C which allows us to use the continuation method mentioned in Section 5.2
The complete pseudo code is given in Algorithm 1.
5.1. Initialization
Number of elements/samples
For efficiency reasons we maintain the number of elements/samples constant during the main optimization loop: Adding or removing elements would reset the gradient information accumulated by the MMA optimizer over the previous iterations.
The user can thus specify the desired number of elements directly, or indirectly through a ratio of the output volume that should be occupied by elements (the overall volume does not change during optimization as the number of elements is kept constant and deformation constraints limit volume increase).
Placing elements/samples
We randomly place the initial elements and spread them out by computing a centroidal Voronoi tessellation (CVT) of their centers. As this ignores the actual shape of the elements, some may lie partially outside the domain. We improve the initial distribution by next minimizing the CVT objective on the samples composing the elements, under the boundary constraint (Equation 4
). This pulls all samples back inside the output domain boundary, while preserving a uniform distribution. An illustration of the resulting initial distribution is presented in
Figure 23.5.2. Iterations
After initialization the optimizer enters an iterative loop that computes the compliance gradient with respect to element parameters (Section 4.4) and performs a descent step, following the MMA method. We however embed two additional mechanisms. The first is a continuation scheme that affords for improved convergence in our setting (Section 5.2). The other is a connectivity step that explicitly encourages connections between neighboring elements (Section 5.2). Finally, we provide some details on how 3D rotations are handled during optimization (Section 5.3).
The total number of iterations is determined by the continuation scheme (see Section 5.2, and pseudocode Algorithm 1).
Continuation
A crucial question when filling the density grid from the elements is how large their region of influence should be. This choice impacts convergence significantly.
The parameters controlling the spatial influence of the elements are and (to a lesser extent) in the sample’s influence function given in Equation 5. directly controls the width of the RBF. If , then the material density given by tightly corresponds to the density of the physical element represented by the sample.
Setting seems to be an obvious choice. However, elements can only be attracted towards regions where they have a non zero contribution. Thus, using a tight RBF will impede the ability of the optimizer to pull elements towards regions of high compliance. In the worst case, some elements can end up floating in regions of low compliance while never be attracted where they would be most useful, as illustrated in Figure (a)a. By setting we can artificially enlarge the region an element spans. This unfortunately creates a situation where the RBFs are no longer representative of the physical elements.
To mitigate these effects while enabling a good convergence, we introduce a continuation method on the radius of influence multiplier , as illustrated in Figure 15. This is akin to multiresolution optimization, where a coarse solution is produced and then refined. We start the gradientbased optimization with a high value of (3 in our implementation), and progressively decrease it every 30 iterations, multiplying by a factor of until it reaches a minimum value. We set the minimum value to to encourage small overlaps between adjacent elements – which is necessary in order to print aggregate geometry. We treat – that controls the smoothness of the density change – in a similar fashion (see Algorithm 1).
The effect of the continuation parameter can be observed in Figure 11. Without it, some elements end up stranded away from the main structure.
Connectivity
Even though the compliance optimization will naturally discourage small gaps between elements, the gaps might still occur due to the discretization. Specifically, a gap sufficiently smaller than a grid cell can get filtered away and go unnoticed by the optimizer. To reduce this issue we embed in the optimization loop the following mechanism that explicitly encourages contacts.
The connectivity improvement algorithm works as follows: 1) detect pairs of samples from distinct elements that we would like to move towards each other, and define for each pair a target length so that elements overlap, 2) minimize the potential energy of a springmass system on the graph defined by those edges (we do 10 steps of BFGS every time at every continuation step, see Algorithm 1).
We compute a graph following Algorithm 2, inspired by Kruskal’s mincovering tree algorithm. Then, we interpret as a springmass network with rest length , and minimize its potential energy via steps of a BGFS solver.
5.3. Parameterization of 3D Rotations
We have to differentiate the sample positions with respect to the element parameters , including 3D rotations. We parameterize them using exponential maps, and in particular the formulation given by [Grassia, 1998] which computes the exponential map from to via an intermediate quaternion representation. The authors provide a C implementation for computing the partial derivatives of the rotation matrix with respect to the exponential map vector in . Since the exponential maps can become illconditioned if the rotation is too high, we perform a check every 30 iterations of our algorithm, and we reparameterize rotations that have become too large (angle ). More specifically, we set and .









6. Results
In this section we first discuss the influence of various parameters, provide timings and statistics on our results, and then show a number of 3D printed results (on a ZCorp 450 color powder printer). We also provide renderings of additional results. The accompanying video shows animations of the optimization process as well as rotating views of results. The inputs used throughout this section are visible in Figure 5 and Figure 34.
Analysis
The main parameter the user has to choose is the number of elements, and thus, indirectly, the solid volume percentage. Figure 19 shows the influence on the result for the chair case. As can be seen, the optimizer successfully arranges the elements in rigid connected structures, even on the result with the smallest number of elements. Of course, more elements affords for a more sturdy structure as more material is available.
Figure 23 visualizes the behavior of the optimization algorithm as it iterates. Animations are available in the accompanying video.
Figure 11 shows the benefits of our continuation and connectivity methods (Sections 5.2 and 5.2). In particular, a few stranded elements are visible when continuation is disabled, while the connectivity step encourages small overlaps between adjacent elements.
Timings are reported in Table 1. All results were obtained on an Intel^{®} Core™ i75930K @ 3.50GHz, 64 GB RAM. Note that the first iterations are slower because RBFs have a larger support due to our continuation scheme (our evaluation procedure uses a singlethread on the CPU). As decreases, the iterations become faster, and are dominated by the cost of solving Equation 10.
Synthesized shapes
As synthesis is fully automatic, it is very simple for the user to produce a variety of results by combining different elements and structural problems. We synthesized and 3D printed several such examples.
The chair is a common daily object with different geometry components, including thin legs and planar seats and backs. Our method automatically assembles a chair from its structural definition, out of different elements including mixture of fruits (Figure 45, top), wood sticks (Figure (a)a), pebbles (Figure 5), long flexible noodles (Figure 41, Figure 45 bottom) and even swords (Figure 42). This later case is challenging as it uses a large number of sharp, elongated swords. Our method succeeds in maintaining a sound structure where the elements remain easily identifiable.
We also synthesize and 3D print table structures from elongated sticks (Figure (b)b) and pebbles (Figure (c)c, Figure 44), and created a bookend made of fruits, rendered in Figure 43.
All these results were printed or rendered without any change after optimization. However, it is worth mentioning that since the output is made of elements, it would be simple to create a tool allowing the user to select, move, scale, delete or add some elements.
Exemplar  # Vertices  # Elements  # Samples  Max (s)  Min (s)  Mean (s)  Total ()  # Iter  Grid size 

Chair  
Fruits  510  
Pebbles  420  
Swords  420  
Wood  420  
Helix  420  
Table  
Fruits  510  
Pebbles  420  
Swords  420  
Wood  420 
7. Limitations and future work
In this work we propose a novel approach for the modeling of complex aggregates of elements, which uses rigidity as a design tool, which supports rigid and deformable elements with arbitrary shapes, and synthesizes results that can be 3D printed.
One limitation of our approach is that there is a difference between what the optimizer sees (the underlying density grid) and the actual final geometry (aggregate of elements). In particular, even
though our results typically have a single connected component there may be fragile connections between elements that barely touch each others. For this reason, it might be necessary to further reinforce the structure.
One possibility is to use the graph built during connectivity improvement (Section 5.2) to add struts between neighboring elements. An illustrative example is shown in the inset figure (struts in red) for reinforcing a noodle chair result. Another possibility would be to slightly scale and move elements to enforce a minimal crosssection for all contacts, or to add struts similarly to [Stava et al., 2012] – however a full FEM simulation would be expensive on our models.
Our method is currently an offline synthesizer with online preview (each iteration takes around 10 seconds). We attempted to achieve a good balance between precision (to capture intricate geometries) and speed. However the optimization is not yet interactive.
The size of the objects we can 3D print in one piece is limited, and thus we can only produce miniatures of e.g. the chair. It would be interesting to consider printing such shapes in several parts that can be later assembled, the contacts between elements being a natural location to embed connectors.
Finally, as future work we would like to explore more user controls, possibly including pausing the optimization, making a few changes, and restarting after these additional user edits. Other controls would include direction and scale fields, as well as encouraging symmetries, to further refine the aesthetics of the results.
References
 [1]
 Bendsøe and Sigmund [2004] M. P. Bendsøe and O. Sigmund. 2004. Topology Optimization: Theory, Methods and Applications. Springer Science + Business Media. https://doi.org/10.1007/9783662050866
 Chen et al. [2008] J. Chen, M. Freytag, and V. Shapiro. 2008. Shape Sensitivity of Constructively Represented Geometric Models. Computer Aided Geometric Design 25, 7 (10 2008), 470–488. https://doi.org/10.1016/j.cagd.2008.01.005
 Chen et al. [2007] J. Chen, V. Shapiro, K. Suresh, and I. Tsukanov. 2007. Shape Optimization With Topological Changes and Parametric Control. Internat. J. Numer. Methods Engrg. 71, 3 (2007), 313–346. https://doi.org/10.1002/nme.1943
 Chen et al. [2016] W. Chen, X. Zhang, S. Xin, Y. Xia, S. Lefebvre, and W. Wang. 2016. Synthesis of Filigrees for Digital Fabrication. ACM Transactions on Graphics 35, 4 (07 2016), 1–13. https://doi.org/10.1145/2897824.2925911
 Christiansen et al. [2015] A. N. Christiansen, J. A. Bærentzen, M. NobelJørgensen, N. Aage, and O. Sigmund. 2015. Combined Shape and Topology Optimization of 3D Structures. Computer and Graphics 46, C (2015), 25–35.
 Cook et al. [2007] R. L. Cook, J. Halstead, M. Planck, and D. Ryu. 2007. Stochastic Simplification of Aggregate Detail. ACM Transactions on Graphics 26, 3 (07 2007), 79. https://doi.org/10.1145/1276377.1276476
 Deng and Chen [2016] J. Deng and W. Chen. 2016. Design for Structural Flexibility Using Connected Morphable Components Based Topology Optimization. Science China Technological Sciences 59, 6 (03 2016), 839–851. https://doi.org/10.1007/s1143101660270
 Dumas et al. [2015] J. Dumas, A. Lu, S. Lefebvre, J. Wu, and C. Dick. 2015. ByExample Synthesis of Structurally Sound Patterns. ACM Trans. Graph. 34, 4, Article 137 (07 2015), 12 pages. https://doi.org/10.1145/2766984
 Gal et al. [2007] R. Gal, O. Sorkine, T. Popa, A. Sheffer, and D. CohenOr. 2007. 3D Collage: Expressive Nonrealistic Modeling. In Proceedings of the 5th International Symposium on Nonphotorealistic Animation and Rendering (NPAR ’07). ACM, New York, NY, USA, 7–14. https://doi.org/10.1145/1274871.1274873
 Gao et al. [2015] H.H. Gao, J.H. Zhu, W.H. Zhang, and Y. Zhou. 2015. An Improved Adaptive Constraint Aggregation for Integrated Layout and Topology Optimization. Computer Methods in Applied Mechanics and Engineering 289 (06 2015), 387–408. https://doi.org/10.1016/j.cma.2015.02.022
 Grassia [1998] F. S. Grassia. 1998. Practical Parameterization of Rotations Using the Exponential Map. Journal of Graphics Tools 3, 3 (01 1998), 29–48. https://doi.org/10.1080/10867651.1998.10487493
 Guest [2015] J. K. Guest. 2015. Optimizing the Layout of Discrete Objects in Structures and Materials: A ProjectionBased Topology Optimization Approach. Computer Methods in Applied Mechanics and Engineering 283 (01 2015), 330–351. https://doi.org/10.1016/j.cma.2014.09.006
 Guo et al. [2014] X. Guo, W. Zhang, and W. Zhong. 2014. Doing Topology Optimization Explicitly and Geometrically—A New Moving Morphable Components Based Framework. Journal of Applied Mechanics 81, 8 (05 2014), 081009. https://doi.org/10.1115/1.4027609
 Guérin et al. [2016] E. Guérin, E. Galin, F. Grosbellet, A. Peytavie, and J.D. Génevaux. 2016. Efficient modeling of entangled details for natural scenes. Computer Graphics Forum 35, 7 (2016), 257–267. https://doi.org/10.1111/cgf.13023
 Hsu and Keyser [2010] S.W. Hsu and J. Keyser. 2010. Piles of Objects. In ACM SIGGRAPH Asia 2010 papers on  SIGGRAPH ASIA '10. Association for Computing Machinery (ACM). https://doi.org/10.1145/1882262.1866181
 Jacobson et al. [2016] A. Jacobson, D. Panozzo, and others. 2016. libigl: A simple C++ geometry processing library. (2016). http://libigl.github.io/libigl/.
 Kang and Wang [2013] Z. Kang and Y. Wang. 2013. Integrated Topology Optimization With Embedded Movable Holes Based on Combined Description by Material Density and Level Sets. Computer Methods in Applied Mechanics and Engineering 255 (03 2013), 1–13. https://doi.org/10.1016/j.cma.2012.11.006
 Kang et al. [2016] Z. Kang, Y. Wang, and Y. Wang. 2016. Structural Topology Optimization With Minimum Distance Control of Multiphase Embedded Components by Level Set Method. Computer Methods in Applied Mechanics and Engineering 306 (07 2016), 299–318. https://doi.org/10.1016/j.cma.2016.04.001
 Kaufman et al. [2008] D. M. Kaufman, S. Sueda, D. L. James, and D. K. Pai. 2008. Staggered Projections for Frictional Contact in Multibody Systems. ACM Transactions on Graphics 27, 5 (12 2008), 1. https://doi.org/10.1145/1409060.1409117
 Landes et al. [2013] P.E. Landes, B. Galerne, and T. Hurtut. 2013. A ShapeAware Model for Discrete Texture Synthesis. Computer Graphics Forum 32, 4 (07 2013), 67–76. https://doi.org/10.1111/cgf.12152
 Lazarov et al. [2016] B. S. Lazarov, F. Wang, and O. Sigmund. 2016. Length Scale and Manufacturability in DensityBased Topology Optimization. Archive of Applied Mechanics 86, 12 (01 2016), 189–218. https://doi.org/10.1007/s0041901511064
 Liu et al. [2014] T. Liu, S. Wang, B. Li, and L. Gao. 2014. A LevelSetBased Topology and Shape Optimization Method for Continuum Structure Under Geometric Constraints. Structural and Multidisciplinary Optimization 50, 2 (03 2014), 253–273. https://doi.org/10.1007/s0015801410457
 Liu et al. [2009] Y. Liu, W. Wang, B. Lévy, F. Sun, D.M. Yan, L. Lu, and C. Yang. 2009. On Centroidal Voronoi Tessellation—energy Smoothness and Fast Computation. ACM Transactions on Graphics 28, 4 (08 2009), 1–17. https://doi.org/10.1145/1559755.1559758
 Luo et al. [2015] S.J. Luo, Y. Yue, C.K. Huang, Y.H. Chung, S. Imai, T. Nishita, and B.Y. Chen. 2015. Legolization: Optimizing LEGO Designs. ACM Transactions on Graphics 34, 6 (10 2015), 1–12. https://doi.org/10.1145/2816795.2818091
 Ma et al. [2011] C. Ma, L.Y. Wei, and X. Tong. 2011. Discrete Element Textures. In ACM SIGGRAPH 2011 papers on  SIGGRAPH '11. Association for Computing Machinery (ACM). https://doi.org/10.1145/1964921.1964957
 Martínez et al. [2015] J. Martínez, J. Dumas, S. Lefebvre, and L.Y. Wei. 2015. Structure and Appearance Optimization for Controllable Shape Design. ACM Transactions on Graphics 34, 6 (10 2015), 1–11. https://doi.org/10.1145/2816795.2818101
 Norato et al. [2015] J. A. Norato, B. Bell, and D. Tortorelli. 2015. A Geometry Projection Method for ContinuumBased Topology Optimization With Discrete Elements. Computer Methods in Applied Mechanics and Engineering 293 (08 2015), 306–327. https://doi.org/10.1016/j.cma.2015.05.005
 Overvelde [2012] J. T. Overvelde. 2012. The Moving Node Approach in Topology Optimization. Master’s thesis. TU Delft, Delft University of Technology.
 Peytavie et al. [2009] A. Peytavie, E. Galin, J. Grosjean, and S. Merillou. 2009. Procedural Generation of Rock Piles using Aperiodic Tiling. Computer Graphics Forum 28, 7 (2009), 1801–1809. https://doi.org/10.1111/j.14678659.2009.01557.x
 Ritchie et al. [2015] D. Ritchie, S. Lin, N. D. Goodman, and P. Hanrahan. 2015. Generating Design Suggestions Under Tight Constraints With GradientBased Probabilistic Programming. Computer Graphics Forum 34, 2 (05 2015), 515–526. https://doi.org/10.1111/cgf.12580
 Roveri et al. [2015] R. Roveri, A. C. Öztireli, S. Martin, B. Solenthaler, and M. Gross. 2015. Example Based Repetitive Structure Synthesis. Computer Graphics Forum 34, 5 (08 2015), 39–52. https://doi.org/10.1111/cgf.12695
 Rvachev [1982] V. L. Rvachev. 1982. Theory of RFunctions and Some Applications. Naukova Dumka 552 (1982). (In Russian).
 Sakurai and Miyata [2014] K. Sakurai and K. Miyata. 2014. Modelling of NonPeriodic Aggregates Having a Pile Structure. Computer Graphics Forum 33, 1 (02 2014), 190–198. https://doi.org/10.1111/cgf.12266
 Schumacher et al. [2016] C. Schumacher, B. Thomaszewski, and M. Gross. 2016. Stenciling: Designing StructurallySound Surfaces With Decorative Patterns. Computer Graphics Forum 35, 5 (08 2016), 101–110. https://doi.org/10.1111/cgf.12967
 Sorkine and Alexa [2007] O. Sorkine and M. Alexa. 2007. Asrigidaspossible Surface Modeling. In Proceedings of the Fifth Eurographics Symposium on Geometry Processing (SGP ’07). Eurographics Association, AirelaVille, Switzerland, Switzerland, 109–116. http://dl.acm.org/citation.cfm?id=1281991.1282006
 Stava et al. [2012] O. Stava, J. Vanek, B. Benes, N. Carr, and R. Měch. 2012. Stress Relief: Improving Structural Strength of 3D Printable Objects. ACM Transactions on Graphics 31, 4 (07 2012), 1–11. https://doi.org/10.1145/2185520.2185544
 Svanberg [1987] K. Svanberg. 1987. The Method of Moving Asymptotes—a New Method for Structural Optimization. Internat. J. Numer. Methods Engrg. 24, 2 (02 1987), 359–373. https://doi.org/10.1002/nme.1620240207
 Wang et al. [2014] Y. Wang, Z. Luo, X. Zhang, and Z. Kang. 2014. Topological Design of Compliant Smart Structures With Embedded Movable Actuators. Smart Materials and Structures 23, 4 (03 2014), 045024. https://doi.org/10.1088/09641726/23/4/045024
 Wu et al. [2016] J. Wu, C. Dick, and R. Westermann. 2016. A System for HighResolution Topology Optimization. IEEE Transactions on Visualization and Computer Graphics 22, 3 (03 2016), 1195–1208. https://doi.org/10.1109/tvcg.2015.2502588
 Xia et al. [2013] L. Xia, J. Zhu, W. Zhang, and P. Breitkopf. 2013. An Implicit Model for the Integrated Optimization of Component Layout and Structure Topology. Computer Methods in Applied Mechanics and Engineering 257 (04 2013), 87–102. https://doi.org/10.1016/j.cma.2013.01.008
 Yoshida et al. [2015] H. Yoshida, S. Igarashi, T. Igarashi, Y. Obuchi, Y. Takami, J. Sato, M. Araki, M. Miki, K. Nagata, and K. Sakai. 2015. ArchitectureScale HumanAssisted Additive Manufacturing. ACM Transactions on Graphics 34, 4 (07 2015), 88:1–88:8. https://doi.org/10.1145/2766951
 Zehnder et al. [2016] J. Zehnder, S. Coros, and B. Thomaszewski. 2016. Designing StructurallySound Ornamental Curve Networks. ACM Transactions on Graphics 35, 4 (07 2016), 1–10. https://doi.org/10.1145/2897824.2925888
 Zhang et al. [2012] J. Zhang, W. Zhang, J. Zhu, and L. Xia. 2012. Integrated Layout Design of MultiComponent Systems Using XFEM and Analytical Sensitivity Analysis. Computer Methods in Applied Mechanics and Engineering 245246 (10 2012), 75–89. https://doi.org/10.1016/j.cma.2012.06.022
 Zhang et al. [2011] W. Zhang, L. Xia, J. Zhu, and Q. Zhang. 2011. Some Recent Advances in the Integrated Layout Design of Multicomponent Systems. Journal of Mechanical Design 133, 10 (2011), 104503. https://doi.org/10.1115/1.4005083
 Zhang et al. [2016] W. Zhang, J. Yuan, J. Zhang, and X. Guo. 2016. A New Topology Optimization Approach Based on Moving Morphable Components (MMC) and the Ersatz Material Model. Structural and Multidisciplinary Optimization 53, 6 (06 2016), 1243–1260. https://doi.org/10.1007/s0015801513723
 Zhang et al. [2015] W. Zhang, W. Zhong, and X. Guo. 2015. Explicit Layout Control in Optimal Design of Structural Systems With Multiple Embedding Components. Computer Methods in Applied Mechanics and Engineering 290 (06 2015), 290–313. https://doi.org/10.1016/j.cma.2015.03.007