Voro++

find_voro_cell.cc – Demonstrating the find_voronoi_cell routine

Voro++ contains a routine that can take a given position vector, and find the particle whose Voronoi cell the vector lies within. For the regular Voronoi tessellation, this corresponds to finding the closest particle to the vector, since this is how the Voronoi tessellation is defined. For the radical Voronoi tessellation, this corresponds to finding the closest particle, but with an additional weighting due to the particle radii.

This example code demonstrates this functionality in two ways. First, it creates a unit cube and adds twenty particles to it, in the same manner as the random_points.cc example. It then scans a grid of points in the plane z = 0.5, and makes calls to the find_voronoi_cell routine for each. The routine returns the position of the nearest particle, plus its ID. The code saves a file called “find_voro_cell.vec” that contains each of the sample points, plus a displacement vector to the closest particle. Voronoi cells and particles are also outputted. The graph above can be reproduced in Gnuplot using the command

The find_voronoi_cell routine can also be used to estimate the size of a Voronoi cell, which is demonstrated in the second part of the code. The code scans over the unit cube using a grid spacing of h = 0.05, and records how many times the sample points are within each Voronoi cell. Multiplying the number of samples in each Voronoi cell by the cube of h gives an estimate of the cell volumes. It should be noted that when carrying out the sampling, the find_voronoi_cell routine may return false in fails to find a particle – this may occur if the container has no particles or if the sample point lies outside the container geometry. The routine handles periodic containers correctly, and may sometimes return a particle position in a periodic image. However, it does not consider any walls that may have been added to a container.

In the final part of the code, a c_loop_all class is used to loop over all of the particles in the container. It computes their actual Voronoi cell volumes, and saves both these values and the sampled Voronoi volumes to the file “find_voro_cell.vol”. These can be compared in Gnuplot using the command

A close correlation between the two is seen. Decreasing the grid spacing would further improve the accuracy of the sampled volumes.

Code listing

 1: // Example code demonstrating find_voronoi_cell function
 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: // The sampling distance for the grids of find_voronoi_cell calls
11: const double h=0.05;
12: 
13: // The cube of the sampling distance, corresponding the amount of volume
14: // associated with a sample point
15: const double hcube=h*h*h;
16: 
17: // Set the number of particles that are going to be randomly introduced
18: const int particles=20;
19: 
20: // This function returns a random double between 0 and 1
21: double rnd() {return double(rand())/RAND_MAX;}
22: 
23: int main() {
24:         int i;
25:         double x,y,z,r,rx,ry,rz;
26: 
27:         // Create a container with the geometry given above, and make it
28:         // non-periodic in each of the three coordinates. Allocate space for
29:         // eight particles within each computational block
30:         container con(0,1,0,1,0,1,5,5,5,false,false,false,8);
31: 
32:         // Randomly add particles into the container
33:         for(i=0;i<particles;i++) {
34:                 x=rnd();
35:                 y=rnd();
36:                 z=rnd();
37:                 con.put(i,x,y,z);
38:         }
39: 
40:         // Output the particle positions in gnuplot format
41:         con.draw_particles("find_voro_cell_p.gnu");
42: 
43:         // Scan a 2D slice in the container, and for each point in the slice,
44:         // find the Voronoi cell that the point is in. Store a vector
45:         FILE *f1=safe_fopen("find_voro_cell.vec","w");
46:         for(x=0.5*h;x<1;x+=h) for(y=0.5*h;y<1;y+=h) {
47:                 if(con.find_voronoi_cell(x,y,0.5,rx,ry,rz,i))
48:                         fprintf(f1,"%g %g %g %g %g %g %g\n",x,y,0.5,rx-x,ry-y,rz-0.5,
49:                                 sqrt((rx-x)*(rx-x)+(ry-y)*(ry-y)+(rz-0.5)*(rz-0.5)));
50:                 else fprintf(stderr,"# find_voronoi_cell error for %g %g 0.5\n",x,y);
51:         }
52:         fclose(f1);
53: 
54:         // Create a blank array for storing the sampled Voronoi volumes
55:         int samp_v[particles];
56:         for(i=0;i<particles;i++) samp_v[i]=0;
57: 
58:         // Scan over a grid covering the entire container, finding which
59:         // Voronoi cell each point is in, and tallying the result as a method
60:         // of sampling the volume of each Voronoi cell
61:         for(z=0.5*h;z<1;z+=h) for(y=0.5*h;y<1;y+=h) for(x=0.5*h;x<1;x+=h) {
62:                 if(con.find_voronoi_cell(x,y,z,rx,ry,rz,i)) samp_v[i]++;
63:                 else fprintf(stderr,"# find_voronoi_cell error for %g %g %g\n",x,y,z);
64:         }
65: 
66:         // Output the Voronoi cells in gnuplot format and a file with the
67:         // comparisons between the Voronoi cell volumes and the sampled volumes
68:         f1=safe_fopen("find_voro_cell.vol","w");
69:         FILE *f2=safe_fopen("find_voro_cell_v.gnu","w");
70:         c_loop_all cla(con);
71:         voronoicell c;
72:         if(cla.start()) do if (con.compute_cell(c,cla)) {
73: 
74:                 // Get the position and ID information for the particle
75:                 // currently being considered by the loop. Ignore the radius
76:                 // information.
77:                 cla.pos(i,x,y,z,r);
78: 
79:                 // Save and entry to the .vol file, storing both the computed
80:                 // Voronoi cell volume, and the sampled volume based on the
81:                 // number of grid points that were inside the cell
82:                 fprintf(f1,"%d %g %g %g %g %g\n",i,x,y,z,c.volume(),samp_v[i]*hcube);
83: 
84:                 // Draw the Voronoi cell
85:                 c.draw_gnuplot(x,y,z,f2);
86:         } while (cla.inc());
87:         fclose(f1);
88:         fclose(f2);
89: }