Voronoi Grid-Shell Structures

We introduce a framework for the generation of grid-shell structures that is based on Voronoi diagrams and allows us to design tessellations that achieve excellent static performances. We start from an analysis of stress on the input surface and we use the resulting tensor field to induce an anisotropic non-Euclidean metric over it. Then we compute a Centroidal Voronoi Tessellation under the same metric. The resulting mesh is hex-dominant and made of cells with a variable density, which depends on the amount of stress, and anisotropic shape, which depends on the direction of maximum stress. This mesh is further optimized taking into account symmetry and regularity of cells to improve aesthetics. We demonstrate that our grid-shells achieve better static performances with respect to quad-based grid shells, while offering an innovative and aesthetically pleasing look.


page 3

page 4

page 5

page 8

page 9

page 10


Collapse capacity of masonry domes under horizontal loads: A static limit analysis approach

A static limit analysis approach is proposed for assessing the collapse ...

SuperMeshing: A Novel Method for Boosting the Mesh Density in Numerical Computation within 2D Domain

Due to the limit of mesh density, the improvement of the spatial precisi...

A new DG method for a pure–stress formulation of the Brinkman problem with strong symmetry

A strongly symmetric stress approximation is proposed for the Brinkman e...

Stress-constrained topology optimization of lattice-like structures using component-wise reduced order models

Lattice-like structures can provide a combination of high stiffness with...

A Streamline-guided De-Homogenization Approach for Structural Design

We present a novel de-homogenization approach for efficient design of hi...

Computational Design of Cold Bent Glass Façades

Cold bent glass is a promising and cost-efficient method for realizing d...

Correlation structure in the elasticity tensor for short fiber-reinforced composites

The present work provides a profound analytical and numerical analysis o...

1 Introduction

Grid-shells, such as steel-glass structures, have been used for about forty years in architecture [Otto and Rash 1995]. While triangle-based grid-shells seem unbeatable from the point of view of strength, quad-based 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 grid-shell 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 grid-shell 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 Voronoi-like 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 grid-shells, 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 grid-shell that approximates this surface closely. We first perform a static analysis of the surface, from which we obtain an anisotropic, non-Euclidean 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 grid-shells, 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 quad-based grid-shells.

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 free-form 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 hex-dominant 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 torsion-free 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 grid-shell structures from either CVT or ACVT to obtain hex-dominant meshes; they do not further investigate the underlying design principles, though.

Statics of grid-shell structures.

Grid-shell 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 grid-shell can be obtained only through a form-finding process, aimed at finding the funicular surface (surface which stands under compression-only 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 form-finding 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 grid-shell 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 grid-shells obtained with this latter method.

The connectivity of the mesh is directly related to the load bearing capacity of the grid-shell. 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 torsion-free structures. Some comparative parametric analyses [Malek and Williams 2013] have been carried out about the influence of the remeshing pattern on the grid-shell 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 grid-shell. For our comparative experiments of different grid-shells 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].


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 pin-restrained. 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.

Figure 1: Smoothing and saturating the density (left), the anisotropy (middle), and the two orthogonal line fields (right) of the Botanic model: upper side original, lower side smoothed.

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 unit-length 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 trade-off 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 grid-shell. 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 cross-parametrization that maps a point of region onto its symmetric mate in region . Cross-parametrizations are computed between adjacent regions in pairs and propagated about the center of symmetry. For two adjacent regions and , we first cross-map 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 component-wise at all the corresponding points in the various regions (see Figure 2).

Figure 2: Symmetrization of for the Lilium dataset: on the left the cross parametrization defined by two symmetry planes; on the center/right the density field (top) and the line field (bottom), before/after symmetrization.

4 Statics-aware 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.

Figure 3: Density, anisotropy and directional field of the British Quad dataset (left); the resulting deformed domain mesh (top right) and the corresponding undeformed domain (bottom right) with the seeds of the CVT and their distance field.

The density and anisotropy fields in input may span large intervals which are not always desirable for designing a grid-shell. 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 grid-shell.

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

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

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

  3. 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 .

Figure 4: A single face of the ACVT with the eigenvectors resulting form PCA (left); the un-stretched polygon with the aligned target polygon (middle); and the computed displacement vectors in the original space (right).
Figure 5: The effects of the regularization process on planarity (left) and regularity (right). The initial ACVT (top), optimized for planarity with Shape-Up (middle) and regularized using our procedure (bottom).

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 per-polygon and per-vertex fitting steps.

In the per-polygon 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., un-stretched) 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 per-vertex 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 semi-regular 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 Shape-Up planarization. As expected, Shape-Up 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 grid-shell 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
Figure 6: Comparison of non-optimized versus symmetrized and optimized tessellation.

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 grid-shells 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 grid-shell 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 grid-shell, depending on the number of iterations in the refinement-and-deformation step, which is the bottleneck of the pipeline. While the non-linear 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 grid-shells in terms of statics, as well as some comparisons with grid-shells that are obtained with state-of-the-art methods in architectural geometry, or correspond to real-world 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 non-linear 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 grid-shell 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
Table 1: Statistics on datasets and results: for each dataset we show statistics on models taken for comparison and models built by us. Models from [Tang et al. 2014, Vouga et al. 2012] are quad meshes. Note that British Tri and British Quad refer to different surfaces. Quad and Voronoi refer to our models of anisotropic quad meshes and ACVT, respectively, computed with parameters . For each model we report: the number of vertices, faces and edges; the total length of beams in the model; the buckling factor ; and the nodal displacement .

A=1… … A=4

Figure 7: 4x4 test on the Neumünster model.

A=1… … A=4

Figure 8: 4x4 test on the Shell model.

A=1… … A=4

Figure 9: 4x4 test on the Paraboloid model.
[Tang et al. 2014] Anisotropic Quad Voronoi Grid-Shell
[Tang et al. 2014] Anisotropic Quad Voronoi Grid-Shell
[Vouga et al. 2012] Anisotropic Quad Voronoi Grid-Shell
[Vouga et al. 2012] Anisotropic Quad Voronoi Grid-Shell
Figure 10: Comparison with quad meshes, top views. From the top: British Quad, Botanic, Lilium, Aquadom.

5.2 Comparison with Quadrilateral meshes

We compared our grid-shells 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 non-uniform 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 grid-shells always achieve better performances, in terms of both buckling and displacement, than isotropic quad meshes obtained with state-of-the-art methods. Voronoi grid-shells 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 non-uniform anisotropic meshing play a role in improving performance.

Figure 11

shows the effect of tessellation on the structural behavior of the grid-shell. 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).

Figure 11: Comparisons of distribution of AxialForces on the Lilium (top) and British Quad (bottom) datasets. Red corresponds to compression, blue corresponds to traction.
Voronoi Grid-Shell Original
Figure 12: Comparison of our meshing of the British Museum versus the original triangulated mesh.

5.3 Comparison with Triangle meshes

We also compared our structures with real examples of triangulated grid-shells. 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 triangle-based grid-shell 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 grid-shell 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, hex-dominant 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 torsion-free 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 shell-grids 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.

Figure 13: Fabricating a model of the Shell surface.

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 laser-cut. 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.

Figure 14: Top: Setup for load tests over the fabricated model; Bottom: displacement / external load plot.
Figure 15: The models used for parameter tuning tests. From the left: Neumünster, Shell, Paraboloid.
Figure 16: Some of the models used for comparison with quadrilateral meshing. From the top: Aquadom, Botanic and Lilium.

6 Concluding remarks

We have presented a practical and physically sound framework for the generation of grid-shell 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 non-Euclidean metric over it. Then we compute an Anisotropic Centroidal Voronoi Tessellation under the same metric. The resulting mesh is hex-dominant 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 non-linear analysis simulations, their behavior in terms of non-linear 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 grid-shells achieve better static performances with respect to quad-based grid-shells, while offering an innovative and aesthetically pleasing look.


  • [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 3-D shapes. IEEE transactions on pattern analysis and machine intelligence 14, 2, 239–256.
  • [Block 2009] Block, P. 2009. Thrust network analysis: exploring three-dimensional equilibrium. PhD thesis, Massachusetts Institute of Technology.
  • [Bommes et al. 2009] Bommes, D., Zimmer, H., and Kobbelt, L. 2009. Mixed-integer 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. Shape-Up: 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 Cohen-Or, D. 2010. K-set 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., Cohen-Or, D., and Chen, B. 2014. Build-to-last: 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.oasys-software.com.
  • [Ogawa et al. 2008] Ogawa, T., Kato, S., and Fujimoto, M. 2008. Buckling load of elliptic and hyperbolic paraboloidal steel single-layer 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 Sorkine-Hornung, O. 2014. Frame fields: Anisotropic and non-orthogonal 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 multi-layer 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. Computer-Aided Design (Mar.), 1–14.
  • [Schiftner and Balzer 2010] Schiftner, A., and Balzer, J. 2010. Statics-sensitive 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. Form-finding 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 metric-dependent 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 self-supporting 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 quad-dominant 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. Particle-based 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 triangle-based point-folding structures. Comp. Graph. Forum 31, 2pt3, 611–620.