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