loops.cc – A demonstration of the loop classes

Two tori demonstrating the loop classes in Voro++

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    : chr@alum.mit.edu
 5: // Date     : August 30th 2011
 7: #include "voro++.hh"
 8: using namespace voro;
10: // Constants determining the configuration of the tori
11: const double dis=1.25,mjrad=2.5,mirad=0.95,trad=mjrad+mirad;
13: // Set the number of particles that are going to be randomly introduced
14: const int particles=100000;
16: // This function returns a random double between 0 and 1
17: double rnd() {return double(rand())/RAND_MAX;}
19: int main() {
20:         int i;
21:         double x,y,z,r;
22:         voronoicell c;
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;
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;
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:         }
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)) {
48:                 // Get the position of the current particle under consideration
49:                 clo.pos(x,y,z);
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);
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 {
68:                 // Get the position of the current particle under consideration
69:                 cls.pos(x,y,z);
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: }