trunk/examples/interface/loops.cc

   1 // Example code demonstrating the loop classes
   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 // Constants determining the configuration of the tori
  11 const double dis=1.25,mjrad=2.5,mirad=0.95,trad=mjrad+mirad;
  12
  13 // Set the number of particles that are going to be randomly introduced
  14 const int particles=100000;
  15
  16 // This function returns a random double between 0 and 1
  17 double rnd() {return double(rand())/RAND_MAX;}
  18
  19 int main() {
  20         int i;
  21         double x,y,z,r;
  22         voronoicell c;
  23
  24         // Create a container as a non-periodic 10 by 10 by 10 box
  25         container con(-5,5,-5,5,-5,5,26,26,26,false,false,false,8);
  26         particle_order po;
  27
  28         // Randomly add particles into the container
  29         for(i=0;i<particles;i++) {
  30                 x=10*rnd()-5;
  31                 y=10*rnd()-5;
  32                 z=10*rnd()-5;
  33
  34                 // If the particle lies within the first torus, store it in the
  35                 // ordering class when adding to the container
  36                 r=sqrt((x-dis)*(x-dis)+y*y);
  37                 if((r-mjrad)*(r-mjrad)+z*z<mirad) con.put(po,i,x,y,z);
  38                 else con.put(i,x,y,z);
  39         }
  40
  41         // Compute Voronoi cells for the first torus. Here, the points
  42         // previously stored in the ordering class are looped over.
  43         FILE *f1=safe_fopen("loops1_m.pov","w");
  44         FILE *f2=safe_fopen("loops1_v.pov","w");
  45         c_loop_order clo(con,po);
  46         if(clo.start()) do if(con.compute_cell(c,clo)) {
  47
  48                 // Get the position of the current particle under consideration
  49                 clo.pos(x,y,z);
  50
  51                 // Save a POV-Ray mesh to one file and a cylinder/sphere
  52                 // representation to the other file
  53                 c.draw_pov_mesh(x,y,z,f1);
  54                 c.draw_pov(x,y,z,f2);
  55         } while (clo.inc());
  56         fclose(f1);
  57         fclose(f2);
  58
  59         // Compute Voronoi cells for the second torus. Here, the subset loop is
  60         // used to search over the blocks overlapping the torus, and then each
  61         // particle is individually tested.
  62         f1=safe_fopen("loops2_m.pov","w");
  63         f2=safe_fopen("loops2_v.pov","w");
  64         c_loop_subset cls(con);
  65         cls.setup_box(-dis-trad,-dis+trad,-mirad,mirad,-trad,trad,false);
  66         if(cls.start()) do {
  67
  68                 // Get the position of the current particle under consideration
  69                 cls.pos(x,y,z);
  70
  71                 // Test whether this point is within the torus, and if so,
  72                 // compute and save the Voronoi cell
  73                 r=sqrt((x+dis)*(x+dis)+z*z);
  74                 if((r-mjrad)*(r-mjrad)+y*y<mirad&&con.compute_cell(c,cls)) {
  75                         c.draw_pov_mesh(x,y,z,f1);
  76                         c.draw_pov(x,y,z,f2);
  77                 }
  78         } while (cls.inc());
  79         fclose(f1);
  80         fclose(f2);
  81 }