1 Introduction
Gridshells, such as steelglass structures, have been used for about forty years in architecture [Otto and Rash 1995]. While trianglebased gridshells seem unbeatable from the point of view of strength, quadbased structures have become popular in the last decade, because of their improved aesthetics and nice mathematical properties. Conversely, there exist fewer studies on more general polygonal structures, most of which are focused on improving mesh geometry for a given topology [Pottmann et al. 2014].
In this paper, we introduce a framework for the generation of gridshell structures that is based on Anisotropic Centroidal Voronoi Tessellations. Our method is driven by the statics of the input surface and it is aimed at improving the strength of the gridshell as well as its aesthetics. Voronoi diagrams appear in nature in many forms, and in several cases they are related to light and strong structures. For instance, bones have a Voronoilike porous structure, with a higher concentration of material where the bone undergoes more stress. This natural principle has also been applied recently to object design for 3D printing [Lu et al. 2014]. We follow a similar approach to design our gridshells, by concentrating more cells of smaller size in zones subject to higher stress, while aligning the elements of our grid to the maximum stress direction.
We start at an input surface and we aim at producing a gridshell that approximates this surface closely. We first perform a static analysis of the surface, from which we obtain an anisotropic, nonEuclidean metric described by the stress tensor. Next we deform the surface, similarly to [Panozzo et al. 2014], in order to transform this anisotropic metric into an Euclidean metric on the deformed surface. We perform Poisson sampling on the deformed surface and we compute a Centroidal Voronoi Tessellation of sampled points. This diagram is mapped back to the original surface to obtain an Anisotropic Centroidal Voronoi Tessellation.
We can control the variation of density and anisotropy of our meshing through two simple parameters. We apply geometric optimization to follow surface symmetries and to improve the local shape of faces of our mesh, making them closer to the faces of Archimedeal solids; this geometric optimization phase greatly contributes to improve the aesthetics of our gridshells, and it also slightly improves its static performances.
We show that the structures generated with our approach, thanks to the great flexibility of Voronoi diagrams, are able to adapt well to the needs of architects and to the designed shapes, while achieving better static behaviour with respect to quadbased gridshells.
2 Related Work
Voronoi diagrams.
The Voronoi Diagram (VD) [Aurenhammer 1991] is a fundamental geometric data structure; within the scope of this work, we discuss only the literature about remeshing techniques based on Centroidal Voronoi Tessellation (CVT) [Du et al. 1999]. VD’s can be defined over 3D surfaces using geodesic distance. This can be done in various ways, either using discrete approximation of geodesic distance [Peyré and Cohen 2006], or using parametrization techniques to bring the problem onto a 2D domain [Alliez et al. 2005]. Valette and Chassery [Valette and Chassery 2004] compute an approximated CVT, based on a discrete global minimization approach. There exist several proposals for computing an Anisotropic CVT (ACVT) under a Riemaniann metric [Du and Wang 2005, Sun et al. 2011, Valette et al. 2008]. Some recent efficient techniques are based on projection of the domain to a 6D space in which the metric becomes Euclidean [Lévy and Bonneel 2012, Zhong et al. 2013]. Panozzo et al. Panozzo:2014 show that a simpler deformation of the surface in 3D is sufficient to get an approximated Euclidean metric, provided that the changes of scale and anisotropy induced by the original metric are not too high.
Architectural geometry.
Most contributions in this field are concerned with the optimization of geometric properties of polygonal meshes approximating a freeform surface. Many works address the planarity of faces, such the construction of PQ (planar quad) meshes [Liu et al. 2006, Liu et al. 2011, Tang et al. 2014, Schiftner and Balzer 2010, Yang et al. 2011, Zadravec et al. 2010], CP (circle packing) meshes [Schiftner et al. 2009], and polygonal hexdominant meshes [Cutler and Whiting 2007, Pottmann et al. 2014, Schiftner et al. 2009, Troche 2008]. Others try to build meshes from a restricted number of tiles or molds [Eigensatz et al. 2010, Fu et al. 2010, Singh and Schaefer 2010, Zimmer et al. 2012]. A few works address the realization of support structures, parallel meshes and torsionfree meshes [Pottmann et al. 2007, Pottmann et al. 2014, Tang et al. 2014]. Among these works, only few focus on the design of a grid topology [Cutler and Whiting 2007, Liu et al. 2011, Schiftner and Balzer 2010, Zadravec et al. 2010] and just Schiftner and Balzer Schiftner:2010 take into account statics. Pottmann et al. Pottmann:2014 mention the possibility of building gridshell structures from either CVT or ACVT to obtain hexdominant meshes; they do not further investigate the underlying design principles, though.
Statics of gridshell structures.
Gridshell structures are a modern response to the ancient need of covering long span spaces. They are compressive structures, i.e. the principal stress comes mainly from axial forces, and this explains the deep interconnection between them and masonry structures. A robust as well as light gridshell can be obtained only through a formfinding process, aimed at finding the funicular surface (surface which stands under compressiononly stresses) that fits the given boundary constraints [Bulenda and Knippers 2001, Ogawa et al. 2008, Otto and Rash 1995]. It is well known [Bulenda and Knippers 2001]
that the form of quad meshes is maintained only if the joints are able to develop bending moments, while triangular meshes maintain their form even if the joints are hinges.
Thrust Network Analysis [Block 2009], a recent formfinding method derived from graphics statics, is specific to masonry. An extension of this method was recently introduced by Tang et al. Tang:2014, which directly allows for gridshell form finding: not only it computes the target funicular surface, but it also optimizes the positions of edges. In Section 5, we compare some of our results with gridshells obtained with this latter method.
The connectivity of the mesh is directly related to the load bearing capacity of the gridshell. While triangular meshes are more rigid and stronger than any other competitor, polygonal meshes have some advantages in terms of ease of construction and lend themselves to the design of torsionfree structures. Some comparative parametric analyses [Malek and Williams 2013] have been carried out about the influence of the remeshing pattern on the gridshell load bearing capacity. There exist surprisingly few studies about the optimal (in a structural sense) connectivity and distribution of edges [Schiftner and Balzer 2010]
, although probably these are – in conjunction with the surface shape – the most influential parameters that govern the structural behavior of the gridshell. For our comparative experiments of different gridshells for a given shape, we adopt the equivalence criterion of simultaneous equal total mass and equal total length of edges, as in
[Malek and Williams 2013].jpg
3 Surface Metric from Static Analysis
The first step of our pipeline is to perform a linear static analysis of the input surface. More precisely we analyze a continuous shell subject to uniform projected load and with all boundary nodes pinrestrained. This analysis returns a tensor field, whose eigenvectors and eigenvalues represent the principal directions and the principal stresses at each point, respectively. In structural mechanics,
isostatic lines are pairs of curves on the shell, which are always tangent to a principal direction, hence always orthogonal to each other. Concentrating the material along the isostatic lines is a good way to improve the structural performance of a structure: in a nutshell, this is what we try to do with our contribution.3.1 Representation
We treat the principal directions and stresses as a double orthogonal line field where and define the minimum and maximum stress at each point of the surface, respectively. Note that only the directions and sizes of are relevant to , not their orientations. Since and are orthogonal, we decouple the scalar and directional information and represent as a triple , where
is a unitlength vector parallel to
, is the maximum stress intensity (henceforth called density), and is the anisotropy (see Figure 1). This representation allows us to better control the influence of over the mesh generation process.3.2 Smoothing
In most cases, the result of static analysis is not directly usable: the signal computed is often irregular with spikes of high stress and abrupt changes of direction in the line field, which are hard to handle during mesh generation. We smooth the line field following the approach of Bommes et al. Bommes:2009, modified as in [Panozzo et al. 2012], see Formula 5. In short, we tradeoff smoothness and faithfulness to the original line field, weighting the second term with anisotropy: we preserve those portions of field where there is a significant difference between the magnitudes of the two stress vectors, while obtaining a smoother field elsewhere.
We also enforce that the two scalar signals and satisfy Lipschitz condition, i.e., , with approximately equal to the diameter of the smallest Voronoi region we expect to obtain. This corresponds to a form of smoothing of the two scalar signals, which is performed through an upper saturation process that preserves the maxima of the function. The results of smoothing are depicted on the lower side of Figure 1.
3.3 Symmetrization
Many architectural models present a few, sometime approximate, symmetry planes that should be preserved in the generated gridshell. Assume we have one or more symmetry planes (shown in red in Fig.2) that partition the mesh into regions. We cross parameterize each symmetry region so that be a crossparametrization that maps a point of region onto its symmetric mate in region . Crossparametrizations are computed between adjacent regions in pairs and propagated about the center of symmetry. For two adjacent regions and , we first crossmap corresponding points on their boundaries, exploiting the common boundary along the symmetry plane, plus symmetric corners that appear along intersections with other planes of symmetry and/or sharp corners on the boundary of the object. Then we compute a harmonic map for each region onto the same parametric domain, in such a way that symmetric points are mapped to the same point in parameter space. Finally, we compute a symmetric field by averaging it componentwise at all the corresponding points in the various regions (see Figure 2).
4 Staticsaware ACVT
Let be a finite set of points, called seeds, sampled in a metric space . The Voronoi Diagram is the partition of into regions such that is the portion of space closer to than to any other seed, with respect to the given metric on . The ACVT is the particular case of VD where the barycenter of each region is coincident with the seed itself. We are interested to the specific case of an ACVT defined over a bounded surface embedded in , with the metric induced by the stress tensor defined in the previous section.
4.1 From general metric to Euclidean metric
We interpret as a frame field and we apply the method described in [Panozzo et al. 2014] to transform the metric induced by into a Euclidean metric on a deformed surface . The metric induced by on is given by symmetric tensor , with
where and are the density and anisotropy described previously, and matrix is expressed at each point in a local coordinate system aligned with . The metric becomes locally Euclidean if the underlying space in the neighborhood of each point is deformed by computed at . We evaluate at each triangle of the input mesh , and we resolve an optimization problem that tends to deform each triangle to its ideal shape to make the metric Euclidean over . See [Panozzo et al. 2014] for further details, and Figure 3 for an example.
The density and anisotropy fields in input may span large intervals which are not always desirable for designing a gridshell. We let the user adjust the desired variation of density and the desired amount of anisotropy over the surface, by introducing two parameters and rescaling the and fields in the intervals and , respectively, prior to computing deformation. We show in Section 5 how such parameters can be used to fine tune the statics as well as the aesthetics of the gridshell.
In order to improve the accuracy of subsequent computations, we refine the input mesh as follows, by subdividing edges that become too elongated under deformation. We set a threshold for the maximum allowed length of an edge of (see next subsection about setting the value of ). After deformation, we split all edges whose length exceeds
, together with their incident triangles, by midpoint subdivision. We estimate
at the centers of new triangles by interpolation, and we deform
back to obtain a refined version of . We iterate between deformation and refinement until all edge lengths are below .4.2 Seed Sampling
We initialize the placement of seeds on by Poisson sampling [Corsini et al. 2012], using a given radius of Poisson disks, which sets a userdefined sampling density, and placing seeds at vertices of . We adapt the refinement of to the desired sampling density by setting
in the previous step. This value has been found experimentally to allow for a rather uniform distribution of seeds and good approximation of the VD.
Since we are dealing with a bordered domain and we want some seeds to remain on the border, we proceed as follows:

We first insert the vertices corresponding to sharp corners on the border of into the set of seeds;

Then we sample just the border of , constrained to the sharp corners;

Finally we sample the interior of , constrained to the seeds inserted in the previous two steps.
4.3 Lloyd relaxation
A CVT is computed through a standard iterative process known as Lloyd relaxation: given a set of seeds, their VD is computed, then each seed is displaced to the centroid of its cell, and the process is repeated until convergence.
We compute a discrete approximated VD, which is sufficient to our purposes: the region of each seed is in fact the collection of vertices of that lie closer to then to the other seeds, according to an approximated geodesic distance computed with the fast method of Campen et al. Campen:2013.
A crucial step of relaxation is the computation of centroids. Given seed and its Voronoi region , at each iteration we choose the vertex inside that minimizes the sum of the squared distances from all the other vertices in the region. Our approach is similar to the one in [Valette and Chassery 2004], but it is based on a simpler, direct, linear computation: for each region we compute the quadric function returning the sum of the squared distances from all the vertices in the region [Garland and Heckbert 1998] and we evaluate for the minimum over .
In order to preserve the boundary, seeds at sharp corners remain still during relaxation; while the other seeds on the boundary are displaced only along the boundary itself, moving a seed each time at the midpoint of its 1D Voronoi region; the boundary is relaxed at each iteration before relaxing the internal seeds.
The CVT is extracted easily from the discrete VD as follows: we set a Voronoi vertex at each triangle of whose vertices belong to three different Voronoi regions, by locating with barycentric coordinates on weighted through the distances of vertices of from their related seeds; and we connect pairs of Voronoi vertices that belong to the border of the same pair of regions. Finally, the ACVT of the original surface is obtained by applying the reverse deformation to the vertices of the CVT of , in order to bring them back to the surface of .
4.4 Regularization
In order to improve the aesthetics, as well as the planarity of faces of the ACVT, we optimize their shape to make them as similar as possible to stretched regular polygons. To this aim, we adopt a framework similar to [Bouaziz et al. 2012], where we alternate perpolygon and pervertex fitting steps.
In the perpolygon step, for each face
, we first perform a Principal Component Analysis to evaluate how much
is stretched with respect to a regular polygon. Then we compute a new polygonal region corresponding to deformed (i.e., unstretched) according to the two lowest rank eigenvectors of the PCA. Next, we define a target regular polygon having the same number of edges and equal perimeter as ; then, using [Besl and McKay 1992], we rigidly align with ; finally, we stretch the oriented polygon back through the reverse deformation that was applied to , and we use the vertices of this stretched regular polygon as target positions to displace the vertices of . Figure 4 shows the steps of this process for a single face.In the pervertex fitting step, for each vertex independently, we evaluate the position minimizing the sum of squared distances from all the target positions specified for by its incident faces. We use a damping factor for improving convergence of this procedure.
An interesting side effect of this regularization procedure is that it tends to make the length of edges more uniform, so that the areas of faces will vary according to the number of sides of polygons. From an aesthetic point of view, this situation matches the look of the Archimedean class of semiregular polyhedra.
Given the similarity of this optimization approach with [Bouaziz et al. 2012], we have also compared our results with the planarization approach presented in that paper. Figure 5 shows our approach in comparison with the initial ACVT and with the result of ShapeUp planarization. As expected, ShapeUp achieves better planar faces, while our algorithm achieves a much better regularity of faces, hence better aesthetics. Planarity is usually measured as distance between diagonals divided by average edge [Tang et al. 2014]
. Unfortunately, this measure is not directly applicable in our case, since diagonals are ambiguous for polygons with an odd number of edges. Then, we generalize the measure of planarity as the average distance of vertices to the best fitting plane divided by half perimeter. Regularity of a face is measured as the sum of squared distance of its vertices to their target positions, divided by its area.
Regularization also slightly improves the overall structural properties of the gridshell structure (10% on average in our experiments). The more uniform length of edges resulting from regularization is also an advantage during production.
4.5 Symmetry
In order to improve the aesthetics of symmetrical structures, we use the same symmetry planes considered in Section 3 and we compute the ACVT just in one of the symmetric sectors. Then we reflect the resulting tessellation to the other sectors, welding them at the regions of seeds placed along symmetry lines. See Figure 6 for a comparison between original ACVT and the optimized and symmetrized tessellation on the Shell dataset.
Non Optimized  Optimized 
5 Results
Our method has been implemented in C++; static analysis has been performed by using the GSA Finite Element Analysis software [Oasys 2014], both on the input surface to obtain the stress tensor, and on the various gridshells to test their behavior. We have tested our method on several surfaces. A summary of the datasets and related results are presented in Table 1.
Overall, the running times for computing a gridshell are negligible with respect to the times required to analyze results. For instance, the models we have analyzed always took less than ten seconds to generate the stress tensor (with GSA); and between one and ten minutes to build the gridshell, depending on the number of iterations in the refinementanddeformation step, which is the bottleneck of the pipeline. While the nonlinear analysis of the result (again with GSA) took over one day for the largest model analyzed.
We present experiments that show the characteristics of our gridshells in terms of statics, as well as some comparisons with gridshells that are obtained with stateoftheart methods in architectural geometry, or correspond to realworld architectures. All structures are assumed to be made of steel, consisting of solid bars with a diameter of 37 millimeters; and the load is distributed uniformly on the whole surface, i.e. each node gets a load that is proportional to the area of its incident faces. We evaluate the following measures, which are most relevant in structural engineering [Meek 1991]:

The nonlinear buckling multiplier , which measures the ability of a structure to support a load equal to a multiple of its weight before collapsing; this attribute measures the robustness of the structure and it should be ideally maximized;

The nodal displacement , which measures the maximum distance with respect to the reference shape when the structure is standing under serviceability load; this attribute measures the degree of deformability of the structure and it should be ideally minimized.
5.1 Tuning parameters
Our method works on three parameters that must be set by the user: the threshold for density , the threshold for anisotropy , and the radius of sampling disks . Comparative tests of static analysis require that different structures have the same weight and total length of beams [Malek and Williams 2013]. We tolerate a 5% of total length variation. In order to keep total length fixed, the number of faces must be decremented as density or anisotropy are incremented, and this is indeed possible by tuning parameter . In the following experiments, we thus take and as free parameters and we set as a constrained variable. Then we test how the variation of density and anisotropy influence the buckling multiplier and nodal displacement.
We have analyzed two funicular surfaces (datasets Shell and Paraboloid) and one light gridshell surface (the dome designed by RFR Paris for the the Abbey of Neumünster in Belgium). We vary density and anisotropy within a range that goes from 1 to 4, with unit step, for a total of 16 test cases per model. Due to the random sampling of seeds, the final meshes may result slightly different even for the same parameters. To disambiguate this randomness, we have performed three experiments for each parameter setting and we have averaged the results of static analysis (so the total number of experiments is in fact 48 per model).
Some pictures illustrating the experiment are shown in Figures 7, 8 and 9 (top views and graphs), while rendering of the best performing models for each dataset are presented in Figure 15.
The Neumünster dataset achieves the highest buckling for and the lowest displacement for , and the latter setting gives the best compromise. This suggests that for this kind of dataset, which is supported on the whole perimeter, density plays a relevant role in improving statics, while anisotropy does not help. The Shell and Paraboloid datasets, on the contrary, rest on a small number of points. The Shell achieves the highest buckling for and the lowest displacement for , and the former setting gives the best compromise. The Paraboloid achieves the highest buckling for and the lowest displacement , but give the best compromise. This suggests that for this class of surfaces variations in both density and anisotropy help improving statics.
Dataset  Model  # Vertices  # Faces  # Edges  Total length (m)  

Aquadom  [Vouga et al. 2012]  1078  1004  2074  3906  1.51  144.90 
Quad (3,3)  1293  1052  2352  3938  2.86  54.04  
Voronoi (3,3)  2382  1177  3752  3898  3.44  48.01  
Botanic  [Tang et al. 2014]  1121  1076  2196  1989  1.00  271.49 
Quad (3,3)  1194  1039  2232  2006  1.80  59.69  
Voronoi (3,3)  2436  1202  3654  2018  1.76  91.59  
British Quad  [Tang et al. 2014]  1648  1568  3216  4286  2.44  19.57 
Quad (2,3)  1987  1585  3583  4145  2.32  24.23  
Voronoi (2,3)  3812  1974  5868  4314  5.78  7.65  
British Tri  Tri (real)  1746  3312  4878  10267  9.62  2.6 
Voronoi (2,3)  3024  5110  15332  10799  6.82  15.09  
Lilium  [Vouga et al. 2012]  1648  636  3216  4286  1.73  61.43 
Quad (3,4)  1987  660  3583  4145  1.86  52.67  
Voronoi (3,4)  3812  695  5868  4314  6.54  16.51  
Neumünster  Voronoi (2,2)  1252  571  1602  975  1.71  10.35 
Paraboloid  Voronoi (4,2)  1424  745  2064  2156  6.59  27.22 
Shell  Voronoi (3,3)  988  477  1441  415  3.05  46.65 









[Tang et al. 2014]  Anisotropic Quad  Voronoi GridShell 
[Tang et al. 2014]  Anisotropic Quad  Voronoi GridShell 
[Vouga et al. 2012]  Anisotropic Quad  Voronoi GridShell 
[Vouga et al. 2012]  Anisotropic Quad  Voronoi GridShell 
5.2 Comparison with Quadrilateral meshes
We compared our gridshells with some quadrilateral meshes obtained with [Tang et al. 2014, Vouga et al. 2012]. As for the previous experiments, we set our parameters to match the total length of edges of the structures we compare with. In order to evaluate how much benefit comes from the Voronoi approach, and how much from allowing for anisotropic and nonuniform meshing, we have also computed anisotropic quadrilateral meshes guided by the same stress tensor that we use for our Voronoi structures. To this aim, we have used the quadrangulation method in [Panozzo et al. 2014] by taking in input (rescaled with the same and parameters we use for the ACVT) as guiding frame field. The experiments are summarized in Table 1 and the related meshes are shown in Figure 10. Our Voronoi gridshells always achieve better performances, in terms of both buckling and displacement, than isotropic quad meshes obtained with stateoftheart methods. Voronoi gridshells have also either better or comparable performances with respect to our anisotropic quad meshes, with the only exception of the Botanic dataset, where the anisotropic quad mesh achieves a smaller displacement. This suggests that both the Voronoi approach and the nonuniform anisotropic meshing play a role in improving performance.
Figure 11
shows the effect of tessellation on the structural behavior of the gridshell. In the Lilium dataset, the forces flow from the top to the restraints along the red paths of structural elements: in our model, such paths are better distributed, thus reducing the elastic strain energy W, as well as the maximal displacement. In the British Quad dataset, almost all the beams of our model undergo the same axial force, whereas in the quad model there is a strong variance of axial forces, including compressions (red) and traction (blue).
Voronoi GridShell  Original 
5.3 Comparison with Triangle meshes
We also compared our structures with real examples of triangulated gridshells. In Figure 12, we show a comparison with the original meshing of the British Museum coverage (dataset British Tri, which is different from the British Quad surface considered in the previous experiment). Related parameters can be found in Table 1. As expected, the triangular mesh has a better static behavior than our structure. The difference in terms of robustness is not dramatic, while the trianglebased gridshell achieves a much better performance in terms of maximum displacement.
By considering this experiment, one may be tempted to deduce that architects should always rely on triangular gridshell structures. However, triangular meshes are considered obsolete nowadays by architects both from an aesthetic and from a manufacturing point of view, while our Voronoi meshes offer an innovative design. More generally, hexdominant structures have several manufacturing advantages: due to the lower valence of nodes, the joints are simpler to manufacture and assemble; besides, it is possible, with further geometric optimization that slightly perturbs the original geometry, to obtain torsionfree structures [Pottmann et al. 2014].
Moreover our patterns have a better perimeter/area ratio, therefore the average size of voronoi panels is significantly lower than the one of triangular meshes for the same total length of beams.
This can be clearly seen in the inset, where the two structures shown in Figure 12 are superimposed. Currently, one of the factors affecting the costs of the shellgrids is the size/radius of the used panels: the larger the size the higher the costs; our structures use smaller panels for the same overall length, allowing for possible economic savings.


5.4 Physical replica
We have fabricated a reduced scale model of the Shell structure composed of 465 joints, 697 beams and 462 panels. The side of this reproduction is 2.4 meters. Each joint has been produced independently using a FDM printer; sticks of wood simulate beams; axternal panels are made of PET (Polyethylene terephthalate) and they have been lasercut. Each component of the structure has a physical label (3D printed on joints, carved by laser on panels or glued paper on sticks), to help us following a map to build the structure. The panels have been fixed by screwing a flat washer at each joint. Some images of the replica are shown in Figure 13.
We have performed load tests on the physical structure, by incrementally applying weights and measuring the displacement of the structure with a proper sensor (see Figure 14). We have have monitored the displacement of the corners of the structure while gradually incrementing the external load. The result of this experiment is shown in the graphs of Figure 14. Obviously we have relied on different materials (wood and ABS, rather then steel), however the general trend of deformation is similar to the simulated mesh.
6 Concluding remarks
We have presented a practical and physically sound framework for the generation of gridshell structures whose topology is based on optimized Anisotropic Centroidal Voronoi Tessellations.
We use the tensor field resulting from FEM stress analysis on the input surface to induce an anisotropic nonEuclidean metric over it. Then we compute an Anisotropic Centroidal Voronoi Tessellation under the same metric. The resulting mesh is hexdominant and made of cells with variable density, depending on the amount of stress, and anisotropic shape, oriented along directions of maximum stress. This mesh is further optimized taking into account symmetry and regularity of cells to improve aesthetics.
We have tested the generated structures evaluating, by means of industrial standard nonlinear analysis simulations, their behavior in terms of nonlinear buckling multiplier and nodal displacement. We have built a reduced scale model and we have performed physical tests on it to verify the soundness of the behavior predicted by the simulation. The result of our experiments demonstrate that our gridshells achieve better static performances with respect to quadbased gridshells, while offering an innovative and aesthetically pleasing look.
References
 [Alliez et al. 2005] Alliez, P., Verdière, E., Devillers, O., and Isenburg, M. 2005. Centroidal Voronoi diagrams for isotropic surface remeshing. Graphical Models 67, 3, 204–231.
 [Aurenhammer 1991] Aurenhammer, F. 1991. Voronoi diagrams: a survey of a fundamental geometric data structure. ACM Computing Surveys (CSUR) 23, 3, 345–405.
 [Besl and McKay 1992] Besl, P. J., and McKay, N. D. 1992. A method for registration of 3D shapes. IEEE transactions on pattern analysis and machine intelligence 14, 2, 239–256.
 [Block 2009] Block, P. 2009. Thrust network analysis: exploring threedimensional equilibrium. PhD thesis, Massachusetts Institute of Technology.
 [Bommes et al. 2009] Bommes, D., Zimmer, H., and Kobbelt, L. 2009. Mixedinteger quadrangulation. ACM Transactions on Graphics 28, 3 (July), 1.
 [Bouaziz et al. 2012] Bouaziz, S., Deuss, M., Schwartzburg, Y., Thibaut, W., and Pauly, M. 2012. ShapeUp: Shaping Discrete Geometry with Projections. Computer Graphics Forum (proc of SGP 2012) 31, 5.
 [Bulenda and Knippers 2001] Bulenda, T., and Knippers, J. 2001. Stability of grid shells. Computers and Structures 79.
 [Campen et al. 2013] Campen, M., Heistermann, M., and Kobbelt, L. 2013. Practical Anisotropic Geodesy. Computer Graphics Forum 32, 5 (Aug.), 63–71.
 [Corsini et al. 2012] Corsini, M., Cignoni, P., and Scopigno, R. 2012. Efficient and flexible sampling with blue noise properties of triangular meshes. IEEE Transactions on Visualization and Computer Graphics 18, 6, 914–924.
 [Cutler and Whiting 2007] Cutler, B., and Whiting, E. 2007. Constrained planar remeshing for architecture. In Graphics Interface, 11–18.
 [Du and Wang 2005] Du, Q., and Wang, D. 2005. Anisotropic centroidal voronoi tessellations and their applications. SIAM J. Sci. Comput. 26, 3 (Mar.), 737–761.
 [Du et al. 1999] Du, Q., Faber, V., and Gunzburger, M. 1999. Centroidal Voronoi Tessellations: Applications and Algorithms. SIAM Review 41, 4 (Jan.), 637–676.
 [Eigensatz et al. 2010] Eigensatz, M., Kilian, M., Schiftner, A., Mitra, N. J., Pottmann, H., and Pauly, M. 2010. Paneling architectural freeform surfaces. ACM Trans. Graph. 29, 4, 45:1–45:10.
 [Fu et al. 2010] Fu, C.W., Lai, C.F., He, Y., and CohenOr, D. 2010. Kset tilable surfaces. ACM Trans. Graph. 29, 4, 44:1–44:6.
 [Garland and Heckbert 1998] Garland, M., and Heckbert, P. 1998. Simplifying surfaces with color and texture using quadric error metrics. In Proc. of the conference on Visualization’98, Ieee, 263–269,.
 [Lévy and Bonneel 2012] Lévy, B., and Bonneel, N. 2012. Variational anisotropic surface meshing with voronoi parallel linear enumeration. In Proceedings of the 21st International Meshing Roundtable, 349–366.
 [Liu et al. 2006] Liu, Y., Pottmann, H., Wallner, J., Yang, Y.L., and Wang, W. 2006. Geometric modeling with conical meshes and developable surfaces. ACM Trans. Graph. 25, 3, 681–689.
 [Liu et al. 2011] Liu, Y., Xu, W., Wang, J., Zhu, L., Guo, B., Chen, F., and Wang, G. 2011. General planar quadrilateral mesh design using conjugate direction field. ACM Trans. Graph. 30, 6, 140:1–140:10.
 [Lu et al. 2014] Lu, L., Sharf, A., Zhao, H., Wei, Y., Fan, Q., Chen, X., Savoye, Y., Tu, C., CohenOr, D., and Chen, B. 2014. Buildtolast: Strength to weight 3d printed objects. ACM Trans. Graph.. to appear.
 [Malek and Williams 2013] Malek, S., and Williams, C. 2013. Structural implications of using cairo tiling and hexagons in gridshells. In Proceedings of the International Association for Shell and Spatial Structures (IASS) Symposium 2013.
 [Meek 1991] Meek, J. L. 1991. Computer Methods in Structural Analysis.
 [Oasys 2014] Oasys, 2014. Gsa analysis. http://http://www.oasyssoftware.com.
 [Ogawa et al. 2008] Ogawa, T., Kato, S., and Fujimoto, M. 2008. Buckling load of elliptic and hyperbolic paraboloidal steel singlelayer reticulated shells of rectangular plan. Journal of the International Association for Shell and Spatial Structures.
 [Otto and Rash 1995] Otto, F., and Rash, B. 1995. Finding Form. Edition Alex Menges, Stuttgart.
 [Panozzo et al. 2012] Panozzo, D., Lipman, Y., Puppo, E., and Zorin, D. 2012. Fields on symmetric surfaces. ACM Trans. Graph. 31, 4, 111:1–111:12.
 [Panozzo et al. 2014] Panozzo, D., Puppo, E., Tarini, M., and SorkineHornung, O. 2014. Frame fields: Anisotropic and nonorthogonal cross fields. ACM Trans. Graph.. to appear.

[Peyré and Cohen 2006]
Peyré, G., and Cohen, L. D.
2006.
Geodesic Remeshing Using Front Propagation.
International Journal of Computer Vision 69
, 1 (May), 145–156.  [Pottmann et al. 2007] Pottmann, H., Liu, Y., Wallner, J., Bobenko, A., and Wang, W. 2007. Geometry of multilayer freeform structures for architecture. ACM Trans. Graphics 26, 3.
 [Pottmann et al. 2014] Pottmann, H., Jiang, C., Höbinger, M., Wang, J., Bompas, P., and Wallner, J. 2014. Cell packing structures. ComputerAided Design (Mar.), 1–14.
 [Schiftner and Balzer 2010] Schiftner, A., and Balzer, J. 2010. Staticssensitive layout of planar quadrilateral meshes. In Advances in Architectural Geometry 2010, C. Ceccato, L. Hesselgren, M. Pauly, H. Pottmann, and J. Wallner, Eds. Springer Vienna, 221–236.
 [Schiftner et al. 2009] Schiftner, A., Höbinger, M., Wallner, J., and Pottmann, H. 2009. Packing circles and spheres on surfaces. ACM Trans. Graph. 28, 5, 139:1–139:8.
 [Singh and Schaefer 2010] Singh, M., and Schaefer, S. 2010. Triangle surfaces with discrete equivalence classes. ACM Trans. Graph. 29, 4, 46:1–46:7.
 [Sun et al. 2011] Sun, F., Choi, Y.K., Wang, W., Yan, D.M., Liu, Y., and Lévy, B. 2011. Obtuse triangle suppression in anisotropic meshes. Comput. Aided Geom. Des. 28, 9 (Dec.), 537–548.
 [Tang et al. 2014] Tang, C., Sun, X., Gomes, A., Wallner, J., and Pottmann, H. 2014. Formfinding with polyhedral meshes made simple. ACM Trans. Graph.. to appear.
 [Troche 2008] Troche, C. 2008. Planar hexagonal meshes by tangent plane intersection. Advances in Architectural Geometry 1, 57–60.
 [Valette and Chassery 2004] Valette, S., and Chassery, J.M. 2004. Approximated Centroidal Voronoi Diagrams for Uniform Polygonal Mesh Coarsening. Computer Graphics Forum 23, 3 (Sept.), 381–389.
 [Valette et al. 2008] Valette, S., Chassery, J. M., and Prost, R. 2008. Generic remeshing of 3d triangular meshes with metricdependent discrete voronoi diagrams. IEEE Transactions on Visualization and Computer Graphics 14, 2 (Mar.), 369–381.
 [Vouga et al. 2012] Vouga, E., Höbinger, M., Wallner, J., and Pottmann, H. 2012. Design of selfsupporting surfaces. ACM Trans. Graph. 31, 4, 87:1–87:11.
 [Yang et al. 2011] Yang, Y.L., Yang, Y.J., Pottmann, H., and Mitra, N. J. 2011. Shape space exploration of constrained meshes. ACM Trans. Graph. 30, 6, 124:1–124:12.
 [Zadravec et al. 2010] Zadravec, M., Schiftner, A., and Wallner, J. 2010. Designing quaddominant meshes with planar faces. Computer Graphics Forum 29, 5, 1671–1679. Proc. Symp. Geometry Processing.
 [Zhong et al. 2013] Zhong, Z., Guo, X., Wang, W., Lévy, B., Sun, F., Liu, Y., and Mao, W. 2013. Particlebased anisotropic surface meshing. ACM Trans. Graph. 32, 4, 99:1–99:14.
 [Zimmer et al. 2012] Zimmer, H., Campen, M., Bommes, D., and Kobbelt, L. 2012. Rationalization of trianglebased pointfolding structures. Comp. Graph. Forum 31, 2pt3, 611–620.