branches/exact/examples/no_release/lloyd_box.cc

   1 // Voronoi calculation example code
   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 // Set up constants for the container geometry
  11 const double boxl=1;
  12
  13 // Set up the number of blocks that the container is divided into
  14 const int bl=10;
  15
  16 // Set the number of particles that are going to be randomly introduced
  17 const int particles=4000;
  18
  19 // Set the number of Voronoi faces to bin
  20 const int nface=40;
  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,l;
  27         double x,y,z,r,dx,dy,dz;
  28         int faces[nface],*fp;
  29         double p[3*particles];
  30
  31         // Create a container with the geometry given above, and make it
  32         // non-periodic in each of the three coordinates. Allocate space for
  33         // eight particles within each computational block
  34         container con(-boxl,boxl,-boxl,boxl,-boxl,boxl,bl,bl,bl,false,false,false,8);
  35
  36         // Randomly add particles into the container
  37         for(i=0;i<particles;i++) {
  38                 x=boxl*(2*rnd()-1);
  39                 y=boxl*(2*rnd()-1);
  40                 z=boxl*(2*rnd()-1);
  41                 con.put(i,x,y,z);
  42         }
  43
  44         for(l=0;l<=200;l++) {
  45                 c_loop_all vl(con);
  46                 voronoicell c;
  47                 for(fp=faces;fp<faces+nface;fp++) *fp=0;
  48                 if(vl.start()) do if(con.compute_cell(c,vl)) {
  49                         vl.pos(i,x,y,z,r);
  50                         c.centroid(dx,dy,dz);
  51                         p[3*i]=x+dx;
  52                         p[3*i+1]=y+dy;
  53                         p[3*i+2]=z+dz;
  54
  55                         i=c.number_of_faces()-4;
  56                         if(i<0) i=0;if(i>=nface) i=nface-1;
  57                         faces[i]++;
  58                 } while (vl.inc());
  59                 con.clear();
  60                 for(i=0;i<particles;i++) con.put(i,p[3*i],p[3*i+1],p[3*i+2]);
  61                 printf("%d",l);
  62                 for(fp=faces;fp<faces+nface;fp++) printf(" %d",*fp);
  63                 puts("");
  64         }
  65
  66         // Output the particle positions in gnuplot format
  67         con.draw_particles("sphere_mesh_p.gnu");
  68
  69         // Output the Voronoi cells in gnuplot format
  70         con.draw_cells_gnuplot("sphere_mesh_v.gnu");
  71
  72         // Output the neighbor mesh in gnuplot format
  73         FILE *ff=safe_fopen("sphere_mesh.net","w");
  74         vector<int> vi;
  75         voronoicell_neighbor c;
  76         c_loop_all vl(con);
  77         if(vl.start()) do if(con.compute_cell(c,vl)) {
  78                 i=vl.pid();
  79                 c.neighbors(vi);
  80                 for(l=0;l<(signed int) vi.size();l++) if(vi[l]>i)
  81                         fprintf(ff,"%g %g %g\n%g %g %g\n\n\n",
  82                                 p[3*i],p[3*i+1],p[3*i+2],
  83                                 p[3*vi[l]],p[3*vi[l]+1],p[3*vi[l]+2]);
  84         } while (vl.inc());
  85         fclose(ff);
  86 }