branches/exact/examples/interface/polygons.cc

   1 // Direct C++ interface 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 void draw_polygon(FILE *fp,vector<int> &f_vert,vector<double> &v,int j);
  11
  12 int main() {
  13         unsigned int i,j;
  14         int id,nx,ny,nz;
  15         double x,y,z;
  16         voronoicell_neighbor c;
  17         vector<int> neigh,f_vert;
  18         vector<double> v;
  19
  20         // Create a pre-container class to import the input file and guess the
  21         // best computational grid size to use.
  22         pre_container pcon(-3,3,-3,3,0,6,false,false,false);
  23         pcon.import("pack_six_cube");
  24         pcon.guess_optimal(nx,ny,nz);
  25
  26         // Set up the container class and import the particles from the
  27         // pre-container
  28         container con(-3,3,-3,3,0,6,nx,ny,nz,false,false,false,8);
  29         pcon.setup(con);
  30
  31         // Open the output files
  32         FILE *fp4=safe_fopen("polygons4_v.pov","w"),
  33              *fp5=safe_fopen("polygons5_v.pov","w"),
  34              *fp6=safe_fopen("polygons6_v.pov","w");
  35
  36         // Loop over all particles in the container and compute each Voronoi
  37         // cell
  38         c_loop_all cl(con);
  39         if(cl.start()) do if(con.compute_cell(c,cl)) {
  40                 cl.pos(x,y,z);id=cl.pid();
  41
  42                 // Gather information about the computed Voronoi cell
  43                 c.neighbors(neigh);
  44                 c.face_vertices(f_vert);
  45                 c.vertices(x,y,z,v);
  46
  47                 // Loop over all faces of the Voronoi cell
  48                 for(i=0,j=0;i<neigh.size();i++) {
  49
  50                         // Draw all quadrilaterals, pentagons, and hexagons.
  51                         // Skip if the neighbor information is smaller than
  52                         // this particle's ID, to avoid double counting. This
  53                         // also removes faces that touch the walls, since the
  54                         // neighbor information is set to negative numbers for
  55                         // these cases.
  56                         if(neigh[i]>id) {
  57                                 switch(f_vert[j]) {
  58                                         case 4: draw_polygon(fp4,f_vert,v,j);
  59                                                 break;
  60                                         case 5: draw_polygon(fp5,f_vert,v,j);
  61                                                 break;
  62                                         case 6: draw_polygon(fp6,f_vert,v,j);
  63                                 }
  64                         }
  65
  66                         // Skip to the next entry in the face vertex list
  67                         j+=f_vert[j]+1;
  68                 }
  69         } while (cl.inc());
  70
  71         // Close the output files
  72         fclose(fp4);
  73         fclose(fp5);
  74         fclose(fp6);
  75
  76         // Draw the outline of the domain
  77         con.draw_domain_pov("polygons_d.pov");
  78 }
  79
  80 void draw_polygon(FILE *fp,vector<int> &f_vert,vector<double> &v,int j) {
  81         static char s[6][128];
  82         int k,l,n=f_vert[j];
  83
  84         // Create POV-Ray vector strings for each of the vertices
  85         for(k=0;k<n;k++) {
  86                 l=3*f_vert[j+k+1];
  87                 sprintf(s[k],"<%g,%g,%g>",v[l],v[l+1],v[l+2]);
  88         }
  89
  90         // Draw the interior of the polygon
  91         fputs("union{\n",fp);
  92         for(k=2;k<n;k++) fprintf(fp,"\ttriangle{%s,%s,%s}\n",s[0],s[k-1],s[k]);
  93         fputs("\ttexture{t1}\n}\n",fp);
  94
  95         // Draw the outline of the polygon
  96         fputs("union{\n",fp);
  97         for(k=0;k<n;k++) {
  98                 l=(k+1)%n;
  99                 fprintf(fp,"\tcylinder{%s,%s,r}\n\tsphere{%s,r}\n",
 100                         s[k],s[l],s[l]);
 101         }
 102         fputs("\ttexture{t2}\n}\n",fp);
 103 }
 104