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