trunk/examples/walls/frustum.cc

   1 // Frustum 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 double rnd() {return 0.001*double(rand())/RAND_MAX;}
  13
  14 int main() {
  15         int i=0;
  16         double x,y,z,evol,vvol;
  17
  18         // Create a container with the geometry given above, and make it
  19         // non-periodic in each of the three coordinates. Allocate space for
  20         // eight particles within each computational block.
  21         container con(-1.2,1.2,-1.2,1.2,0,1,14,14,7,
  22                         false,false,false,8);
  23
  24         // Add a cylindrical wall to the container
  25         wall_cone cone(0,0,2,0,0,-1,atan(0.5));
  26         con.add_wall(cone);
  27
  28         // Place particles in a regular grid within the frustum, for points
  29         // which are within the wall boundaries
  30         for(z=0.1;z<1;z+=0.2) for(y=-0.85;y<1;y+=0.2) for(x=-0.95;x<1;x+=0.2) {
  31                 if (con.point_inside(x,y,z)) {
  32                         con.put(i,x,y,z);i++;
  33                 }
  34         }
  35
  36         // Output the particle positions and Voronoi cells in Gnuplot format
  37         con.draw_particles("frustum_p.gnu");
  38         con.draw_cells_gnuplot("frustum_v.gnu");
  39
  40         // Output the particle positions and Voronoi cells in POV-Ray format
  41         con.draw_particles_pov("frustum_p.pov");
  42         con.draw_cells_pov("frustum_v.pov");
  43
  44         // Compute the volume of the Voronoi cells and compare it to the
  45         // exact frustum volume
  46         evol=pi*1*(0.5*0.5+0.5*1+1*1)/3;
  47         vvol=con.sum_cell_volumes();
  48         printf("Exact frustum volume : %g\n"
  49                "Voronoi cell volume  : %g\n"
  50                "Difference           : %g\n",evol,vvol,vvol-evol);
  51 }