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
povray +H600 +W600 +A0.1 -J +Opolygons4.png polygons4.povpovray +H600 +W600 +A0.1 -J +Opolygons5.png polygons5.pov
povray +H600 +W600 +A0.1 -J +Opolygons6.png polygons6.pov
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: