trunk/examples/no_release/cylinder_inv.cc

   1 // Cylindrical wall 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 // Set up constants for the container geometry
  11 const double x_min=-6,x_max=6;
  12 const double y_min=-6,y_max=6;
  13 const double z_min=-6,z_max=6;
  14
  15 // Set the computational grid size
  16 const int n_x=6,n_y=6,n_z=6;
  17
  18 struct wall_cylinder_inv : public wall {
  19         public:
  20                 wall_cylinder_inv(double xc_,double yc_,double zc_,double xa_,double ya_,double za_,double rc_,int w_id_=-99)
  21                         : w_id(w_id_), xc(xc_), yc(yc_), zc(zc_), xa(xa_), ya(ya_), za(za_),
  22                         asi(1/(xa_*xa_+ya_*ya_+za_*za_)), rc(rc_) {}
  23                 bool point_inside(double x,double y,double z) {
  24                         double xd=x-xc,yd=y-yc,zd=z-zc;
  25                         double pa=(xd*xa+yd*ya+zd*za)*asi;
  26                         xd-=xa*pa;yd-=ya*pa;zd-=za*pa;
  27                         return xd*xd+yd*yd+zd*zd>rc*rc;
  28                 }
  29                 template<class v_cell>
  30                 bool cut_cell_base(v_cell &c,double x,double y,double z) {
  31                         double xd=x-xc,yd=y-yc,zd=z-zc;
  32                         double pa=(xd*xa+yd*ya+zd*za)*asi;
  33                         xd-=xa*pa;yd-=ya*pa;zd-=za*pa;
  34                         pa=xd*xd+yd*yd+zd*zd;
  35                         if(pa>1e-5) {
  36                                 pa=2*(sqrt(pa)*rc-pa);
  37                                 return c.nplane(-xd,-yd,-zd,-pa,w_id);
  38                         }
  39                         return true;
  40                 }
  41                 bool cut_cell(voronoicell &c,double x,double y,double z) {return cut_cell_base(c,x,y,z);}
  42                 bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) {return cut_cell_base(c,x,y,z);}
  43         private:
  44                 const int w_id;
  45                 const double xc,yc,zc,xa,ya,za,asi,rc;
  46 };
  47
  48 int main() {
  49         int i;double x,y,z;
  50
  51         // Create a container with the geometry given above, and make it
  52         // non-periodic in each of the three coordinates. Allocate space for
  53         // eight particles within each computational block.
  54         container con(x_min,x_max,y_min,y_max,z_min,z_max,n_x,n_y,n_z,
  55                         false,false,false,8);
  56
  57         // Add a cylindrical wall to the container
  58         wall_cylinder_inv cyl(0,0,0,0,1,0,4.2);
  59         con.add_wall(cyl);
  60
  61         // Place particles in a regular grid within the frustum, for points
  62         // which are within the wall boundaries
  63         for(z=-5.5;z<6;z+=1) for(y=-5.5;y<6;y+=1) for(x=-5.5;x<6;x+=1)
  64                 if (con.point_inside(x,y,z)) {con.put(i,x,y,z);i++;}
  65
  66         // Output the particle positions in POV-Ray format
  67         con.draw_particles_pov("cylinder_inv_p.pov");
  68
  69         // Output the Voronoi cells in POV-Ray format
  70         con.draw_cells_pov("cylinder_inv_v.pov");
  71
  72         // Output the particle positions in gnuplot format
  73         con.draw_particles("cylinder_inv.par");
  74
  75         // Output the Voronoi cells in gnuplot format
  76         con.draw_cells_gnuplot("cylinder_inv.gnu");
  77 }