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
rnd()21 double rnd() {return double(rand())/RAND_MAX;}
22
main()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 }
90