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