cut_region.cc – The region 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 : [email protected] 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: }