Voro++
voro++.hh
Go to the documentation of this file.
1 // Voro++, a 3D cell-based Voronoi library
2 //
3 // Author : Chris H. Rycroft (LBL / UC Berkeley)
4 // Email : chr@alum.mit.edu
5 // Date : August 30th 2011
6 
7 /** \file voro++.hh
8  * \brief A file that loads all of the Voro++ header files. */
9 
10 /** \mainpage Voro++ class reference manual
11  * \section intro Introduction
12  * Voro++ is a software library for carrying out three-dimensional computations
13  * of the Voronoi tessellation. A distinguishing feature of the Voro++ library
14  * is that it carries out cell-based calculations, computing the Voronoi cell
15  * for each particle individually, rather than computing the Voronoi
16  * tessellation as a global network of vertices and edges. It is particularly
17  * well-suited for applications that rely on cell-based statistics, where
18  * features of Voronoi cells (eg. volume, centroid, number of faces) can be
19  * used to analyze a system of particles.
20  *
21  * Voro++ is written in C++ and can be built as a static library that can be
22  * linked to. This manual provides a reference for every function in the class
23  * structure. For a general overview of the program, see the Voro++ website at
24  * http://math.lbl.gov/voro++/ and in particular the example programs at
25  * http://math.lbl.gov/voro++/examples/ that demonstrate many of the library's
26  * features.
27  *
28  * \section class C++ class structure
29  * The code is structured around several C++ classes. The voronoicell_base
30  * class contains all of the routines for constructing a single Voronoi cell.
31  * It represents the cell as a collection of vertices that are connected by
32  * edges, and there are routines for initializing, making, and outputting the
33  * cell. The voronoicell_base class form the base of the voronoicell and
34  * voronoicell_neighbor classes, which add specialized routines depending on
35  * whether neighboring particle ID information for each face must be tracked or
36  * not. Collectively, these classes are referred to as "voronoicell classes"
37  * within the documentation.
38  *
39  * There is a hierarchy of classes that represent three-dimensional particle
40  * systems. All of these are derived from the voro_base class, which contains
41  * constants that divide a three-dimensional system into a rectangular grid of
42  * equally-sized rectangular blocks; this grid is used for computational
43  * efficiency during the Voronoi calculations.
44  *
45  * The container_base, container, and container_poly are then derived from the
46  * voro_base class to represent a particle system in a specific
47  * three-dimensional rectangular box using both periodic and non-periodic
48  * boundary conditions. In addition, the container_periodic_base,
49  * container_periodic, and container_periodic_poly classes represent
50  * a particle system in a three-dimensional non-orthogonal periodic domain,
51  * defined by three periodicity vectors that represent a parallelepiped.
52  * Collectively, these classes are referred to as "container classes" within
53  * the documentation.
54  *
55  * The voro_compute template encapsulates all of the routines for computing
56  * Voronoi cells. Each container class has a voro_compute template within
57  * it, that accesses the container's particle system, and computes the Voronoi
58  * cells.
59  *
60  * There are several wall classes that can be used to apply certain boundary
61  * conditions using additional plane cuts during the Voronoi cell compution.
62  * The code also contains a number of small loop classes, c_loop_all,
63  * c_loop_subset, c_loop_all_periodic, and c_loop_order that can be used to
64  * iterate over a certain subset of particles in a container. The latter class
65  * makes use of a special particle_order class that stores a specific order of
66  * particles within the container. The library also contains the classes
67  * pre_container_base, pre_container, and pre_container_poly, that can be used
68  * as temporary storage when importing data of unknown size.
69  *
70  * \section voronoicell The voronoicell classes
71  * The voronoicell class represents a single Voronoi cell as a convex
72  * polyhedron, with a set of vertices that are connected by edges. The class
73  * contains a variety of functions that can be used to compute and output the
74  * Voronoi cell corresponding to a particular particle. The command init()
75  * can be used to initialize a cell as a large rectangular box. The Voronoi cell
76  * can then be computed by repeatedly cutting it with planes that correspond to
77  * the perpendicular bisectors between that particle and its neighbors.
78  *
79  * This is achieved by using the plane() routine, which will recompute the
80  * cell's vertices and edges after cutting it with a single plane. This is the
81  * key routine in voronoicell class. It begins by exploiting the convexity
82  * of the underlying cell, tracing between edges to work out if the cell
83  * intersects the cutting plane. If it does not intersect, then the routine
84  * immediately exits. Otherwise, it finds an edge or vertex that intersects
85  * the plane, and from there, traces out a new face on the cell, recomputing
86  * the edge and vertex structure accordingly.
87  *
88  * Once the cell is computed, there are many routines for computing features of
89  * the the Voronoi cell, such as its volume, surface area, or centroid. There
90  * are also many routines for outputting features of the Voronoi cell, or
91  * writing its shape in formats that can be read by Gnuplot or POV-Ray.
92  *
93  * \subsection internal Internal data representation
94  * The voronoicell class has a public member p representing the
95  * number of vertices. The polyhedral structure of the cell is stored
96  * in the following arrays:
97  *
98  * - pts: a one-dimensional array of floating point numbers, that represent the
99  * position vectors x_0, x_1, ..., x_{p-1} of the polyhedron vertices.
100  * - nu: the order of each vertex n_0, n_1, ..., n_{p-1}, corresponding to
101  * the number of other vertices to which each is connected.
102  * - ed: a two-dimensional table of edges and relations. For the ith vertex,
103  * ed[i] has 2n_i+1 elements. The first n_i elements are the edges e(j,i),
104  * where e(j,i) is the jth neighbor of vertex i. The edges are ordered
105  * according to a right-hand rule with respect to an outward-pointing normal.
106  * The next n_i elements are the relations l(j,i) which satisfy the property
107  * e(l(j,i),e(j,i)) = i. The final element of the ed[i] list is a back
108  * pointer used in memory allocation.
109  *
110  * In a very large number of cases, the values of n_i will be 3. This is because
111  * the only way that a higher-order vertex can be created in the plane()
112  * routine is if the cutting plane perfectly intersects an existing vertex. For
113  * random particle arrangements with position vectors specified to double
114  * precision this should happen very rarely. A preliminary version of this code
115  * was quite successful with only making use of vertices of order 3. However,
116  * when calculating millions of cells, it was found that this approach is not
117  * robust, since a single floating point error can invalidate the computation.
118  * This can also be a problem for cases featuring crystalline arrangements of
119  * particles where the corresponding Voronoi cells may have high-order vertices
120  * by construction.
121  *
122  * Because of this, Voro++ takes the approach that it if an existing vertex is
123  * within a small numerical tolerance of the cutting plane, it is treated as
124  * being exactly on the plane, and the polyhedral topology is recomputed
125  * accordingly. However, while this improves robustness, it also adds the
126  * complexity that n_i may no longer always be 3. This causes memory management
127  * to be significantly more complicated, as different vertices require a
128  * different number of elements in the ed[][] array. To accommodate this, the
129  * voronoicell class allocated edge memory in a different array called mep[][],
130  * in such a way that all vertices of order k are held in mep[k]. If vertex
131  * i has order k, then ed[i] points to memory within mep[k]. The array ed[][]
132  * is never directly initialized as a two-dimensional array itself, but points
133  * at allocations within mep[][]. To the user, it appears as though each row of
134  * ed[][] has a different number of elements. When vertices are added or
135  * deleted, care must be taken to reorder and reassign elements in these
136  * arrays.
137  *
138  * During the plane() routine, the code traces around the vertices of the cell,
139  * and adds new vertices along edges which intersect the cutting plane to
140  * create a new face. The values of l(j,i) are used in this computation, as
141  * when the code is traversing from one vertex on the cell to another, this
142  * information allows the code to immediately work out which edge of a vertex
143  * points back to the one it came from. As new vertices are created, the l(j,i)
144  * are also updated to ensure consistency. To ensure robustness, the plane
145  * cutting algorithm should work with any possible combination of vertices
146  * which are inside, outside, or exactly on the cutting plane.
147  *
148  * Vertices exactly on the cutting plane create some additional computational
149  * difficulties. If there are two marginal vertices connected by an existing
150  * edge, then it would be possible for duplicate edges to be created between
151  * those two vertices, if the plane routine traces along both sides of this
152  * edge while constructing the new face. The code recognizes these cases and
153  * prevents the double edge from being formed. Another possibility is the
154  * formation of vertices of order two or one. At the end of the plane cutting
155  * routine, the code checks to see if any of these are present, removing the
156  * order one vertices by just deleting them, and removing the order two
157  * vertices by connecting the two neighbors of each vertex together. It is
158  * possible that the removal of a single low-order vertex could result in the
159  * creation of additional low-order vertices, so the process is applied
160  * recursively until no more are left.
161  *
162  * \section container The container classes
163  * There are four container classes available for general usage: container,
164  * container_poly, container_periodic, and container_periodic_poly. Each of
165  * these represent a system of particles in a specific three-dimensional
166  * geometry. They contain routines for importing particles from a text file,
167  * and adding particles individually. They also contain a large number of
168  * analyzing and outputting the particle system. Internally, the routines that
169  * compute Voronoi cells do so by making use of the voro_compute template.
170  * Each container class contains routines that tell the voro_compute template
171  * about the specific geometry of this container.
172  *
173  * \section voro_compute The voro_compute template
174  * The voro_compute template encapsulates the routines for carrying out the
175  * Voronoi cell computations. It contains data structures suchs as a mask and a
176  * queue that are used in the computations. The voro_compute template is
177  * associated with a specific container class, and during the computation, it
178  * calls routines in the container class to access the particle positions that
179  * are stored there.
180  *
181  * The key routine in this class is compute_cell(), which makes use of a
182  * voronoicell class to construct a Voronoi cell for a specific particle in the
183  * container. The basic approach that this function takes is to repeatedly cut
184  * the Voronoi cell by planes corresponding neighboring particles, and stop
185  * when it recognizes that all the remaining particles in the container are too
186  * far away to possibly influence cell's shape. The code makes use of two
187  * possible methods for working out when a cell computation is complete:
188  *
189  * - Radius test: if the maximum distance of a Voronoi cell
190  * vertex from the cell center is R, then no particles more than a distance
191  * 2R away can possibly influence the cell. This a very fast computation to
192  * do, but it has no directionality: if the cell extends a long way in one
193  * direction then particles a long distance in other directions will still
194  * need to be tested.
195  * - Region test: it is possible to test whether a specific region can
196  * possibly influence the cell by applying a series of plane tests at the
197  * point on the region which is closest to the Voronoi cell center. This is a
198  * slower computation to do, but it has directionality.
199  *
200  * Another useful observation is that the regions that need to be tested are
201  * simply connected, meaning that if a particular region does not need to be
202  * tested, then neighboring regions which are further away do not need to be
203  * tested.
204  *
205  * For maximum efficiency, it was found that a hybrid approach making use of
206  * both of the above tests worked well in practice. Radius tests work well for
207  * the first few blocks, but switching to region tests after then prevent the
208  * code from becoming extremely slow, due to testing over very large spherical
209  * shells of particles. The compute_cell() routine therefore takes the
210  * following approach:
211  *
212  * - Initialize the voronoicell class to fill the entire computational domain.
213  * - Cut the cell by any wall objects that have been added to the container.
214  * - Apply plane cuts to the cell corresponding to the other particles which
215  * are within the current particle's region.
216  * - Test over a pre-computed worklist of neighboring regions, that have been
217  * ordered according to the minimum distance away from the particle's
218  * position. Apply radius tests after every few regions to see if the
219  * calculation can terminate.
220  * - If the code reaches the end of the worklist, add all the neighboring
221  * regions to a new list.
222  * - Carry out a region test on the first item of the list. If the region needs
223  * to be tested, apply the plane() routine for all of its particles, and then
224  * add any neighboring regions to the end of the list that need to be tested.
225  * Continue until the list has no elements left.
226  *
227  * The compute_cell() routine forms the basis of many other routines, such as
228  * store_cell_volumes() and draw_cells_gnuplot() that can be used to calculate
229  * and draw the cells in a container.
230  *
231  * \section walls Wall computation
232  * Wall computations are handled by making use of a pure virtual wall class.
233  * Specific wall types are derived from this class, and require the
234  * specification of two routines: point_inside() that tests to see if a point
235  * is inside a wall or not, and cut_cell() that cuts a cell according to the
236  * wall's position. The walls can be added to the container using the
237  * add_wall() command, and these are called each time a compute_cell() command
238  * is carried out. At present, wall types for planes, spheres, cylinders, and
239  * cones are provided, although custom walls can be added by creating new
240  * classes derived from the pure virtual class. Currently all wall types
241  * approximate the wall surface with a single plane, which produces some small
242  * errors, but generally gives good results for dense particle packings in
243  * direct contact with a wall surface. It would be possible to create more
244  * accurate walls by making cut_cell() routines that approximate the curved
245  * surface with multiple plane cuts.
246  *
247  * The wall objects can used for periodic calculations, although to obtain
248  * valid results, the walls should also be periodic as well. For example, in a
249  * domain that is periodic in the x direction, a cylinder aligned along the x
250  * axis could be added. At present, the interior of all wall objects are convex
251  * domains, and consequently any superposition of them will be a convex domain
252  * also. Carrying out computations in non-convex domains poses some problems,
253  * since this could theoretically lead to non-convex Voronoi cells, which the
254  * internal data representation of the voronoicell class does not support. For
255  * non-convex cases where the wall surfaces feature just a small amount of
256  * negative curvature (eg. a torus) approximating the curved surface with a
257  * single plane cut may give an acceptable level of accuracy. For non-convex
258  * cases that feature internal angles, the best strategy may be to decompose
259  * the domain into several convex subdomains, carry out a calculation in each,
260  * and then add the results together. The voronoicell class cannot be easily
261  * modified to handle non-convex cells as this would fundamentally alter the
262  * algorithms that it uses, and cases could arise where a single plane cut
263  * could create several new faces as opposed to just one.
264  *
265  * \section loops Loop classes
266  * The container classes have a number of simple routines for calculating
267  * Voronoi cells for all particles within them. However, in some situations it
268  * is desirable to iterate over a specific subset of particles. This can be
269  * achieved with the c_loop classes that are all derived from the c_loop_base
270  * class. Each class can iterate over a specific subset of particles in a
271  * container. There are three loop classes for use with the container and
272  * container_poly classes:
273  *
274  * - c_loop_all will loop over all of the particles in a container.
275  * - c_loop_subset will loop over a subset of particles in a container that lie
276  * within some geometrical region. It can loop over particles in a
277  * rectangular box, particles in a sphere, or particles that lie within
278  * specific internal computational blocks.
279  * - c_loop_order will loop over a specific list of particles that were
280  * previously stored in a particle_order class.
281  *
282  * Several of the key routines within the container classes (such as
283  * draw_cells_gnuplot and print_custom) have versions where they can be passed
284  * a loop class to use. Loop classes can also be used directly and there are
285  * some examples on the library website that demonstrate this. It is also
286  * possible to write custom loop classes.
287  *
288  * In addition to the loop classes mentioned above, there is also a
289  * c_loop_all_periodic class, that is specifically for use with the
290  * container_periodic and container_periodic_poly classes. Since the data
291  * structures of these containers differ considerably, it requires a different
292  * loop class that is not interoperable with the others.
293  *
294  * \section pre_container The pre_container classes
295  * Voro++ makes use of internal computational grid of blocks that are used to
296  * configure the code for maximum efficiency. As discussed on the library
297  * website, the best performance is achieved for around 5 particles per block,
298  * with anything in the range from 3 to 12 giving good performance. Usually
299  * the size of the grid can be chosen by ensuring that the number of blocks is
300  * equal to the number of particles divided by 5.
301  *
302  * However, this can be difficult to choose in cases when the number of
303  * particles is not known a priori, and in thes cases the pre_container classes
304  * can be used. They can import an arbitrary number of particle positions from
305  * a file, dynamically allocating memory in chunks as necessary. Once particles
306  * are imported, they can guess an optimal block arrangement to use for the
307  * container class, and then transfer the particles to the container. By
308  * default, this procedure is used by the command-line utility to enable it to
309  * work well with arbitrary sizes of input data.
310  *
311  * The pre_container class can be used when no particle radius information is
312  * available, and the pre_container_poly class can be used when radius
313  * information is available. At present, the pre_container classes can only be
314  * used with the container and container_poly classes. They do not support
315  * the container_periodic and container_periodic_poly classes. */
316 
317 #ifndef VOROPP_HH
318 #define VOROPP_HH
319 
320 #include "config.hh"
321 #include "common.hh"
322 #include "cell.hh"
323 #include "v_base.hh"
324 #include "rad_option.hh"
325 #include "container.hh"
326 #include "unitcell.hh"
327 #include "container_prd.hh"
328 #include "pre_container.hh"
329 #include "v_compute.hh"
330 #include "c_loops.hh"
331 #include "wall.hh"
332 
333 #endif
Header file for the classes encapsulating functionality for the regular and radical Voronoi tessellat...
Header file for the derived wall classes.
Header file for the loop classes.
Header file for the voronoicell and related classes.
Master configuration file for setting various compile-time options.
Header file for the voro_compute template and related classes.
Header file for the unitcell class.
Header file for the container_periodic_base and related classes.
Header file for the container_base and related classes.
Header file for the pre_container and related classes.
Header file for the base Voronoi container class.
Header file for the small helper functions.