Log In Sign Up

Printable Aggregate Elements

by   Jeremie Dumas, et al.

Aggregating base elements into rigid objects such as furniture or sculptures is a great way for designers to convey a specific look and feel. Unfortunately, there is no existing solution to help model structurally sound aggregates. The challenges stem from the fact that the final shape and its structural properties emerge from the arrangements of the elements, whose sizes are large so that they remain easily identifiable. Therefore there is a very tight coupling between the object shape, structural properties, and the precise layout of the elements. We present the first method to create aggregates of elements that are structurally sound and can be manufactured on 3D printers. Rather than having to assemble an aggregate shape by painstakingly positioning elements one by one, users of our method only have to describe the structural purpose of the desired object. This is done by specifying a set of external forces and attachment points. The algorithm then automatically optimizes a layout of user-provided elements that answers the specified scenario. The elements can have arbitrary shapes: convex, concave, elongated, and can be allowed to deform. Our approach creates connections between elements through small overlaps preserving their appearance, while optimizing for the global rigidity of the resulting aggregate. We formulate a topology optimization problem whose design variables are the positions and orientations of individual elements. Global rigidity is maximized through a dedicated gradient descent scheme. Due to the challenging setting -- number of elements, arbitrary shapes, orientation, and constraints in 3D -- we propose several novel steps to achieve convergence.


page 1

page 7

page 8

page 9

page 10

page 11


The Topology of Shapes Made with Points

In architecture, city planning, visual arts, and other design areas, sha...

Explicit topology optimization through moving node approach: beam elements recognition

Structural optimization (topology, shapes, sizing) is an important tool ...

Structural Design Using Laplacian Shells

We introduce a method to design lightweight shell objects that are struc...

Separating Variables in Bivariate Polynomial Ideals

We present an algorithm which for any given ideal I⊆𝕂 [x,y] finds all el...

Volumetric Procedural Models for Shape Representation

This article describes a volumetric approach for procedural shape modeli...

Organic Narrative Charts

Storyline visualizations display the interactions of groups and entities...

Learning Shapes by Convex Composition

We present a mathematical and algorithmic scheme for learning the princi...

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 111, ropes 222, or metal ornaments 333 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 user-specified 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 gradient-descent 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.


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], data-driven 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 level-of-detail 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 user-specified 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 r-functions [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

(1b) subject to

where the constraint forces elements to stay confined to a user-defined 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 gradient-based method. To do so, we have to define the compliance and its gradient in terms of the element parameters .

Figure 6. Our method pipeline. The element/sample parameters decide the element distribution , from which we compute the material densities for compliance analysis .

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 density-based 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.

We next detail our problem formulation (Section 4), and the numerical scheme we design to follow the compliance gradient and maximize the system’s rigidity (Section 5).

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 world-space 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:


, 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 non-rigid 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 non-linear 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.

Figure 7. Deformable element parameterization. The position of a sample dof is expressed hierarchically as in Equation 3. The root sample position is fixed in the local coordinate system of the element.

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 :


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 user-given 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:


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.


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:


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:


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.


The total material density at any point in is defined as the max of the densities induced by all the sample points :


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 non-differentiable 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.


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:


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:


where is the -th evaluation point of the quadrature rule for the cell , and is the associated weight. In our experiments, we used a 2-point 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:


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]:


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


Chain Rule

The current pipeline for computing the compliance from the element parameters can be summarized as follows:


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 row-vector, this can be expressed as:


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 th-sample th-coordinate is




When the in Equation 16 is reached by a single sample , then partial derivative exists, and is non-zero 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:


If sample belongs to element , then we can write:


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 word-space 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 non-linear and non-convex, with constraints that are also non-linear and non-convex. 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 gradient-based optimizer, which can handle such non-linear constraints and does not require a line-search 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.

[linenos, frame=none, framesep=0mm, numbersep=1mm, mathescape]lua function optimize() local num_iter = 30; local alpha = 3.0; local beta = 1.0; local alpha_min = 0.9; local beta_max = 2.0; for k = 1, num_iter do – Main MMA update Update(); end while alpha > alpha_min do alpha = math.max(alpha_min, alpha * 0.9); SetContinuationParameters(alpha, beta); ReparameterizeRotations(); for k = 1, num_iter do – Main MMA update Update(); end if alpha < 2.0 then – Sub-solver step: 10 BGFS updates BuildConnectivityGraph(); – Described in Algorithm 2 OptimizeConnectivity(10); end end while beta < beta_max do beta = math.max(beta_max, beta * 2.0); SetContinuationParameters(alpha, beta); ReparameterizeRotations(); for k = 1, num_iter do – Main MMA update Update(); end if alpha < 2.0 then – Sub-solver step: 10 BGFS updates BuildConnectivityGraph(); OptimizeConnectivity(10); end end end
Algorithm 1 Pseudo-code of our solver.

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.

(a) no continuation nor connectivity
(b) with continuation but no connectivity
(c) with continuation and connectivity
Figure 11. The effects of continuation and connectivity improvement steps. Our continuation method helps avoid stranded elements, while connectivity improvement reduces gaps between elements. LABEL:sub@fig:continuation:none With neither continuation nor connectivity some elements tend to isolate from the rest. LABEL:sub@fig:continuation:no_connectivity With continuation only the results improve but some elements tend to stretch out in different directions. LABEL:sub@fig:continuation:both Enabling both continuation and connectivity resolves those problems.

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 pseudo-code Algorithm 1).


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 multi-resolution optimization, where a coarse solution is produced and then refined. We start the gradient-based 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.

(a) element samples
(b) density grid,
(c) density grid,
Figure 15. Effect of the continuation parameter on the material densities. In the first iterations of the optimization, each element is set to occupy more physical space than it actually covers, to help the optimizer attract elements towards weak regions. Then, as we progressively reduce (by a factor of every iterations in our algorithm), the regular grid densities will match the actual physical object, up to the discretization error.


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 spring-mass 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 min-covering tree algorithm. Then, we interpret as a spring-mass network with rest length , and minimize its potential energy via steps of a BGFS solver.

Input: Element parameters , sample positions and their radii .
Output: Connectivity graph .
 // candidate connections
 // current distance between samples
 // union-find data structure for element ids
 // final set of edges for
1 foreach  do
2       ;
3       if  then
4             if  or or  then
5                   ;
                   // and are now neighbors in
// finalize: set target distances between samples
Algorithm 2 BuildConnectivityGraph

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 ill-conditioned 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 .

(a) 100 elements
(b) 150 elements
(c) 200 elements
Figure 19. Comparison of results with different number of elements under the same loading scenario (chair case). Ten different elements (shown in Figure 34) are used to produce results with respectively , , and elements.
(a) iteration 0
(b) iteration 60
(c) iteration 500
Figure 23. The effect of iterations. Initially the elements tend to be floating and detached from one another. With more iterations they gradually form a stronger structure.
(a) bananas
(b) lemons
(c) mangos
(d) oranges
(e) pears
(f) swords
(g) woodstick
(h) helix
(i) table problem domain
(j) bookend problem domain
Figure 34. Input elements and structural problem definitions.
(a) woodstick chair
(b) woodstick table
(c) pebble table
Figure 38. 3D printed results. All results are manufactured on a powder printer, using colors for LABEL:sub@fig:result:pebble_bridge and without for LABEL:sub@fig:result:woodstick_chair and LABEL:sub@fig:result:woodstick_bridge.
(a) small deformation,
(b) larger deformation,
Figure 41. A chair made of flexible noodles. Different results can be achieved by allowing different amounts of deformation.
Figure 42. A 3D printed iron throne made of swords. Cushion advised.
Figure 43. Bookend made of fruits. The loading scenario is shown in Figure (j)j.
Figure 44. A 3D-printed table made of pebbles supporting a teapot.
Figure 45. Different views of 3D printed chairs, the top one made of fruits printed in color, the bottom one made of deformable noodles, printed in white.

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.


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™ i7-5930K @ 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 single-thread 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
Fruits 510
Pebbles 420
Swords 420
Wood 420
Helix 420
Fruits 510
Pebbles 420
Swords 420
Wood 420
Table 1. Timings. Number mesh vertices in the output, number of elements, number of samples per element, time per iteration (max, min and mean), total time, number of iterations and grid size. The number of iterations depends on the continuation parameter (we used for the fruits, for the rest).

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 cross-section 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 off-line synthesizer with on-line 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.


  • [1]
  • Bendsøe and Sigmund [2004] M. P. Bendsøe and O. Sigmund. 2004. Topology Optimization: Theory, Methods and Applications. Springer Science + Business Media.
  • 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.
  • 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.
  • 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.
  • Christiansen et al. [2015] A. N. Christiansen, J. A. Bærentzen, M. Nobel-Jø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.
  • 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.
  • Dumas et al. [2015] J. Dumas, A. Lu, S. Lefebvre, J. Wu, and C. Dick. 2015. By-Example Synthesis of Structurally Sound Patterns. ACM Trans. Graph. 34, 4, Article 137 (07 2015), 12 pages.
  • Gal et al. [2007] R. Gal, O. Sorkine, T. Popa, A. Sheffer, and D. Cohen-Or. 2007. 3D Collage: Expressive Non-realistic Modeling. In Proceedings of the 5th International Symposium on Non-photorealistic Animation and Rendering (NPAR ’07). ACM, New York, NY, USA, 7–14.
  • 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.
  • Grassia [1998] F. S. Grassia. 1998. Practical Parameterization of Rotations Using the Exponential Map. Journal of Graphics Tools 3, 3 (01 1998), 29–48.
  • Guest [2015] J. K. Guest. 2015. Optimizing the Layout of Discrete Objects in Structures and Materials: A Projection-Based Topology Optimization Approach. Computer Methods in Applied Mechanics and Engineering 283 (01 2015), 330–351.
  • 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.
  • 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.
  • 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).
  • Jacobson et al. [2016] A. Jacobson, D. Panozzo, and others. 2016. libigl: A simple C++ geometry processing library. (2016).
  • 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.
  • 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.
  • 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.
  • Landes et al. [2013] P.-E. Landes, B. Galerne, and T. Hurtut. 2013. A Shape-Aware Model for Discrete Texture Synthesis. Computer Graphics Forum 32, 4 (07 2013), 67–76.
  • Lazarov et al. [2016] B. S. Lazarov, F. Wang, and O. Sigmund. 2016. Length Scale and Manufacturability in Density-Based Topology Optimization. Archive of Applied Mechanics 86, 1-2 (01 2016), 189–218.
  • Liu et al. [2014] T. Liu, S. Wang, B. Li, and L. Gao. 2014. A Level-Set-Based Topology and Shape Optimization Method for Continuum Structure Under Geometric Constraints. Structural and Multidisciplinary Optimization 50, 2 (03 2014), 253–273.
  • 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.
  • 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.
  • 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).
  • 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.
  • Norato et al. [2015] J. A. Norato, B. Bell, and D. Tortorelli. 2015. A Geometry Projection Method for Continuum-Based Topology Optimization With Discrete Elements. Computer Methods in Applied Mechanics and Engineering 293 (08 2015), 306–327.
  • 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.
  • Ritchie et al. [2015] D. Ritchie, S. Lin, N. D. Goodman, and P. Hanrahan. 2015. Generating Design Suggestions Under Tight Constraints With Gradient-Based Probabilistic Programming. Computer Graphics Forum 34, 2 (05 2015), 515–526.
  • 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.
  • Rvachev [1982] V. L. Rvachev. 1982. Theory of R-Functions and Some Applications. Naukova Dumka 552 (1982). (In Russian).
  • Sakurai and Miyata [2014] K. Sakurai and K. Miyata. 2014. Modelling of Non-Periodic Aggregates Having a Pile Structure. Computer Graphics Forum 33, 1 (02 2014), 190–198.
  • Schumacher et al. [2016] C. Schumacher, B. Thomaszewski, and M. Gross. 2016. Stenciling: Designing Structurally-Sound Surfaces With Decorative Patterns. Computer Graphics Forum 35, 5 (08 2016), 101–110.
  • Sorkine and Alexa [2007] O. Sorkine and M. Alexa. 2007. As-rigid-as-possible Surface Modeling. In Proceedings of the Fifth Eurographics Symposium on Geometry Processing (SGP ’07). Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, 109–116.
  • 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.
  • 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.
  • 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.
  • Wu et al. [2016] J. Wu, C. Dick, and R. Westermann. 2016. A System for High-Resolution Topology Optimization. IEEE Transactions on Visualization and Computer Graphics 22, 3 (03 2016), 1195–1208.
  • 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.
  • 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. Architecture-Scale Human-Assisted Additive Manufacturing. ACM Transactions on Graphics 34, 4 (07 2015), 88:1–88:8.
  • Zehnder et al. [2016] J. Zehnder, S. Coros, and B. Thomaszewski. 2016. Designing Structurally-Sound Ornamental Curve Networks. ACM Transactions on Graphics 35, 4 (07 2016), 1–10.
  • Zhang et al. [2012] J. Zhang, W. Zhang, J. Zhu, and L. Xia. 2012. Integrated Layout Design of Multi-Component Systems Using XFEM and Analytical Sensitivity Analysis. Computer Methods in Applied Mechanics and Engineering 245-246 (10 2012), 75–89.
  • 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.
  • 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.
  • 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.