branches/2d/boundary/voro++_bdry.cc

   1 // Voro++, a cell-based Voronoi library
   2 //
   3 // Authors  : Chris H. Rycroft (LBL / UC Berkeley)
   4 //            Cody Robert Dance (UC Berkeley)
   5 // Email    : chr@alum.mit.edu
   6 // Date     : August 30th 2011
   7
   8 #include "voro++_2d.hh"
   9 using namespace voro;
  10
  11 #include <vector>
  12 #include <cstring>
  13 using namespace std;
  14
  15 const double pad=0.01;
  16
  17 int main(int argc,char **argv) {
  18
  19         if(argc<2) {
  20                 fprintf(stderr,"Syntax: voro++_bdry <input_file>\n");
  21                 return 1;
  22         }
  23
  24         // Read file to determine the size of the system and the number
  25         // of particles
  26         unsigned int i=0,j;
  27         int id;
  28         double x,y,minx=large_number,maxx=-minx,miny=minx,maxy=maxx;
  29         vector<int> vid;
  30         vector<double> vpos;
  31         vector<char> vbd;
  32         int mid=0;
  33         bool neg_label=false,boundary_track=false,start=false;
  34         char *buf(new char[512]);
  35
  36         FILE *fp=safe_fopen(argv[1],"r");
  37
  38         while(fgets(buf,512,fp)!=NULL) {
  39
  40                 if(strcmp(buf,"#Start\n")==0||strcmp(buf,"# Start\n")==0) {
  41
  42                         // Check that two consecutive start tokens haven't been
  43                         // encountered
  44                         if(boundary_track) voro_fatal_error("File import error - two consecutive start tokens found",VOROPP_FILE_ERROR);
  45                         start=true;boundary_track=true;
  46
  47                 } else if(strcmp(buf,"#End\n")==0||strcmp(buf,"# End\n")==0||
  48                           strcmp(buf,"#End")==0||strcmp(buf,"# End")==0) {
  49
  50                         // Check that two consecutive end tokens haven't been
  51                         // encountered
  52                         if(start) voro_fatal_error("File import error - end token immediately after start token",VOROPP_FILE_ERROR);
  53                         if(!boundary_track) voro_fatal_error("File import error - found end token without start token",VOROPP_FILE_ERROR);
  54                         vbd[i-1]|=2;boundary_track=false;
  55                 } else {
  56
  57                         // Try and read three entries from the line
  58                         if(sscanf(buf,"%d %lg %lg",&id,&x,&y)!=3) voro_fatal_error("File import error #1",VOROPP_FILE_ERROR);
  59                         vid.push_back(id);
  60                         vpos.push_back(x);
  61                         vpos.push_back(y);
  62                         vbd.push_back(start?1:0);
  63                         i++;
  64
  65                         // Determine bounds
  66                         if(id<0) neg_label=true;
  67                         if(id>mid) mid=id;
  68                         if(x<minx) minx=x;
  69                         if(x>maxx) maxx=x;
  70                         if(y<miny) miny=y;
  71                         if(y>maxy) maxy=y;
  72
  73                         start=false;
  74                 }
  75         }
  76
  77         if(boundary_track) voro_fatal_error("File import error - boundary not finished",VOROPP_FILE_ERROR);
  78         if(!feof(fp)) voro_fatal_error("File import error #2",VOROPP_FILE_ERROR);
  79         delete [] buf;
  80
  81         // Add small amount of padding to container bounds
  82         double dx=maxx-minx,dy=maxy-miny;
  83         minx-=pad*dx;maxx+=pad*dx;dx+=2*pad*dx;
  84         miny-=pad*dy;maxy+=pad*dy;dy+=2*pad*dy;
  85
  86         // Guess the optimal computationl grid, aiming at eight particles per
  87         // grid square
  88         double lscale=sqrt(8.0*dx*dy/i);
  89         int nx=int(dx/lscale)+1,ny=int(dy/lscale)+1;
  90
  91         // Print diagnostic information
  92         printf("Container bounds : [%g:%g] [%g:%g]\n"
  93                "Total particles  : %d\n"
  94                "Compute grid     : %d by %d\n",minx,maxx,miny,maxy,i,nx,ny);
  95
  96         // Create container
  97         container_boundary_2d con(minx,maxx,miny,maxy,nx,ny,false,false,16);
  98
  99         // Import data
 100         for(j=0;j<vid.size();j++) {
 101                 if(vbd[j]&1) con.start_boundary();
 102                 con.put(vid[j],vpos[2*j],vpos[2*j+1]);
 103                 if(vbd[j]&2) con.end_boundary();
 104         }
 105
 106         // Carry out all of the setup prior to computing any Voronoi cells
 107         con.setup();
 108
 109         // Save the boundary in a format that can be read by Gnuplot
 110         char *outfn(new char[strlen(argv[1])+5]);
 111         sprintf(outfn,"%s.bd",argv[1]);
 112
 113         con.draw_boundary_gnuplot(outfn);
 114
 115         //Test Global Connectivity Info
 116 //      cout << "beginning" << endl;
 117 //      char *connectfn(new char[strlen(argv[1])+5]);
 118 //      sprintf(connectfn,"%s.connect",argv[1]);
 119 //      con.print_custom(argv[2],connectfn );
 120 //      cout << "finished" << endl;
 121
 122         // Compute the Voronoi cells and save them to file
 123         sprintf(outfn,"%s.gnu",argv[1]);
 124         con.draw_cells_gnuplot(outfn);
 125
 126         // Output the neighbor mesh in gnuplot format
 127         if(mid>16777216) puts("Network output disabled due to too high Ids");
 128         else if(neg_label) puts("Network output disabled due to negative IDs");
 129         else {
 130
 131                 int *mp=new int[mid+1];
 132                 for(j=0;j<i;j++) mp[vid[j]]=j;
 133
 134                 sprintf(outfn,"%s.net",argv[1]);
 135                 FILE *ff=safe_fopen(outfn,"w");
 136                 int l1,l2;
 137                 vector<int> vi;
 138                 voronoicell_nonconvex_neighbor_2d c;
 139                 c_loop_all_2d cl(con);
 140                 if(cl.start()) do if(con.compute_cell(c,cl)) {
 141                         id=cl.pid();l1=2*mp[id];
 142                         c.neighbors(vi);
 143                         for(j=0;j<vi.size();j++) if(vi[j]>id) {
 144                                 l2=2*mp[vi[j]];
 145                                 fprintf(ff,"%g %g\n%g %g\n\n\n",
 146                                         vpos[l1],vpos[l1+1],vpos[l2],vpos[l2+1]);
 147                         }
 148                 } while (cl.inc());
 149                 fclose(ff);
 150
 151                 delete [] mp;
 152         }
 153 //      delete [] connectfn;
 154         delete [] outfn;
 155         return 0;
 156 }