1. Introduction
In structural design, shells are considered to be one of the most efficient structures because they can be simultaneously lightweight and robust. Shells are common in additive fabrication applications because using shells instead of solids reduces the cost of material and decreases the fabrication time.
An optimallyshaped shell can carry its load relying only on tensile/compression forces, with no bending involved, which is very efficient in terms of the required material. These types of shells are commonly found in architecture (domes). However, the shape of the shell may be determined by considerations other than its loadcarrying properties. For example, the shape of the airplane is determined by its aerodynamics; the shape of the car body both by the aerodynamics and aesthetics; the top of a table or a shelf is flat, as objects need to be placed on it; the artistic intent primarily determines the shape of a lamp or a statue.
Shell structures with shapes fixed by considerations other than loading are often reinforced by additional means, most commonly increasing thickness in critical areas or adding ribs (Fig. 2). Formally, a common optimization goal for a shell reinforcement structure is to minimize the weight of the added material while keeping the maximal stress of the structure bounded. The first ensures the structure remains lightweight, while the second prevents structural failure.
This problem has been well studied for two dimensions in the limit of low volumes. In 2D, the optimal layout forms a pattern of orthogonal lines (HenckyPrandtl net) and, for a given layout, the minimum weight structure has all members fully stressed. These structures form, in the limit of low volumes, classical Michell structures and can be obtained by solving a convex optimization problem. It is also well understood for pure bending problems for plates, i.e., the special case of flat shells with loads orthogonal to the surface; in this case, it also reduces to a different convex problem.
The situation is far more complex for the reinforcement of fixed curved shells embedded in three dimensions that we consider. For these shells, the weightoptimal structure may be locally either beamlike, forming ribs aligned with stress directions, or membranelike, forming variablethickness walls with no perforation [sigmund2016non]. The first case typically corresponds to bendingdominated regions while the second to areas dominated by inplane forces. The optimal local structure is determined by the surface shape, the supports, and the loads.
In this general case, the problem is no longer convex and cannot be optimally solved either by methods that assume that result is only a variable thickness shell or by Michelltruss type methods.
In this paper, we propose a novel efficient computational method for constructing optimized reinforcement structures for shells naturally producing a full range of behaviors spanning the space between variablethickness shell and ribtype reinforcement.
We partition the problem into three steps: (1) determine the field of (approximately) optimal stress directions; (2) construct the skeleton of the reinforcement structure that follows these directions, forming polygonal (predominantly) quad cells aligned with the field; (3) optimize how material is distributed inside the cells.
Contributions.

For computing the field of optimal stress directions, we developed a generalization of HenckyPrandtl nets which takes bending into account and can still be solved by minimizing a convex energy.

For material distribution optimization, we use a lowparametric structure model for cells to efficiently optimize the distribution of the material. As the global optimization problem is fundamentally nonconvex (we discuss the reasons on Section 3), to solve it we propose an efficient global/local method which shows stable and fast convergence behavior.
We demonstrate our approach by optimizing several 3D shapes. This evaluation shows our method can handle shells with arbitrary curvature, and successfully transitions between membrane and bendingdominated regions, obtaining the expected optimal substructures.
2. Related Work
Our work builds on the ideas from classical structure design for 2D elasticity and plates, with the key ones originating the work of Michell [michell1904l].
We complement these fundamental ideas with quadrangulation techniques which can be reinterpreted as a way to transition from an infinite continuum of fieldaligned beams to its discretization. We use a variation of [bommes2009mixed], but any conforming, fieldaligned method could be used (e.g., [kalberer2007quadcover; Myles:2014; Aigerman:2015; CampenBK15; esck2016]). The optimization method for computing the optimal strain field can be viewed as a specialized crossfield optimization method. Similarly to the recently proposed method of [Knoppel:2015] it has the advantage of being convex.
Shell Optimization
The closest works to ours are the recent works [kilian2017material] and [li2017rib], with which we share a number of ideas. The former describes an elegant connection between curvature and Michell trusses and optimizes the surface shape so principal stress and curvature directions coincide. Only tensile forces are considered, and the volume approximation they use is valid for narrow beams (see Section 3). Similar to our work, [li2017rib] keeps the shell surface fixed. This work considers a network of ribs, aligned with stress lines, and minimizes their volume; similarly to [kilian2017material], this work uses a narrowbeam approximation for the volume, and always produces a thinbeam structure. The crosssection shape of individual beams is optimized, which produces additional weight reduction. We discuss differences to these works in more detail in Section 8.
Our approach also shares similarity with [groen2017homogenization], which uses similar structures for constructing 2D optimized structures.
On the other extreme, [zhao2017stress] considers optimization of variable shell thickness, while keeping the topology fixed, which is suboptimal for bending reinforcement. Our work aims to bridge the gap between these extremes.
The approach of [pietroni2015statics] aligns a network of beams to an input stress field. Another recent related work, [jiang2017design] considers structures made out of beams with a small number of distinct crosssections. Both methods are suitable for architectural design; instead, we focus on applications, like 3D printing, which allow greater flexibility of structures.
Structural Optimization
The literature on structural optimization is quite extensive, and there is no chance that we can do justice to all of it. The main types of approaches found in the literature include topology optimization methods (SIMP or ESObased), analytic methods for optimal structures directly based on Michelltype theories, and methods based on shape derivatives (using an explicit or implicit evolving surface representation). Important books, which include reviews of many other works are [rozvany1976optimal], [allaire2002shape], [bendsoe2013topology], as well as recent reviews, [munk2015topology] and [sigmund2013topology].
The most prevalent methods in topology optimization of structures are based on SIMPtype methods (see [bendsoe2013topology] for a review), which relax the problem to optimizing a density over a domain, which is then converted to a structure by thresholding. This approach has many advantages, including simplicity of implementation [sigmund200199], connection to homogenization theory, flexibility in integrating functionals, and ease of scalable implementation ([Aage2015toptop] and [wu2016system]). Nevertheless, the result will typically depend on initialization: for complex topology to emerge, the domain needs to be discretized using a fine grid. The parameters of the result (e.g., the sparsity of the structure, or minimal thickness) need to be controlled indirectly through algorithm parameters. Finally, the result is a voxelized structure, which then needs to be converted in some way to a form more suitable for manufacturing. In comparison, our method directly produces solutions based on a globally optimal field (in lowvolume limit) and a beam skeleton for the optimized structure, which can be directly adjusted by the user in a variety of ways (e.g., converted to a splinebased CAD model if desired). We compare in more detail in Section 8.
Ground structure methods are among the oldest methods for optimizing topology of truss/beam structures. These methods start with a structure consisting of a large number of redundant beams and optimize it to determine the crosssections, which automatically eliminates some of the beams. Recent examples of applying these type of methods include [sokol201199; zegard2014grand; zegard2015grand3]. Compared to our approach, ground structure methods have to restrict the directions of beams to a small set, which affects both optimality and flexibility of the design. The larger the initial set of beams, the closer they may approximate the optimal result. In computer graphics, the ground structure method was used early in [smith2002creating] for truss structure design. [panetta2015elastic] used a version of a ground structure method to obtain initial topologies for computing microstructures with prescribed material properties, followed by shape optimization.
To a great extent, our work was inspired by the beautiful structures explored in the literature on analytic or semianalytic structure design, e.g. [rozvany2012structural], which includes many examples of exact problem solutions, such as HenckyPrandtl nets. Our goal is to use this type of ideas in the general setting of surface domains, taking advantage of the optimality criteria and insights into the structure of the solutions. A concise exposition of the theory underlying Michelltype optimal layouts can be found in [strang1983hencky]. We note that the application of Michelltype structures in 3D is only appropriate for certain problem settings: e.g., with no lowerbound constraints on shell thickness, variable thickness shells are likely to emerge as a solution [sigmund2016non].
Shapederivative based optimization techniques (e.g., [Allaire2008909]) can obtain very good results when one needs to improve an existing design, by evolving the shape to a local minimum. However, while levelset methods of this type allow for topology changes, the result does vary considerably depending on the starting point. In contrast, our goal is to obtain a starting point that is close to the global optimum, as long as the desired structure has a relatively low volume.
Digital fabrication.
The works closest to ours in this domain are [li2010beam], [tam2015slam], and [tam2015principal]. These methods are based on constructing structures from stress lines on surfaces which, while different from the optimal fields we compute, are often a close approximation. The overall pipeline of the method of [li2010beam] is similar to ours: they start with a field, and construct trusses following the field by tracing lines from supports to loads. The method is limited to two dimensions and demonstrated only for relatively simple structures.
[tam2015slam] uses FDM to add material directly along the principal stress lines, on 2.5D surfaces. The main problems they solve are stress line generation and selection. They minimize strain energy subject to a maximal total print length (i.e max material) and a consistent maximum spacing between lines.
Applying topology optimization (SIMP and ground structure methods) to 3D printing applications is discussed in [zegard2016bridging].
3. Motivating examples
To motivate our method, we start with simple examples of qualitatively different behavior of optimized structures.
The two key behaviors of shelllike optimal structures, observed in specialcase analytic solutions and topology optimization (cf. [sigmund2016non]), are the formation of discrete narrow ribs protruding from the surface in bendingdominated cases (most forces are perpendicular to the surface), and relatively smooth variation in shell thickness in the pure tensile/compression case (inplane forces), as shown in Figure 3.
These behaviors are not observed in the simplified models of beam networks approximating a surface: beams in a typical network do not expand in the direction parallel to the surface, to merge into a variablethickness shell optimal in such cases. It turns out that this is due to the qualitatively inaccurate volume computation with the volume of the beam network approximated by the sum of individual beam volumes. We now consider two simple examples showing why this is the case.
First, we consider optimization of a single horizontal beam of width and height , clamped at one end, and loaded at the other, at an angle to the beam direction (Figure 4, left). This example clarifies optimal behavior when there is stress in only one direction, both for bending () and tension (). The second example involves two intersecting beams (Figure 4, right). It is a simple model for a piece of a surface with there is stress in two directions.
For a single beam, the force along the beam is proportional to the crosssection . The force perpendicular to the beam (bending) is proportional to . The total force is proportional to , where is the beam length, and and are due to projection to the direction of . For simplicity, assume the proportionality constant to be 1. We optimize the beam volume , keeping the force balance constant: . Eliminating , yields unconstrained minimization of , with . Clearly, the solution is to maximize the relative thickness , if : increasing thickness has a higher payoff as forces grow as vs. only for width. However, in the case of a vertical beam, ; the constraint takes the form , i.e., the volume value is fixed and the choice of makes no difference. We conclude that the solution can be always taken to be the ”thick but narrow” beams (we refer to them as ribs), with maximal possible .
The situation is more complicated for surfaces. In the second example, we approximate the surface locally by beams aligned with the perpendicular stress directions (see figure 4 right). For simplicity, we assume the forces, widths , lengths and thicknesses to be the same for both beams. If we view the intersecting beams individually and approximate the total volume as , then the reasoning above applies to each beam: for , only the crosssection matters, but even a small bending component will prioritize maximal solutions, so both beams will be thick and narrow; for inplane forces , we get and the optimal volume no matter which we use.
Considering beams in separation ignores the fact that the intersection area of the beams is counted twice: this part of material is performing “double work”, supporting loads in two directions along two beams. The correct combined volume of two beams is given by
assuming the same beam width and thickness for both. For inplane forces, as before, we have the constraint for constant . The functional can be expressed then as , by eliminating . From the expression it is clear that one would want to maximize the width in this case as opposed to thickness. The optimal volume for max thickness is . In comparison, if we use large thickness , then optimal is close to zero, and the volume is close to , two times higher than optimal .
In the case of two intersecting beams with an outofplane load in addition to inplane, requiring a combination of a tensile and bending force to support, there is an optimal tradeoff between and minimizing the volume, as long as loading cannot be supported by pure bending forces.
If we further impose constraints on maximal and minimal surface thickness, even in tensiledominated areas, beams would form, because the optimal solid shell there would be too thin. A general optimization method should be able to smooth between grillagelike structures for bending dominated areas and “thin and wide” structures for the rest of the surface.
We conclude, from these examples, that to reinforce a shell in a manner close to optimal for arbitrary loads and shell shape, both solid variablethickness and riblike structures may be required in different areas of the surface, and for these to emerge, in a beambased optimization problem, a nonconvex volume function accounting for beam intersection areas has to be used.
4. Problem formulation
We start with a description of the discretization of variablethickness perforated surface shell structure that we use, and the optimization problem we aim to solve.
4.1. Perforated variable thickness shells
Our input is a shell of an initial thickness
(constant per triangle) represented by a triangle mesh, with a vector of external forces
applied to its vertices, and a set of fixed vertices (supports).We aim to compute an optimized shape which we call a perforated shell of variable thickness . consists of (a) a partition of the input surface into polygonal faces (typically quads), corresponding to 3D cells, and (b) an extruded shape for each cell, consisting of blocks as described below.
Given an approximate user target for cell size, our goal is to optimize the edge orientation of the cell boundaries, and thicknesses and width of blocks forming each cell to minimize the weight, while maintaining an upper bound on stresses (calculated using a beam model). can be viewed as the reinforcement structure for .
Cell geometry parametrization.
For each edge of a face of we introduce two parameters, width and thickness . With each edge, we associate a hexahedron (block) constructed by creating a strip inside the face at distance from the edge and extrude the resulting trapezoid along the triangle normal (See Figure 5). While this geometry may result in gaps between adjacent cells, this has no mechanical implications as we treat each block as anchored to edge vertices. We model each side of the cell as a beam including tension/compression and bending forces. While this is a very coarse approximation of the shape, it allows us to obtain an approximation of the solution robustly and quickly. This results can be further refined by shape optimization methods (e.g. [panetta2015elastic]).
Volume discretization.
To simplify our problem, we split all polygonal cells into triangular subcells. We refer to the additional edges inserted in this way as diagonals. We treat these in a special way in the optimization, and in the end ensure that the triangular cells can be merged back into the original polygonal cell (Section 7).
We consider triangular cells with sides , , and blocks with rectangular crosssections of width and thickness along the edge (Figure 5). The simplest approximation of the volume which is typically used in lowvolume truss models, when applied to our system would yield simply . However, as discussed in Section 3 this approximation of the volume results in highly suboptimal results for shells regions dominated by tensile forces.
A more precise approximation is the volume of the extrusion of 3 trapezoidal regions by different heights. Denote , where is the height of the triangle with base on side , for any is the area of the triangle. Let be the permutation of for which , the volume is given by the following simple expression:
(1) 
This volume can be written as where is given by (1) for arbitrary permutation . This expression for the volume is useful for the optimization method in Section 7.
The outofplane heights can be constrained not to exceed a userdefined value , and the normalized widths are constrained so that the trapezoidal areas do not overlap:
(2) 
4.2. Elastic deformation discretization
We model the perforated shell structure as a beam network: for each interior edge, there are two beams, corresponding to the blocks of incident cells along the edge.
Notation.
The beam network consists of a set of beams that are joined together at nodes. For a node , is the set of indices of nodes connected to it, and the vector connecting nodes and corresponds to the edge . For each edge, denotes the unit vector along . The edge has length . We assume that all cells are made of uniform material with as Young modulus. We use and notation for the onedimensional (tensile or bending) stresses and strains of beam connecting vertices in a given cell.
Beam linear elasticity discretization.
The tensile strain along each beam is the scalar , same for both blocks at the edge where is the displacement of a vertex . It is related to the stress by .
For our problem, bending discretization is critical. We use a pure displacementbased beam bending approximation, but a more standard beam element could be used. Our beams are clamped to a freely rotating plane at each vertex, i.e., preserve the angle between the beam and the (freely moving) normal to that plane. We use a simple discrete beam model for bending: we neglect the torsion and shear and define the bending strain, i.e., the change in the curvature of a beam as
(3) 
where and are normals meeting at nodes and (Figure 6), and denotes linearized change of the normal. The normal change, in turn, is expressed in terms of the displacements , of the incident vertices.
This leads to the expression for a scalar bending strain on beams:
(4) 
with the expressions for given in the Appendix.
At a distance from the middle surface of a beam, the strain is given by , where we omit the beam index. Based on the standard Bernoulli beam assumptions after integration over the total beam energy can be expressed as
which leads to the following expression for the stiffness matrix of the beam system:
(5) 
Given the expression for strain , clearly both strain and stress are maximal on one of the surfaces, i.e., for for a given beam. This leads to the following stress constraint:
(6) 
Optimization problem.
We use index for triangular cells. Let is the vector of width parameters of all cells, is the vector of corresponding normalized widths , is the vector of all thickness parameters, the diagonal matrix with thicknesses on the diagonal, and and are displacements and forces respectively. Let be the vector of all ones.
We now formulate our general optimization problem:
(7) 
where the absolute value in the stress constraint is taken elementwise, denotes a vector of ones, and the minimum is taken over all partitions of into triangulated cells. Additionally, we enforce the same thickness on two sides of all cell diagonals.
As optimization over all possible partitions into cells is an intractable problem combining combinatorial and continuous aspects, we use an heuristics to fix the partition first, using
beam continuum approximation (Section 6). Once is fixed, we optimize with respect to and only.While the above formula volume differs only by a seemingly simple quadratic term from the simplest approximation, this completely changes the behavior of the problems, and, in particular, the behavior of the solvers. The problem no longer reduces to convex by a change of variables as it is the case for the simplest formula (cf., e.g., [hemp1973optimum]), and a different type of solvers need to be applied. In our experiments, commonly used general purpose nonconvex solvers converge very slowly and often fail to make progress. Our solution is described in Section 7.
5. Overview of the approach
Our pipeline for solving the optimization problem consists of the following steps.

Field optimization. Compute a pertriangle crossfield on surface using stressbased optimization (Section 6). This field corresponds to an idealized system of densely spaced thin beams (beam continuum
) with directions chosen to minimize weight. The problem is formulated in terms of displacements and the desired crossfield is obtained from the symmetric strain tensors. This requires solving a convex optimization problem with inequality constraints.

Quadrangulation. Create a quaddominant mesh aligned to this crossfield with a usercontrolled spacing of edges. This step is done using a version of mixedinteger quadrangulation [bommes2009mixed], although any quad layout method with field alignment can be used (Section 7). The faces of the mesh will correspond to the cells of the perforated shell .

Cell optimization. Optimize shape parameters of the cells (Section 7). This step has a maximal impact on the final outcome. We introduce a substructure for each cell, with a small number of control parameters defining its shape (widths and thicknesses of rectangular beams along each side). We derive the optimal material distribution by efficiently solving a nonlinear, nonconvex problem minimizing the total volume of all cells with respect to and , while keeping stresses below a userdefined maximum. To make the problem tractable we defined an efficient localglobal optimization method.

Final geometry construction. Finally, we derive the final geometry of according to optimized widths and thickness. We obtain a triangle mesh by performing a sequence of boolean operations between meshes representing the beams. The final watertight and manifold mesh can be directly used for 3D printing or decomposed into cells for FEM analysis.
In the following sections, we describe the details of the steps of the pipeline.
6. Weightminimizing field optimization
In this section, we describe our method for constructing a field of directions on the surface approximating optimal directions for weight minimization with bounded stress. The cells in our construction will be aligned with these directions.
The key idea is to solve a version of the beam weight minimization we assume that there is a continuum of infinitely thin infinitely close beams in two orthogonal directions forming the surface. The density of the beams and their orientations are the optimization variables. The idea of using this type of continua goes back to [michell1904l].
In the simplest case of planar elasticity (which is largely equivalent to the case of a fixed thickness shell in this case) the problem was extensively studied and was shown to reduce to a convex problem. We first describe this classical theory (Michell continua) and then extend it to the case of shells with bending forces (Section 6).
6.1. Michell continua
Here, we briefly review the classical solution, following [strang1983hencky]. The best directions are known to be the principal stress directions of the optimized structure. (These are not the same as the stress direction on the original shell, although these fields are often close.)
The force balance for a plate or a shell with no bending is given by the standard equations in terms of inplane stress tensor , strain , possibly varying elasticity tensor , and external force density :
(8) 
where for a 4tensor and a 2tensor.
A Michell continuum is an idealization of a beam network, characterized, at every point, by beam densities and in two directions, in other words, how many beams cross a unitlength segment along one of the coordinate directions. In the limit of small thicknesses, the total fraction of a small area covered by trusses at a point is . The total volume of the trusses in the continuum can be defined as . Note that this approximation of area covered by trusses suffers from the same flaw pointed out in Section 3.
The optimal trusses have to be oriented along stress directions, and be critically stressed, i.e., all (nonaveraged) stresses on the trusses are equal to maximal stress . This leads to the relationship between and corresponding averaged principal stress: , , where denotes the
th singular value.
Then, we obtain the following optimization problem for the volume, formulated entirely in terms of stresses.
(9) 
where are singular values of the stress tensor. This problem is known to be convex [strang1983hencky]
(although it is more difficult than the linear programming formulation for a truss network). Note that principal stress directions are not fixed and are determined by the optimization. We use these directions as the field for orienting beams in
.The problem (9) has a simple dual (see appendix) of the form
(10) 
Where is the strain of the optimal solution. The dual problem is significantly easier to deal with in the case of continua.
6.2. Continuum optimization with bending
Next we generalize problem (10) to include bending forces.
If the thickness of the shell remains fixed, one can add bending to the functional with relative ease without changing convexity of the problem. We set the shell thickness in this case to half of the maximal allowable thickness; while the resulting field is suboptimal, as we show experimentally in Section 8, inaccuracy in the beam direction has less effect on the overall weight reduction, compared to width/thickness optimization of beams.
We make the standard assumption of planar stress for the shell, i.e., no stresses are active in the direction perpendicular to the shell surface. The strain at a distance from the midline of the shell is given by
where is the bending strain tensor, equal to the linearization of the change in the shape operator (Figure 8)
Consistently with the Michell continuum, we seek to minimize the total weight of a beam continuum, bounding the stress everywhere by
. We observe that the eigenvalues of a
symmetric matrix using the already mentioned substitutions , , , are of the form . Note that these are respectively convex and concave functions of the argument, and, therefore, reach their maxima (respectively minima) on the boundary of the shell, for . For this reason, it is sufficient to bound eigenvalues of stress (or strain) only for and , to guarantee the bounds elsewhere.In the case of bending, the dual problem formulated in terms of displacements has a simpler form relative to the primal problem:
(11) 
Note that now there are two sets of constraints, corresponding to two surfaces of the shell.
The strain tensors and can be interpreted as defining crossfields on the surface.We use the angular average of these fields to align the cells. Finally, we describe a discretization of this convex problem, and how to use the resulting field to build a mesh.
6.3. Discretization and field smoothing
As a first step, we solve a discrete version of the problem (11), which yields displacements
at vertices. From these displacements, we compute the pertriangle strain field eigenvectors, forming a crossfield on the surface, i.e., an assignment of 4 unit vectors, aligned with perpendicular principal strain directions, to each triangle.
This field is only partially defined: we discuss the properties of the solution that allow us to identify areas of the mesh where the directions of the field can be used. In the areas where the field is not defined or is numerically unstable, the beam orientation is not relevant (in the limit of infinite beam density), and we extend the field smoothly. We use a crossfield construction procedure of [bommes2009mixed] to complete the field to the whole surface, identifying, along the way, the singularities of the field.
Discretization of the optimization problem.
The optimization problem (11) has a relatively simple discretization that can be readily plugged in into a cone program solver. e.g., [mosek].
We assume that the surface is given as triangle mesh, , and the same notation is used for edge vectors and vertices as we used for beam networks.
The variables in the problem are displacements, which we discretize using standard piecewiselinear functions on the surface, with the vector of unknowns (we use nonbold letters for highdimensional vectors including all components of corresponding threedimensional quantities).
The two quantities that need to be discretized are tensile and bending strains; we define these per triangle.
If are the vectors along triangle edges, for a triangle , we have the following expression for the strain, computed as :
(12) 
where is the triangle area, is the vertex after in CCW order, is before , and is the vector in the triangle plane perpendicular to the side.
If by abuse of notation, we use to denote the vector of three distinct components of the strain tensor, we can write the expression in the form , where
is the vector of all displacement degrees of freedom.
To discretize the bending strain, we use the trianglebased approximation of the shape operator, following the overall idea in [Onate:1994:STE] and [grinspun2006computing], using vertices of the edgeadjacent triangles to compute the changes of the average normals at edge midpoints, and computing the gradient of the normal using the formulas (12).
This leads to the following expressions for the bending strain on triangles, in which we neglect the triangle deformation: in the deformation modes for which the bending strain is high (i.e., if the curvature changes a lot, inplane deformations are small). We use the following formula [grinspun2006computing], Figure 9:
(13) 
where are linearized changes in the angles between normals of adjacent triangles.
Similarly to , we can write . Then the discrete problem takes the form
(14) 
where , similarly to the beam case, denotes the vector of pervertex forces, and is the vector of all vertex displacements. We use eigenvalue formulas defined for (18) to convert the problem to a convex cone problem, which we solve using the MOSEK solver [mosek].
Detecting field zones.
While the output of the previous step defines a tensor for each triangle, not all of these are meaningful. In some cases (if the triangle is not deformed at all, or deformed negligibly) the strain is zero. More generally, some points may have isotropic strains of the form , where is a nonzero constant, for which all vectors are eigenvectors, so the crossfield is not defined uniquely on this triangle. For general fields, such points are usually isolated. However, for the fields corresponding to the solution of the problem we are considering, the situation is different. There are three possible regimes (see, e.g., [strang1983hencky]). Specifically, the possibilities include

, , principal strains are critical and have opposite directions; this corresponds to welldefined two orthogonal beam families;

, , principal stresses are critical and have the same direction; in this case, beams are not defined;

, ; this corresponds to a single family of beams;

, : in this case, stresses (which are dual variables to the inequality constraints) are both zero, which means there are no beams in this area.
In cases 1 and 3, the crossfield directions are welldefined (purple zones on Fig. 10). In cases 2 and 4 these are either not defined or are not relevant, due to the absence of structure. For this reason, for our construction, we use only regions 1, and 3, which we detect by requiring at least one eigenvalue to be close to , and the difference of eigenvalues to be more than a constant.
We call the resulting field salient. The situation is essentially identical to the crossfields used for constructing quadrangulations: typically, a field aligned with curvature directions is used as a starting point, and only in areas where the difference of two principal curvatures is high.
Completing the field.
To complete the field on the whole surface, which is needed for a complete structure, we use crossfield constrained optimization procedure of [bommes2009mixed]. In this algorithm, the crossfield is encoded by a pertriangle angle, with respect to a reference direction in each triangle. The angles on salient triangles are fixed. On the remaining triangles, these are found by a greedy solve of a mixedinteger problem minimizing the energy
where the summation is over all dual edges connecting triangles and , is the angle between reference directions in triangles, and is an integer unknown accounting for the fact that crossfield values represented by angles are the same.
In the resulting field, defined by the angles for all triangles, and integers for all edges, one can easily compute pervertex field index and identify field singularities, which become irregular (valence different from 4) vertices of the quad mesh at the next step. We refer to [bommes2009mixed] for details of the index computation.
6.4. Construction of the quaddominant network
In general, there may be no optimal beam spacing (in the lowvolume limit, the finer the structure, the lower the optimal volume can be for a given stress). For this reason, the beam spacing is defined by a userspecified parameter . The most direct approach for constructing a quad mesh aligned with a field would be to trace it. However, while it was shown [myles2014robust; ray2014robust] that this approach can be implemented robustly, in general it requires Tjoints (i.e. beams joining other beams in the middle), and it is in general hard to ensure uniform spacing over the whole mesh. We choose a more conservative approach based on constructing a conforming quadrangulation, without Tjoints, using a version of the mixedinteger quadrangulation algorithm [bommes2009mixed] at this step.
While the method does not guarantee perfect alignment of the parametric lines to the input field, it minimizes the deviation in leastsquares sense. We refer to [bommes2009mixed] for further details.
7. Optimizing cell geometry
Given an input quaddominant mesh, split into triangles, we aim at deriving the optimal width and thickness for each edge. As we previously stated, to derive an optimal structure it is important the edges follow the directions we derived in the field optimization step. While our material distribution optimization method works for any mesh, the further the edges deviate from their optimum directions the greater would be the total weight.
7.1. Optimization algorithm
We introduce a domaindecompositionstyle algorithm for solving the problem (7). We observe that in the optimization problem (7), all constraints except are localized, i.e., each constraint uses variables related to one cell. Moreover, the functional itself is a sum of local volume terms . expresses the equation of force balance, i.e., that the sum of beam forces at each node is equal to the external force at this node. Our approach is to fix individual beam forces to their values for current values of geometry parameters and , and then solve for an update to and as a set of local volumeminimization problems with variables , replacing the global constraint with local constraints requiring that individual beam forces remain the same.
We start with an outline and then elaborate on how the local step optimization problems are solved.
Initially, we assign sufficiently high values to and , to ensure that max stress constraints are satisfied.

Global step. The global step is just the standard solve of the elastic equilibrium problem, for fixed cell parameters: Solve , for fixed defined by and . Compute beam forces as described below.

Local step. The local step is the key part of our algorithm. Recall that an important feature of optimal structures is that they are critically stressed i.e., the maximal stress on any element is equal to the maximum possible . The reasons for this are straightforward: if a stress on an element is below zero, one is free to remove some material, increasing the stress on the remaining part, and decreasing the weight.
This motivates our approach. For the local step, we keep the displacements computed at the global step fixed and solve for widths and thicknesses, that would result in maximal stresses on blocks reaching the critical value for given displacements. Each system has 6 unknowns, with 3 constraints on stresses.
Block forces.
To formulate our local algorithm, we introduce block forces, for individual blocks of each cell; we determine and for each cell by minimizing the cell volume, while keeping the block critically stressed i.e., with maximal stress and block forces constant.
Let be the element stiffness matrix corresponding to a block . The vector of forces corresponding to a beam is the vector , i.e. the derivative of the block energy with respect to displacements. Most forces in this vector are zero.
We use lowercase, nonbold for the column vector of , and for column vector of , corresponding to the stresses on block . After some rearrangement, the force vector due to elastic forces produced by a block is
Note that this equation implies that is in the span of vectors and .
Let , be the dual pair of vectors to , . Let the magnitudes of tensile and bending stresses be , ; by taking dot products of both sides with the dual vectors , , we arrive at the equation, where we drop beam/cell subsripts:
Local optimization problem.
The stress in the block, under our assumptions, reaches its maximal value at the top or bottom, where its magnitude is equal to . This leads to the critical stress constraint . Using expressions for and above, which we keep fixed at the local step, this is equivalent to
where we have switched to the variable introduced in (7), where is the corresponding cell triangle height. Without loss of generality, we assume to be 1, which can be achieved by renormalizing all forces to be one. The complete local problem in variables , , is:
(15) 
By eliminating the stresses,we arrive at a single constraint per block relating and , which we express as follows:
(16) 
where , the new variable we introduce to simplify the expressions. This allows us to eliminate all variables from the local optimization problem, leaving only three variables , with constraints , , and .
We say that a cell is filled if the equality is satisfied, i.e. the blocks completely fill the cell.
Without loss of generality, we assume that for the solution ; in practice, 6 problems corresponding to 6 permutations of need to be solved and minimal solution picked.
Proposition 1 ().
The function is a concave function of . As a consequence, its minima are reached on the boundary of the constraint domain; specifically, it is reached at one of the five types of configurations:

all three blocks have maximal thickness: , ;

the cell is filled, i.e. , and no inequality constraint reaches equality;

the cell is filled, and two thicker beams have equal thickness, ;

the cell is filled, and two thinner beams have equal thickness, ;

the cell is filled, and the thickest beam has maximal thickness
In the first case, the solution is completely determined. In the second case, there are four possibilities: no inequality constraint is active (a 2variable unconstrained optimization problem needs to be solved, e.g., parametrized by ); the other three cases define oneparametric families of solutions, and onedimensional unconstrained optimization needs to be performed to find exact values, as we explain below. These families can be parametrized by, e.g., , with the values of the remaining and determined from the active constraints. The proposition is proved in the appendix.
This behavior of is in stark contrast to the lowvolume formulation ignoring common areas of beamlike parts of the structure: one can see that in three cases out of four, it creates a completely filled cell.
Solving the optimization problem.
Proposition 1 leads to an efficient algorithm for the local step.
Observe that the constraint has the form
(17) 
i.e., it is quadratic in . This allows us to reduce the problem to a set of unconstrained optimization problems in one or two variables.

Compute , , , for current displacements.

Evaluate , for the case 1 solution with .

For each permutation of solve three onedimensional optimization problems, minimizing , for each of the cases 35, and the twodimensional problem for case 2 of proposition Prop 1. In each case 35, substitute the active constraint for into (17), yielding a quadratic equation in two remaining free variables, one of which is . Solve it to express the other variable in terms of , and solve a one dimensional optimization problem for . This yields a set of solution candidates; the minimal solution is guaranteed to belong to it.

Pick the minimal solution from the set of solutions obtained for all possible permutations and cases on the previous step.

Update using formula , and recompute the global stiffness matrix .
The convergence behavior of the method is considered in Section 8.
Handling polygonal cells and postprocessing.
There are several factors not considered in the solution method above: (1) possible inconsistency of thickness values across diagonal edges inside triangulated polygonal cells; (2) the coherence of block widths and thicknesses along the edge lines of the quad mesh, approximating the optimized stress lines. While jumps in thickness/width along these lines do not affect the stresses in our simplified model, in practice, these are likely to lead to localized stress concentrations close to jumps, and they are aesthetically objectionable. (3) stress values may slightly exceed the maximal stress after a final global step.
We experimentally observe that many of the candidate solutions have close values, especially in areas with no predominant stress direction.
We address (1) primarily in the process of optimization, at step 4, we pick a minimal candidate solution with lowest block thickness on the diagonal, which may not be the most optimal one, as long as it does not deviate above a threshold. Once the optimization is complete, for each subcell of a polygonal cell we increase the lowest block thickness to the maximal minimal thickness over the whole cell. In addition, for each cell, we store a number of candidate solutions with the smallest volumes.
We address (2) in a postprocessing step using stored candidate solutions: for each nondiagonal edge of a cell, we find its continuation edges along quad mesh edge line in both directions, and choose the candidate solution closest in width to the average of the previous and next edge widths along the edge line.
To address (3), we find all blocks with stress exceeding , and increase their thickness and width, while maintaining constraints, to decrease the stress to the bound. This process is repeated iteratively until convergence. We note that all additional steps are designed to ensure that the final result satisfies stress constraints: in all cases, we never decrease the amount of material in cells, so while the resulting solution may be suboptimal, it always satisfies stress constraints. In practice, the effect of these alterations on the resulting weight is small.
Inflation.
Once we have performed the weight optimization, we have one value of thickness and one value of width for each half edge of the optimized polygonal mesh. If we simply extrude the solid block that matches thickness and width for each half edge cell we end up in the situation illustrated in Figure 11
.a, where there are visible discontinuities between adjacent blocks, causing possible structural discontinuities. Instead, for each continuous stream of quad edges, we generate a unique solid block that interpolates thickness and width along its length (Fig.
11.b). In detail, given a sequence of aligned edges, we first derive a thickness value for each vertex by averaging the thickness from its adjacent half edges. Similarly, we interpolate widths, but this time we derive two different values for each vertex, one for each side of the sequence. We then define a tangent vector for each vertex as the cross product between its normal and its direction along the edge sequence (obtained by averaging the direction of the two incident edges). Then, having defined a proper reference frame, a thickness, and a width for each vertex, we have all the information to extrude a proper volumetric block. We perform a boolean operation to merge all the blocks together to a manifold watertight mesh using the approach of [Zhou:2016].8. Evaluation
Topology optimization.
To validate our approach vs. a generalpurpose topology optimization method, we solve a similar problem with topology optimization code [Aage2015toptop] by restricting the volume of the material to a cylinder of fixed small thickness, and choosing the volume grid resolution to be half of the cylinder thickness, leaving little room for shell shape variation. We observe that for small target volume fractions, as expected, structures emerging in topology optimization are similar to the beam structures we construct (Figure 12), with similar volume fractions. For a large volume fraction, the topology optimization method results in variablethickness shells, but in this case, due to severe constraints on the thickness of the shell, this does not make a significant difference. Default parameters were used in [Aage2015toptop].
Figure 13 shows the results of SIMP topology optimization vs. our method, with compliance as a function of the volume fraction. We observe similar behaviors for both methods. We note that in our case we need to choose the beam spacing parameter: when this is chosen too coarse, the performance will deteriorate.
The role of beam direction.
In Figure 14, we explore the dependence of the role of beam direction in structure optimality, by comparing structure volume for several fields, in addition to our optimized field (Eq. 14). As a “worstcase” baseline we use the crossfield at a maximal distance from the original field (d); as one can see, the field makes a significant difference when it is very far away from the optimal directions, i.e., the choice of directions matters. Similarly, the curvature field (c), unrelated to stress directions, produces relatively high values of the volume. On the other hand, the difference between the optimized field and the stress field (a and b), while present, is relatively minor. This justifies the idea to use, as in, e.g., [li2017rib], the elasticity stress field instead of the optimized field. The former has the additional advantage of better performance and higher smoothness.
Convergence and dependence on the starting point.
We have observed that our algorithm almost invariably converges in several iterations, and yields the expected behavior of solutions in the extreme cases (ribs in the case of bendingdominant shells, and wideandthin beams in hightension/lowbending areas.
The plots show the volumes at each iteration for a basic and more complex problem. Note that sometimes the optimal value is approached from below: the local step overshoots the volume reduction and the stress exceeds the maximal allowed level. Nevertheless, the method recovers reliably.
We also compare our method to two generalpurpose constrained optimization methods, SLSQP [kraft1988software] and Ipopt, a barrier interior point method [byrd2000trust]. In this setup, we used an approximation of the beams volume that is smoother (it uses average instead of max) and simple lower and upper bounds on the width and thickness. As a starting point, we have used the solution of the convex problem with volume ignoring the overlaps. Somewhat surprisingly, these methods were not able to change this initial solution by much after a few hundred iterations; although moving in the right direction, in terms of values, it may differ by a factor up to two from the optimal solution.
While the SLSQP solutions exhibited oscillatory behavior, alternating decreasing the functional with decreasing constraints violations, the interior point method solutions mostly stayed closed to the initial. Figure 16 shows comparative results with Ipopt for a small cantilever test case, we observe that Ipopt converges to a volume almost an order of magnitude larger. When Ipopt is initialized with the solution of the convex problem with simplified volume functional (green), the optimization fails to find a better solution. Using a different initialization (orange: maximum width and thickness) does not provide better results. In contrast, our method converges to a much better solution in only a few iterations.
Effects of constraints on structure parameters.
Figure 17 shows the effect of increasing maximal allowable thickness in a bending scenario, with the structure moving from fully solid, to a structure of narrow but tall beams (if the bound is increased to infinity, in principle, any bending force can be realized by a zerovolume infinitely thin rib).
Similarly, Figure 18 shows the effect of increasing minimal thickness in a tensile scenario. In this case, tall beams appear when increasing the lower bound because widths decrease to compensate for the increase of thickness. The usage of beams in tensile scenarios, while common, is often a consequence of constraints on minimum element size and are suboptimal in the absence of these constraints. On the torsion cylinder experiment, the effect of using a minimal thickness of 0.04 vs. 0.00 increases the volume by 60%.
Relation to related work.
A direct applestoapples comparison of different methods for optimization for shells is extremely difficult for a number of reasons, most importantly, because of different models used (all approaches use a variation of a beam model, and the specific choice has a significant impact on stress estimation). Another important reason includes a strong dependence of the results on the maximal width/thickness/aspect ratio constraints. We summarize qualitative differences from the most closely related works
[kilian2017material] and [li2017rib].The key difference from [kilian2017material] is that it changes the shape of the surface, focuses on the case when bending forces can be neglected, and uses the narrowbeam volume approximation. This may be appropriate for targeted architectural applications but different types of solutions, produced by our method, may be more optimal when these can be practically manufactured. In comparison, we deal with a shell of fixed surface shape, take bending forces into account and use a more precise volume approximation, which leads to a nonconvex problem.
[li2017rib] takes the approach closer to ours, in that it constructs a quadrangulation following a field to determine beam directions and takes bending into account in the elasticity model. However, a narrowbeam approximation for the volume is used, and rather than using an optimized stress field as we do, the stress field of the nonoptimized shell is used instead. In some cases, the optimized field provides an advantage, although the results are typically close as Figure 14 demonstrates. More importantly, the volume approximation has a major effect as shown in Figure 19, where we use one of the models we share with [li2017rib]. In this case, the simpler volume approximation can increase the volume used by up to 44%. On the other hand, [li2017rib] introduces Ibeam optimization of beam crosssections which we have not implemented, and may yield a substantial additional benefit. This technique can be easily combined with ours. Finally, [li2017rib] does additional optimization of discrete beam directions, which can also be added to our method.
Experiments.
We have printed several simple structures to validate our optimization experimentally (Figures 20, 21, 22). Due to the highly approximate nature of the physical model used, we did not attempt a quantitative match to the simulated values, but we did closely match the values of the printed models. The comparison in all cases is between an optimized model and a uniformthickness model of the same weight. Observed displacements in all cases differed by a factor more than 3, suggesting a similar difference in stress.
Finally, Figure 24 shows a set of optimized shells obtained for a variety of shapes using our method; In all cases we have preserved a minimal width/thickness beam to indicate the mesh edges, but the load is carried by a relatively small number of beams; with shell shape constraints we found that unless the thickness bound is set to very low values, even in tension areas ribs tend to appear. This is consistent with the observation that with no thickness, the load carried by bending forces is maximized, if the goal is to reduce the volume.
d  

Aquadom  85.60  7742  0.5  0.005  100  0.2  72.330 
Beetle  169.01  7558  1.2  0.012  50  0.5  35.999 
Boat  92.60  9084  1.0  0.05  100  0.2  33.067 
Botanic  42.59  2152  0.5  0.01  100  0.5  35.586 
Bucket  21.57  7860  0.3  0.03  50  0.5  45.744 
Bunny  106.00  7471  2.0  0.9  100  0.25  22.291 
Duct  460.58  6932  1.0  0.01  100  0.2  69.593 
Leaf  99.24  2980  2.0  0.9  46  0.5  1.954 
Shelf  88.54  2400  4.0  0.5  60  0.5  1.663 
Spoon  133.11  6528  4.0  0.9  29  0.5  4.129 
Stevia  90.68  4840  2.0  0.9  50  0.5  23.941 
Vase  69.64  4798  1.0  0.45  100  0.5  21.817 
Aquadom  16.989  90.335  0.032  1.601 
Beetle  170.128  1287.227  0.103  5.739 
Boat  189.997  418.539  0.160  7.099 
Botanic  9.134  80.533  0.098  1.266 
Bucket  15.028  8.325  0.047  26.924 
Bunny  7270.428  1134.806  1.040  1.229 
Duct  1075.645  19777.817  0.194  2.105 
Leaf  2304.430  659.348  1.158  1.105 
Shelf  1176.000  1386.550  1.090  7.252 
Spoon  1544.960  326.733  1.090  18.621 
Stevia  8329.664  2568.409  1.178  1.198 
Vase  3505.556  492.212  1.026  1.953 
9. Conclusions and Future Work
In this paper, we have described a way to approximate and efficiently solve the problem of minimizing the weight of a support structure for a shell. The proposed method separates the construction into three stages, with the first stage optimizing the field the beam directions must follow and creating a corresponding quaddominant mesh, the second stage creating a cell structure with optimized shape parameters, and the third stage creating an actual realization.
This makes it particularly flexible and allows one to integrate a variety of additional user inputs and constraints, e.g., by modifying the field to change the truss directions, or by adding beams in the second stage to support a connection to a separate object. It is also very efficient, with optimization converging in a few iterations, and quite consistently.
While we use a highly simplified model for cell mechanics, the overall approach admits replacing this model with a more advanced finite element formulation. In this case, it is likely that the local step would require numerical optimization; however, as long as the model for a cell stays lowparametric, one is likely to be able to solve it efficiently.
Limitations.
Our work has two main limitations: first, the formulation for the field optimization still has restrictive lowvolume and fixedthickness assumptions. Based on our evaluation of field direction sensitivity of the final design, we do not view this as a major limitation. The second, more significant, limitation is the highly simplified model we used. If exact results are needed, this model can be used to quickly obtain an initial result, which can then be refined using a more advanced mechanical description of cells and shape optimization.
References
Appendix A Convexity of truss continuum optimization
Expressing the functional in terms of entries of the stress matrix, using , , and , we obtain
(18) 
which verifies convexity of the energy.
Appendix B Dual of the continuum shell problem
This derivation follows [strang1983hencky], which we include here for completeness.
Integrating by parts, and assuming either free boundary or fixed boundary , we obtain
To obtain the dual we need to minimize over all possible in a coordinate system aligned with the expression
The minimum of this expression is , if , for either or ; otherwise, it is zero; this leads to the dual problem (10).
Appendix C Proof of the properties of .
We drop the cell superscript in this proof. We assume that . By the constraints of the problem, , .
(1) is a quadratic function of and linear in ; once we express it in terms of , to eliminate the constraint between and , it becomes a cubic function of , with Hessian
where , and .
A direct evaluation shows that evaluated for is negative for positive and . Similar is true for and . We conclude if for some , any critical point in the interior of the domain is a maximum or a saddle, so there are no minima in the interior. If all are zero, the volume is a linear function of , and the optimum is also on the boundary.
The constraint defines, if for some a quadratic surface in coordinates, and the constraints , , define 3 halfspaces. Their intersection is a tetrahedron with a curved face on the ellipsoid. The faces of the feasible domain correspond to one of the inequality constraints becoming an equality. Similarly, by a direct calculation, substituting either , or into and computing Hessians, we observe that these are not positive definite, hence there can be no solution on faces corresponding to the linear constraint. Finally, the same fact holds for pairs of constraints that define edges of the constraint domain, e.g., and . We conclude that either the solution is on the face corresponding to the constraint or at the only vertex not on this face of the feasible domain, , (case 1 of the proposition). If the minimum satisfies , it is either in the interior (case 2), or on one of the three edges of the boundary of that face (cases 35).
Appendix D Discrete bending strain for beams
The common normal at is defined as
where is the edge following in CCW order around .
At each vertex, we define a normal , as the average of crossproducts of pairs of incident edges. For an edge , incident at a vertex , let be the next CCW edge, and be the previous edge. We define
the (unscaled) normal to the triangle formed by . Then , and .
We express the bending strain in the form , where is the vector of all vertex displacements. Define . Then . The part of affecting is the part perpendicular to :
, where is the matrix . By substitution into the expression for we get:
where is defined as follows:

, if both , are in , i.e. is not on the boundary.

, if .

, if .
For a vector , let be the infinitesimal rotation matrix about : .
Then
where is a matrix.
where , a vector of length 3. From this expression, we can immediately obtain .
Comments
There are no comments yet.