cut_region.cc – The region that can cut a Voronoi cell

The volume that can cut a Voronoi cell

This example code computes the region of space where a particle can be placed that will result in a cut of a particular Voronoi cell. This code makes use of the plane_intersects() routine, that returns true or false depending on whether a given plane intersects the cell.

The code works by finding the boundary of the region that results in a plane cut. For each normal vector n, the code finds the distance r such that any plane closer than r will result in a cut, while any plane that is further away will not. In the code, the distances are found by using a bisection algorithm, by making repeated calls to the plane_intersects() routine to find the value of r to within a tolerance specified by tolwidth.

The code can be applied to any Voronoi cell. In this code, a simple cell is constructed in lines 23 to 24, by making an octahedron and applying a single plane cut to it. The resulting surface has a spherical structure that is not convex. The Voronoi cell computation routines make use of the plane_intersects() routine, in order to determine what particles need to be tested when constructing each cell.

Code listing

 1: // Cell cutting region 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 "voro++.hh"
 8: using namespace voro;
 9: 
10: const double pi=3.1415926535897932384626433832795;
11: 
12: // This constant sets the tolerance in the bisection search algorithm
13: const double tolwidth=1e-7;
14: 
15: // This constant determines the density of points to test
16: const double theta_step=pi/200;
17: 
18: int main() {
19:         double x,y,z,r,rmin,rmax;
20:         double theta,phi,phi_step;
21:         voronoicell v;
22:         FILE *fp=safe_fopen("cell_cut_region.gnu","w");
23: 
24:         // Initialize the Voronoi cell to be an octahedron and make a single
25:         // plane cut to add some variation
26:         v.init_octahedron(1);
27:         v.plane(0.4,0.3,1,0.1);
28: 
29:         // Output the cell in gnuplot format
30:         v.draw_gnuplot(0,0,0,"cell.gnu");
31: 
32:         // Now test over direction vectors from the center of the sphere. For
33:         // each vector, carry out a search to find the maximum distance along
34:         // that vector that a plane will intersect with cell, and save it to
35:         // the output file.
36:         for(theta=theta_step*0.5;theta<pi;theta+=theta_step) {
37:                 phi_step=2*pi/(int(2*pi*sin(theta)/theta_step)+1);
38:                 for(phi=phi_step*0.5;phi<2*pi;phi+=phi_step) {
39: 
40:                         // Calculate a direction to look along
41:                         x=sin(theta)*cos(phi);
42:                         y=sin(theta)*sin(phi);
43:                         z=cos(theta);
44: 
45:                         // Now carry out a bisection search. Here, we initialize
46:                         // a minimum and a maximum guess for the distance
47:                         // along this vector. Keep multiplying rmax by two until
48:                         // the plane no longer makes a cut.
49:                         rmin=0;rmax=1;
50:                         while (v.plane_intersects(x,y,z,rmax)) rmax*=2;
51: 
52:                         // We now know that the distance is somewhere between
53:                         // rmin and rmax. Test the point halfway in-between
54:                         // these two. If it intersects, then move rmin to this
55:                         // point; otherwise, move rmax there. At each stage the
56:                         // bounding interval is divided by two. Exit when the
57:                         // width of the interval is smaller than the tolerance.
58:                         while (rmax-rmin>tolwidth) {
59:                                 r=(rmax+rmin)*0.5;
60:                                 if (v.plane_intersects(x,y,z,r)) rmin=r;
61:                                 else rmax=r;
62:                         }
63: 
64:                         // Output this point to file
65:                         r=(rmax+rmin)*0.5;
66:                         x*=r;y*=r;z*=r;
67:                         fprintf(fp,"%g %g %g\n",x,y,z);
68:                 }
69:         }
70: 
71:         fclose(fp);
72: }