Voro++

degenerate2.cc – A complicated degenerate vertex example

This is a further example of creating a Voronoi cell with degenerate high order vertices. It uses a procedure similar to the previous degenerate example, applying many rotated planes around a particular axis in order to create a high order vertex. However, unlike the previous example, the axes are chosen at random rather than along the coordinate directions, resulting in a complicated edge topology.

The code outputs the results in Gnuplot and POV-Ray formats, and the above image was generated using the POV-Ray header file “degenerate2.pov”.

Code listing

 1: // Degenerate vertex 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: // The total number of points to create as degenerate vertices
13: const int points=100;
14: 
15: // The number of planes that will be cut around each point
16: const int n=64;
17: const double step=2*pi/n;
18: 
19: // The angle (in radians) of the cutting planes from horizontal
20: const double theta=0.04;
21: 
22: // This function returns a random floating point number between 0 and 1
23: double rnd() {return double(rand())/RAND_MAX;}
24: 
25: int main() {
26:         double x,y,z,rsq,r,phi;
27:         voronoicell v;
28:         int n=0;
29: 
30:         // Initialize the Voronoi cell to be a cube of side length 2, centered on
31:         // the origin
32:         v.init(-1,1,-1,1,-1,1);
33: 
34:         // Plane cutting
35:         while(n<points) {
36: 
37:                 // Choose a random point
38:                 x=2*rnd()-1;
39:                 y=2*rnd()-1;
40:                 z=2*rnd()-1;
41: 
42:                 // Skip it if it's outside the unit sphere or too close to the
43:                 // origin
44:                 rsq=x*x+y*y+z*z;
45:                 if(rsq<0.01||rsq>1) continue;
46: 
47:                 // Rescale the point so that it has modulus 1, and then apply
48:                 // plane cuts around this point
49:                 r=1/sqrt(rsq);x*=r;y*=r;z*=r;
50:                 rsq=sqrt(x*x+y*y);r=z/rsq;
51:                 for(phi=rnd()*step;phi<2*pi;phi+=step)
52:                         v.plane(x*cos(theta)+sin(theta)*(-y*cos(phi)/rsq-x*r*sin(phi)),
53:                                 y*cos(theta)+sin(theta)*(x*cos(phi)/rsq-y*r*sin(phi)),
54:                                 z*cos(theta)+sin(theta)*rsq*sin(phi),1);
55:                 n++;
56:         }
57: 
58:         // Output the Voronoi cell to a file in Gnuplot format
59:         v.draw_gnuplot(0,0,0,"degenerate2.gnu");
60: 
61:         // Optional POV-Ray output
62:         v.draw_pov(0,0,0,"degenerate2_v.pov");
63:         v.draw_pov_mesh(0,0,0,"degenerate2_m.pov");
64: }