random_points.cc – The Voronoi diagram for random points in a cube

Voronoi cells for 20 particles in a cube

This example code demonstrates a basic use of the container class, that is used to hold a particle system prior to the computation of Voronoi cells. In lines 11 to 14 of the code, the boundaries of the container are defined, and the total volume is stored in the variable cvol. For computational efficiency, the container is divided up into a rectangular grid of regions, each of which store the particles in their part of the container. Line 17 of the code creates six blocks in each coordinate, making 216 in total.

On line 31 the container object is created. The first six arguments are the coordinate boundaries. The next three set the computational grid size. The following three set whether each direction is periodic or not; for this example we use non-periodic boundary conditions. The final entry allocates initial space for eight particles per region, although a precise value of this is typically not necessary, as the code can dynamically allocate more.

In lines 36 to 41, twenty particles are added to the container at random positions using the put() routine. This function computes which region the particle is in, and adds its position and a numerical ID (passed in the variable i) to that region.

On line 44, the sum_cell_volumes() routine is called, that calculates all of the Voronoi cells, and sums their volumes. Since the cells should exactly partition the container, this should exactly match the value of the container volume storedcvol. The output of the code is:

Container volume : 8
Voronoi volume   : 8
Difference       : -1.77636e-15

The difference on the order of 10-15 is due to rounding error, as the double-precision floating point numbers used in the calculation typically store around this much accuracy.

In the remaining part of the code, the particle positions and their Voronoi cells are outputted to files that can be read by Gnuplot, using the command:

splot "random_points_p.gnu" u 2:3:4, "random_points_v.gnu" with lines

Code listing

 1: // Voronoi calculation example code
 2: //
 3: // Author   : Chris H. Rycroft (LBL / UC Berkeley)
 4: // Email    : [email protected]
 5: // Date     : August 30th 2011
 6: 
 7: #include "voro++.hh"
 8: using namespace voro;
 9: 
10: // Set up constants for the container geometry
11: const double x_min=-1,x_max=1;
12: const double y_min=-1,y_max=1;
13: const double z_min=-1,z_max=1;
14: const double cvol=(x_max-x_min)*(y_max-y_min)*(x_max-x_min);
15: 
16: // Set up the number of blocks that the container is divided into
17: const int n_x=6,n_y=6,n_z=6;
18: 
19: // Set the number of particles that are going to be randomly introduced
20: const int particles=20;
21: 
22: // This function returns a random double between 0 and 1
23: double rnd() {return double(rand())/RAND_MAX;}
24: 
25: int main() {
26:         int i;
27:         double x,y,z;
28: 
29:         // Create a container with the geometry given above, and make it
30:         // non-periodic in each of the three coordinates. Allocate space for
31:         // eight particles within each computational block
32:         container con(x_min,x_max,y_min,y_max,z_min,z_max,n_x,n_y,n_z,
33:                         false,false,false,8);
34: 
35:         // Randomly add particles into the container
36:         for(i=0;i<particles;i++) {
37:                 x=x_min+rnd()*(x_max-x_min);
38:                 y=y_min+rnd()*(y_max-y_min);
39:                 z=z_min+rnd()*(z_max-z_min);
40:                 con.put(i,x,y,z);
41:         }
42: 
43:         // Sum up the volumes, and check that this matches the container volume
44:         double vvol=con.sum_cell_volumes();
45:         printf("Container volume : %g\n"
46:                "Voronoi volume   : %g\n"
47:                "Difference       : %g\n",cvol,vvol,vvol-cvol);
48: 
49:         // Output the particle positions in gnuplot format
50:         con.draw_particles("random_points_p.gnu");
51: 
52:         // Output the Voronoi cells in gnuplot format
53:         con.draw_cells_gnuplot("random_points_v.gnu");
54: }