trunk/examples/interface/find_voro_cell.cc

   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 }