## 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

`splot 'find_voro_cell_v.gnu' w l t 'Voronoi cells', \`

'find_voro_cell.vec' w vec t 'Voronoi cell vectors', \

'find_voro_cell_p.gnu' w p u 2:3:4 t 'Particles'

'find_voro_cell.vec' w vec t 'Voronoi cell vectors', \

'find_voro_cell_p.gnu' w p u 2:3:4 t 'Particles'

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

`set xlabel 'Calculated volume'`

set ylabel 'Sampled volume'

plot [0:0.1] [0:0.1] 'find_voro_cell.vol' u 5:6 with points

set ylabel 'Sampled volume'

plot [0:0.1] [0:0.1] 'find_voro_cell.vol' u 5:6 with points

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 : chr@alum.mit.edu 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: }