loops.cc – A demonstration of the loop classes
Voro++ has several loop classes that can be used to iterate over a subset of
particles in a container class in a particular order, and carry out
calculations with them. There are three loop classes that can be used with the
standard container classes. The class c_loop_all is the most
frequently used, and loops over all of the particles in a container. In this
example, the other two loop classes are demonstrated:
c_loop_subset can be used to loop over particles in a specific
region of the container, while the c_loop_order class will loop
over a specific order of particles previously stored in a
particle_order class.
The example creates a container class as cube of side length
10. It then adds 100000 particles with random positions to the container.
Before each particle is added to the container, its position is tested to see
whether it lies within a torus. If it does, then the particle is added, and its
position is also stored in the particle_order class (line 37).
Otherwise the particle is added without storing it in the ordering class (line
38). At the end of this procedure, the container has 100000 particles and the
particle_order class has a record of the subset of particles that
lie within this torus.
On line 45, a c_loop_order class is created that is associated
with the particle_order. A call is then made to the
start function, which sets the loop class to consider the first
particle. Usually this returns true, but it may return false in the case when
there are no particles to consider.
If a true result is returned, then a do/while loop is entered. The
compute_cell routine is called to calculate the Voronoi cell for
the current particle in question. This routine returns true if the Voronoi cell
was computed, and false if the cell could not be computed for any reason (such
as being completely removed by a wall). Assuming a true result, the routine
then queries the position of the current particle (on line 49), and then draws
the Voronoi cell as a POV-Ray mesh, and as POV-Ray collection of spheres and
cylinders. The inc routine is then called to advance the loop to
the next particle to be considered. When the inc routine returns
false, there are no more particles to be considered and the routine exits.
The loop therefore calculates Voronoi cells for all the particles that were
previously stored on the particle_order class because they were
within a torus. In the diagram above, this corresponds to the
irregularly-shaped orange torus, comprising of all of the Voronoi cells for
particles that were within the torus. The jagged appearance is due to the fact
that this is just a subset of Voronoi cells from the 100000 in total that fill
the container.
The rest of the code demonstrates another method of looping over a subset of
particles in a container. On line 64, a v_loop_subset class is
created, and on line 65, it is initialized to scan over particles in
rectangular box within the container – this particular rectangular box is
chosen because it is the bounding box for a second torus. The loop then
proceeds in the same way as for the previous case, making a call to
start to find the first particle, and making calls to
inc to find subsequent particles. For each particle, it determines
whether the particle lies within the second torus, and if so, it computes the
Voronoi cell and prints it. This results in the second, blue torus visible in
the figure above.
Code listing
1: // Example code demonstrating the loop classes 2: // 3: // Author : Chris H. Rycroft (LBL / UC Berkeley) 4: // Email : [email protected] 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: }