trunk/examples/walls/torus.cc

   1 // Custom wall class 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 // Major and minor torus radii
  11 const double arad=9,brad=3.5;
  12
  13 // The outer radius of the torus, that determines how big the container should
  14 // be
  15 const double crad=arad+brad;
  16
  17 // Set up constants for the container geometry
  18 const double x_min=-crad-0.5,x_max=crad+0.5;
  19 const double y_min=-crad-0.5,y_max=crad+0.5;
  20 const double z_min=-brad-0.5,z_max=brad+0.5;
  21
  22 // Set the computational grid size
  23 const int n_x=10,n_y=10,n_z=3;
  24
  25 // This class creates a custom toroidal wall object that is centered on the
  26 // origin and is aligned with the xy plane. It is derived from the pure virtual
  27 // "wall" class. The "wall" class contains virtual functions for cutting the
  28 // Voronoi cell in response to a wall, and for telling whether a given point is
  29 // inside the wall or not. In this derived class, specific implementations of
  30 // these functions are given.
  31 class wall_torus : public wall {
  32         public:
  33
  34                 // The wall constructor initializes constants for the major and
  35                 // minor axes of the torus. It also initializes the wall ID
  36                 // number that is used when the plane cuts are made. This is
  37                 // only tracked with the voronoicell_neighbor class and is
  38                 // ignored otherwise. It can be omitted, and then an arbitrary
  39                 // value of -99 is used.
  40                 wall_torus(double imjr,double imnr,int iw_id=-99)
  41                         : w_id(iw_id), mjr(imjr), mnr(imnr) {};
  42
  43                 // This returns true if a given vector is inside the torus, and
  44                 // false if it is outside. For the current example, this
  45                 // routine is not needed, but in general it would be, for use
  46                 // with the point_inside() routine in the container class.
  47                 bool point_inside(double x,double y,double z) {
  48                         double temp=sqrt(x*x+y*y)-mjr;
  49                         return temp*temp+z*z<mnr*mnr;
  50                 }
  51
  52                 // This template takes a reference to a voronoicell or
  53                 // voronoicell_neighbor object for a particle at a vector
  54                 // (x,y,z), and makes a plane cut to to the object to account
  55                 // for the toroidal wall
  56                 template<class vc_class>
  57                 inline bool cut_cell_base(vc_class &c,double x,double y,double z) {
  58                         double orad=sqrt(x*x+y*y);
  59                         double odis=orad-mjr;
  60                         double ot=odis*odis+z*z;
  61
  62                         // Unless the particle is within 1% of the major
  63                         // radius, then a plane cut is made
  64                         if(ot>0.01*mnr) {
  65                                 ot=2*mnr/sqrt(ot)-2;
  66                                 z*=ot;
  67                                 odis*=ot/orad;
  68                                 x*=odis;
  69                                 y*=odis;
  70                                 return c.nplane(x,y,z,w_id);
  71                         }
  72                         return true;
  73                 }
  74
  75                 // These virtual functions are called during the cell
  76                 // computation in the container class. They call instances of
  77                 // the template given above.
  78                 bool cut_cell(voronoicell &c,double x,
  79                                 double y,double z) {return cut_cell_base(c,x,y,z);}
  80                 bool cut_cell(voronoicell_neighbor &c,double x,
  81                                 double y,double z) {return cut_cell_base(c,x,y,z);}
  82         private:
  83                 // The ID number associated with the wall
  84                 const int w_id;
  85                 // The major radius of the torus
  86                 const double mjr;
  87                 // The minor radius of the torus
  88                 const double mnr;
  89 };
  90
  91 int main() {
  92
  93         // Create a container with the geometry given above, and make it
  94         // non-periodic in each of the three coordinates. Allocate space for
  95         // eight particles within each computational block.
  96         container con(x_min,x_max,y_min,y_max,z_min,z_max,n_x,n_y,n_z,
  97                         false,false,false,8);
  98
  99         // Add the custom toroidal wall to the container
 100         wall_torus tor(arad,brad);
 101         con.add_wall(tor);
 102
 103         // Import the particles from a file
 104         con.import("pack_torus");
 105
 106         // Output the particle positions in POV-Ray format
 107         con.draw_particles_pov("torus_p.pov");
 108
 109         // Output the Voronoi cells in POV-Ray format
 110         con.draw_cells_pov("torus_v.pov");
 111 }