To build Voronoi tessellations, SKIRT uses code taken from the Voro++ library written by Chris H. Rycroft (Harvard University/Lawrence Berkeley Laboratory) at https://github.com/chr1shr/voro (git commit 122531f) with minor adjustments to avoid compiler warnings and to limit the amount of code that needs to be included with SKIRT.
Voro++ is a software library written in C++ for carrying out three-dimensional computations of the Voronoi tessellation. A distinguishing feature of the Voro++ library is that it carries out cell-based calculations, computing the Voronoi cell for each particle individually, rather than computing the Voronoi tessellation as a global network of vertices and edges. In the context of SKIRT, this means that the operation can be easily parallelized.
Refer to the documentation of the VoronoiMeshSnapshot class for more information on how this library is used in SKIRT.
The Voro++ code as used in SKIRT is structured around just a few C++ classes. The cell
class represents a single Voronoi cell as a collection of vertices that are connected by edges, and there are routines for initializing, making, and outputting the cell. The container
class represents a particle system in a three-dimensional rectangular box. It divides this box into a rectangular grid of equally-sized rectangular blocks; this grid is used for computational efficiency during the Voronoi calculations. The compute
class encapsulates all of the routines for actually computing Voronoi cells. Finally, the loop
class allows iterating over the particles in a container.
The cell
class represents a single Voronoi cell as a convex polyhedron, with a set of vertices that are connected by edges. The class contains a variety of functions that can be used to compute and output the Voronoi cell corresponding to a particular particle. The command init() can be used to initialize a cell as a large rectangular box. The Voronoi cell can then be computed by repeatedly cutting it with planes that correspond to the perpendicular bisectors between that particle and its neighbors.
This is achieved by using the plane() routine, which will recompute the cell's vertices and edges after cutting it with a single plane. This is the key routine in the cell
class. It begins by exploiting the convexity of the underlying cell, tracing between edges to work out if the cell intersects the cutting plane. If it does not intersect, then the routine immediately exits. Otherwise, it finds an edge or vertex that intersects the plane, and from there, traces out a new face on the cell, recomputing the edge and vertex structure accordingly.
Once the cell is computed, there are several routines for computing features of the Voronoi cell, such as its volume or centroid.
The polyhedral structure of the cell is stored in the following arrays:
In a very large number of cases, the values of n_i will be 3. This is because the only way that a higher-order vertex can be created in the plane() routine is if the cutting plane perfectly intersects an existing vertex. For random particle arrangements with position vectors specified to double precision this should happen very rarely. The situation is different for cases featuring crystalline arrangements of particles where the corresponding Voronoi cells may have high-order vertices by construction.
Because of this, Voro++ takes the approach that it if an existing vertex is within a small numerical tolerance of the cutting plane, it is treated as being exactly on the plane, and the polyhedral topology is recomputed accordingly. However, while this improves robustness, it also adds the complexity that n_i may no longer always be 3. This causes memory management to be significantly more complicated, as different vertices require a different number of elements in the ed[][] array. To accommodate this, the cell
class allocates edge memory in a different array called mep[][], in such a way that all vertices of order k are held in mep[k]. If vertex i has order k, then ed[i] points to memory within mep[k]. The array ed[][] is never directly initialized as a two-dimensional array itself, but points at allocations within mep[][]. To the user, it appears as though each row of ed[][] has a different number of elements. When vertices are added or deleted, care must be taken to reorder and reassign elements in these arrays.
During the plane() routine, the code traces around the vertices of the cell, and adds new vertices along edges which intersect the cutting plane to create a new face. The values of l(j,i) are used in this computation, as when the code is traversing from one vertex on the cell to another, this information allows the code to immediately work out which edge of a vertex points back to the one it came from. As new vertices are created, the l(j,i) are also updated to ensure consistency. To ensure robustness, the plane cutting algorithm should work with any possible combination of vertices which are inside, outside, or exactly on the cutting plane.
Vertices exactly on the cutting plane create some additional computational difficulties. If there are two marginal vertices connected by an existing edge, then it would be possible for duplicate edges to be created between those two vertices, if the plane routine traces along both sides of this edge while constructing the new face. The code recognizes these cases and prevents the double edge from being formed. Another possibility is the formation of vertices of order two or one. At the end of the plane cutting routine, the code checks to see if any of these are present, removing the order one vertices by just deleting them, and removing the order two vertices by connecting the two neighbors of each vertex together. It is possible that the removal of a single low-order vertex could result in the creation of additional low-order vertices, so the process is applied recursively until no more are left.
The compute
class encapsulates the routines for carrying out the Voronoi cell computations. It contains data structures suchs as a mask and a queue that are used in the computations. During the computation, the class calls routines in the container class to access the particle positions that are stored there.
The key routine in this class is compute_cell(), which makes use of a cell
class to construct a Voronoi cell for a specific particle in the container. The basic approach that this function takes is to repeatedly cut the Voronoi cell by planes corresponding neighboring particles, and stop when it recognizes that all the remaining particles in the container are too far away to possibly influence cell's shape. The code makes use of two possible methods for working out when a cell computation is complete:
Another useful observation is that the regions that need to be tested are simply connected, meaning that if a particular region does not need to be tested, then neighboring regions which are further away do not need to be tested.
For maximum efficiency, it was found that a hybrid approach making use of both of the above tests works well in practice. Radius tests work well for the first few blocks, but switching to region tests after then prevents the code from becoming extremely slow, due to testing over very large spherical shells of particles. The compute_cell() routine therefore takes the following approach: