branches/exact/examples/no_release/finite_sys.cc

   1 // Irregular packing 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=-15,x_max=15;
  12 const double y_min=-7,y_max=7;
  13 const double z_min=-15,z_max=15;
  14
  15 // Golden ratio constants
  16 const double Phi=0.5*(1+sqrt(5.0));
  17 const double phi=0.5*(1-sqrt(5.0));
  18
  19 // Set up the number of blocks that the container is divided
  20 // into.
  21 const int n_x=8,n_y=8,n_z=8;
  22
  23 // ID for dodecahedron faces
  24 const int wid=-10;
  25
  26 // Create a wall class that, whenever called, will replace the Voronoi cell
  27 // with a prescribed shape, in this case a dodecahedron
  28 class wall_initial_shape : public wall {
  29         public:
  30                 wall_initial_shape() {
  31
  32                         // Create a dodecahedron with neighbor information all
  33                         // set to -10
  34                         v.init(-2,2,-2,2,-2,2);
  35                         v.nplane(0,Phi,1,wid);v.nplane(0,-Phi,1,wid);v.nplane(0,Phi,-1,wid);
  36                         v.nplane(0,-Phi,-1,wid);v.nplane(1,0,Phi,wid);v.nplane(-1,0,Phi,wid);
  37                         v.nplane(1,0,-Phi,wid);v.nplane(-1,0,-Phi,wid);v.nplane(Phi,1,0,wid);
  38                         v.nplane(-Phi,1,0,wid);v.nplane(Phi,-1,0,wid);v.nplane(-Phi,-1,0,wid);
  39                 };
  40                 bool point_inside(double x,double y,double z) {return true;}
  41                 bool cut_cell(voronoicell &c,double x,double y,double z) {
  42
  43                         // Just ignore this case
  44                         return true;
  45                 }
  46                 bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) {
  47
  48                         // Set the cell to be equal to the dodecahedron
  49                         c=v;
  50                         return true;
  51                 }
  52         private:
  53                 voronoicell_neighbor v;
  54 };
  55
  56 // Determines whether any of the sides in the neighbor information are from the
  57 // initial dodecahedron
  58 bool has_dodec_sides(vector<int> &vi) {
  59         for(unsigned int i=0;i<vi.size();i++) if(vi[i]==wid) return true;
  60         return false;
  61 }
  62
  63 int main() {
  64
  65         // Create a container with the geometry given above. This is bigger
  66         // than the particle packing itself.
  67         container con(x_min,x_max,y_min,y_max,z_min,z_max,n_x,n_y,n_z,
  68                         false,false,false,8);
  69
  70         // Create the "initial shape" wall class and add it to the container
  71         wall_initial_shape(wis);
  72         con.add_wall(wis);
  73
  74         // Import the irregular particle packing
  75         con.import("pack_semicircle");
  76
  77         // Open files to save the "inside" and "outside" particles
  78         FILE *finside=safe_fopen("finite_sys_in.pov","w"),
  79              *foutside=safe_fopen("finite_sys_out.pov","w");
  80
  81         // Loop over all particles
  82         double x,y,z;
  83         vector<int> vi;
  84         voronoicell_neighbor c;
  85         c_loop_all cl(con);
  86         if(cl.start()) do {
  87
  88                 // Get particle position
  89                 cl.pos(x,y,z);
  90
  91                 // Remove half of the particles to see a cross-section
  92                 if(y<0) continue;
  93
  94                 if(con.compute_cell(c,cl)) {
  95
  96                         // Get the neighboring IDs of all the faces
  97                         c.neighbors(vi);
  98
  99                         // Depending on whether any of the faces are
 100                         // from the original dodecahedron, print to
 101                         // the "inside" or "outside" file
 102                         fprintf(has_dodec_sides(vi)?foutside:finside,
 103                                 "sphere{<%g,%g,%g>,s}\n",x,y,z);
 104                 }
 105         } while(cl.inc());
 106
 107         // Close the output files
 108         fclose(finside);
 109         fclose(foutside);
 110 }