trunk/zeo/network.cc

   1 // Voronoi calculation 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 <cstring>
   8
   9 #include "voro++.hh"
  10 using namespace voro;
  11
  12 #include "v_network.hh"
  13 #include "r_table.cc"
  14
  15 // A guess for the memory allocation per region
  16 const int memory=16;
  17
  18 // A maximum allowed number of regions, to prevent enormous amounts of memory
  19 // being allocated
  20 const int max_regions=16777216;
  21
  22 // A buffer size
  23 const int bsize=2048;
  24
  25 // Output routine
  26 template<class c_class>
  27 void compute(c_class &con,char *buffer,int bp,double vol);
  28
  29 // Commonly used error message
  30 void file_import_error() {
  31         voro_fatal_error("File import error",VOROPP_FILE_ERROR);
  32 }
  33
  34 int main(int argc,char **argv) {
  35         char *farg,buffer[bsize];
  36         bool radial;int i,n,bp;
  37         double bx,bxy,by,bxz,byz,bz,x,y,z,vol;
  38
  39         // Check the command line syntax
  40         if(argc==2) {
  41                 radial=false;farg=argv[1];
  42         } else if(argc==3&&strcmp(argv[1],"-r")==0) {
  43                 radial=true;farg=argv[2];
  44         } else {
  45                 fputs("Syntax: ./network [-r] <filename.v1>\n",stderr);
  46                 return VOROPP_CMD_LINE_ERROR;
  47         }
  48
  49         // Check that the file has a ".v1" extension
  50         bp=strlen(farg);
  51         if(bp+2>bsize) {
  52                 fputs("Filename too long\n",stderr);
  53                 return VOROPP_CMD_LINE_ERROR;
  54         }
  55         if(bp<3||farg[bp-3]!='.'||farg[bp-2]!='v'||farg[bp-1]!='1') {
  56                 fputs("Filename must end in '.v1'\n",stderr);
  57                 return VOROPP_CMD_LINE_ERROR;
  58         }
  59
  60         // Try opening the file
  61         FILE *fp(fopen(farg,"r"));
  62         if(fp==NULL) voro_fatal_error("Unable to open file for import",VOROPP_FILE_ERROR);
  63
  64         // Read header line
  65         if(fgets(buffer,bsize,fp)!=buffer) file_import_error();
  66         if(strcmp(buffer,"Unit cell vectors:\n")!=0)
  67                 voro_fatal_error("Invalid header line",VOROPP_FILE_ERROR);
  68
  69         // Read in the box dimensions and the number of particles
  70         if(fscanf(fp,"%s %lg %lg %lg",buffer,&bx,&x,&x)!=4) file_import_error();
  71         if(strcmp(buffer,"va=")!=0) voro_fatal_error("Invalid first vector",VOROPP_FILE_ERROR);
  72         if(fscanf(fp,"%s %lg %lg %lg",buffer,&bxy,&by,&x)!=4) file_import_error();
  73         if(strcmp(buffer,"vb=")!=0) voro_fatal_error("Invalid second vector",VOROPP_FILE_ERROR);
  74         if(fscanf(fp,"%s %lg %lg %lg",buffer,&bxz,&byz,&bz)!=4) file_import_error();
  75         if(strcmp(buffer,"vc=")!=0) voro_fatal_error("Invalid third vector",VOROPP_FILE_ERROR);
  76         if(fscanf(fp,"%d",&n)!=1) file_import_error();
  77
  78         // Print the box dimensions
  79         printf("Box dimensions:\n"
  80                "  va=(%f 0 0)\n"
  81                "  vb=(%f %f 0)\n"
  82                "  vc=(%f %f %f)\n\n",bx,bxy,by,bxz,byz,bz);
  83
  84         // Check that the input parameters make sense
  85         if(n<1) voro_fatal_error("Invalid number of particles",VOROPP_FILE_ERROR);
  86         if(bx<tolerance||by<tolerance||bz<tolerance)
  87                 voro_fatal_error("Invalid box dimensions",VOROPP_FILE_ERROR);
  88
  89         // Compute the internal grid size, aiming to make
  90         // the grid blocks square with around 6 particles
  91         // in each
  92         double ls=1.8*pow(bx*by*bz,-1.0/3.0);
  93         double nxf=bx*ls+1.5;
  94         double nyf=by*ls+1.5;
  95         double nzf=bz*ls+1.5;
  96
  97         // Check the grid is not too huge, using floating point numbers to avoid
  98         // integer wrap-arounds
  99         if (nxf*nyf*nzf>max_regions) {
 100                 fprintf(stderr,"voro++: Number of computational blocks exceeds the maximum allowed of %d\n"
 101                         "Either increase the particle length scale, or recompile with an increased\nmaximum.\n",
 102                         max_regions);
 103                 return VOROPP_MEMORY_ERROR;
 104         }
 105
 106         // Now that we are confident that the number of regions is reasonable,
 107         // create integer versions of them
 108         int nx=int(nxf);
 109         int ny=int(nyf);
 110         int nz=int(nzf);
 111         printf("Total particles = %d\n\nInternal grid size = (%d %d %d)\n\n",n,nx,ny,nz);
 112
 113         vol=bx*by*bz;
 114         if(radial) {
 115
 116                 // Create a container with the geometry given above
 117                 container_periodic_poly con(bx,bxy,by,bxz,byz,bz,nx,ny,nz,memory);
 118
 119                 // Read in the particles from the file
 120                 for(i=0;i<n;i++) {
 121                         if(fscanf(fp,"%s %lg %lg %lg",buffer,&x,&y,&z)!=4) file_import_error();
 122                         con.put(i,x,y,z,radial_lookup(buffer));
 123                 }
 124                 fclose(fp);
 125
 126                 // Copy the output filename
 127                 for(i=0;i<bp-2;i++) buffer[i]=farg[i];
 128                 compute(con,buffer,bp,vol);
 129         } else {
 130
 131                 // Create a container with the geometry given above
 132                 container_periodic con(bx,bxy,by,bxz,byz,bz,nx,ny,nz,memory);
 133
 134                 // Read in the particles from the file
 135                 for(i=0;i<n;i++) {
 136                         if(fscanf(fp,"%s %lg %lg %lg",buffer,&x,&y,&z)!=4)
 137                                 voro_fatal_error("File import error",VOROPP_FILE_ERROR);
 138                         con.put(i,x,y,z);
 139                 }
 140                 fclose(fp);
 141
 142                 // Copy the output filename
 143                 for(i=0;i<bp-2;i++) buffer[i]=farg[i];
 144                 compute(con,buffer,bp,vol);
 145         }
 146 }
 147
 148 inline void extension(const char *ext,char *bu) {
 149         char *ep((char*) ext);
 150         while(*ep!=0) *(bu++)=*(ep++);*bu=*ep;
 151 }
 152
 153 template<class c_class>
 154 void compute(c_class &con,char *buffer,int bp,double vol) {
 155         char *bu(buffer+bp-2);
 156         int id;
 157         double vvol(0),x,y,z,r;
 158         voronoicell c(con);
 159         voronoi_network vn(con,1e-5),vn2(con,1e-5);
 160
 161         // Compute Voronoi cells and
 162         c_loop_all_periodic vl(con);
 163         if(vl.start()) do if(con.compute_cell(c,vl)) {
 164                 vvol+=c.volume();
 165                 vl.pos(id,x,y,z,r);
 166                 vn.add_to_network(c,id,x,y,z,r);
 167                 vn2.add_to_network_rectangular(c,id,x,y,z,r);
 168         } while(vl.inc());
 169
 170         // Carry out the volume check
 171         printf("Volume check:\n  Total domain volume  = %f\n"
 172                "  Total Voronoi volume = %f\n",vol,vvol);
 173
 174         // Print non-rectangular cell network
 175         extension("nd2",bu);vn.draw_network(buffer);
 176         extension("nt2",bu);vn.print_network(buffer);
 177
 178         // Print rectangular cell network
 179         extension("ntd",bu);vn2.draw_network(buffer);
 180         extension("net",bu);vn2.print_network(buffer);
 181
 182         // Output the particles and any constructed periodic images
 183         extension("par",bu);con.draw_particles(buffer);
 184
 185         // Output the Voronoi cells in gnuplot format
 186         extension("out",bu);con.draw_cells_gnuplot(buffer);
 187
 188         // Output the unit cell in gnuplot format
 189         extension("dom",bu);con.draw_domain_gnuplot(buffer);
 190 }