Voro++

polygons.cc – Detecting polygonal faces with a certain number of sides

Voronoi cells are used frequently to analyze particle packings, and this code demonstrates how the C++ interface can be used to examine information about specific faces within a Voronoi tessellation, as one example of this type of analysis. The code computes a Voronoi tessellation and then extracts all of the faces with four, five, and six sides, to create the images above.

To begin, a pre_container class is created to be a cube with side length six. Usually, when a container class is used, it is necessary to specify how many computational blocks the container is subdivided into in each of the three directions, with an aim being to end up with around 3–12 particles in each block. In many cases, the user knows beforehand how many particles to expect, and they can choose the number of blocks accordingly. However, in situations where the number of particles is not known a priori, the pre_container class can be used. The class can import particles from a file, and dynamically allocates memory as needed to read in the entire file (line 26). Once the total number of particles is known, the guess_optimal routine (line 27) chooses an appropriate size of grid. A container class can then be set up (line 31), and the particles can be transferred to it (line 32).

The code then creates a c_loop_all class to loop over all of the particles in the container. For each particle, it computes the Voronoi cell and then gathers information about the computed cell. On line 46, a vector of the neighboring particle IDs corresponding to each face is returned. On line 47, a call is made to the face_vertices routine, which returns information about which vertices comprise each face. It is a vector of integers with a specific format: the first entry is a number k corresponding to the number of vertices making up a face, and this is followed k additional entries describing which vertices make up this face. For example, the sequence (3, 16, 20, 13) would correspond to a triangular face linking vertices 16, 20, and 13 together. On line 48, the vertex positions are returned: this corresponds to a vector of triplets (x, y, z) describing the position of each vertex.

On line 51, the code loops over all of the faces of the Voronoi cell and determines if the need to be printed. To avoid double counting, it only considers faces whose neighbor ID is larger than the particle ID itself. If the face has four, five, or six vertices, it calls a subroutine to draw that face to a specific file. Finally, on line 80, the outline of the container domain is also outputted. The POV-Ray files can be rendered using

which produces the images shown above.

Code listing

  1: // Direct C++ interface 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 <vector>
  8: using namespace std;
  9: 
 10: #include "voro++.hh"
 11: using namespace voro;
 12: 
 13: void draw_polygon(FILE *fp,vector<int> &f_vert,vector<double> &v,int j);
 14: 
 15: int main() {
 16:         unsigned int i,j;
 17:         int id,nx,ny,nz;
 18:         double x,y,z;
 19:         voronoicell_neighbor c;
 20:         vector<int> neigh,f_vert;
 21:         vector<double> v;
 22: 
 23:         // Create a pre-container class to import the input file and guess the
 24:         // best computational grid size to use.
 25:         pre_container pcon(-3,3,-3,3,0,6,false,false,false);
 26:         pcon.import("pack_six_cube");
 27:         pcon.guess_optimal(nx,ny,nz);
 28: 
 29:         // Set up the container class and import the particles from the
 30:         // pre-container
 31:         container con(-3,3,-3,3,0,6,nx,ny,nz,false,false,false,8);
 32:         pcon.setup(con);
 33: 
 34:         // Open the output files
 35:         FILE *fp4=safe_fopen("polygons4_v.pov","w"),
 36:              *fp5=safe_fopen("polygons5_v.pov","w"),
 37:              *fp6=safe_fopen("polygons6_v.pov","w");
 38: 
 39:         // Loop over all particles in the container and compute each Voronoi
 40:         // cell
 41:         c_loop_all cl(con);
 42:         if(cl.start()) do if(con.compute_cell(c,cl)) {
 43:                 cl.pos(x,y,z);id=cl.pid();
 44: 
 45:                 // Gather information about the computed Voronoi cell
 46:                 c.neighbors(neigh);
 47:                 c.face_vertices(f_vert);
 48:                 c.vertices(x,y,z,v);
 49: 
 50:                 // Loop over all faces of the Voronoi cell
 51:                 for(i=0,j=0;i<neigh.size();i++) {
 52: 
 53:                         // Draw all quadrilaterals, pentagons, and hexagons.
 54:                         // Skip if the neighbor information is smaller than
 55:                         // this particle's ID, to avoid double counting. This
 56:                         // also removes faces that touch the walls, since the
 57:                         // neighbor information is set to negative numbers for
 58:                         // these cases.
 59:                         if(neigh[i]>id) {
 60:                                 switch(f_vert[j]) {
 61:                                         case 4: draw_polygon(fp4,f_vert,v,j);
 62:                                                 break;
 63:                                         case 5: draw_polygon(fp5,f_vert,v,j);
 64:                                                 break;
 65:                                         case 6: draw_polygon(fp6,f_vert,v,j);
 66:                                 }
 67:                         }
 68: 
 69:                         // Skip to the next entry in the face vertex list
 70:                         j+=f_vert[j]+1;
 71:                 }
 72:         } while (cl.inc());
 73: 
 74:         // Close the output files
 75:         fclose(fp4);
 76:         fclose(fp5);
 77:         fclose(fp6);
 78: 
 79:         // Draw the outline of the domain
 80:         con.draw_domain_pov("polygons_d.pov");
 81: }
 82: 
 83: void draw_polygon(FILE *fp,vector<int> &f_vert,vector<double> &v,int j) {
 84:         static char s[6][128];
 85:         int k,l,n=f_vert[j];
 86: 
 87:         // Create POV-Ray vector strings for each of the vertices
 88:         for(k=0;k<n;k++) {
 89:                 l=3*f_vert[j+k+1];
 90:                 sprintf(s[k],"<%g,%g,%g>",v[l],v[l+1],v[l+2]);
 91:         }
 92: 
 93:         // Draw the interior of the polygon
 94:         fputs("union{\n",fp);
 95:         for(k=2;k<n;k++) fprintf(fp,"\ttriangle{%s,%s,%s}\n",s[0],s[k-1],s[k]);
 96:         fputs("\ttexture{t1}\n}\n",fp);
 97: 
 98:         // Draw the outline of the polygon
 99:         fputs("union{\n",fp);
100:         for(k=0;k<n;k++) {
101:                 l=(k+1)%n;
102:                 fprintf(fp,"\tcylinder{%s,%s,r}\n\tsphere{%s,r}\n",
103:                         s[k],s[l],s[l]);
104:         }
105:         fputs("\ttexture{t2}\n}\n",fp);
106: }
107: