1 // Voro++, a 3D cell-based Voronoi library
2 //
3 // Author : Chris H. Rycroft (LBL / UC Berkeley)
4 // Email : chr@alum.mit.edu
5 // Date : August 30th 2011
6
7 /** \file cmd_line.cc
8 * \brief Source code for the command-line utility. */
9
10 #include <cstring>
11 #include "voro++.hh"
12 using namespace voro;
13
14 enum blocks_mode {
15 none,
16 length_scale,
17 specified
18 };
19
20 // A maximum allowed number of regions, to prevent enormous amounts of memory
21 // being allocated
22 const int max_regions=16777216;
23
24 // This message gets displayed if the user requests the help flag
help_message()25 void help_message() {
26 puts("Voro++ version 0.4.4, by Chris H. Rycroft (UC Berkeley/LBL)\n\n"
27 "Syntax: voro++ [options] <x_min> <x_max> <y_min>\n"
28 " <y_max> <z_min> <z_max> <filename>\n\n"
29 "By default, the utility reads in the input file of particle IDs and positions,\n"
30 "computes the Voronoi cell for each, and then creates <filename.vol> with an\n"
31 "additional column containing the volume of each Voronoi cell.\n\n"
32 "Available options:\n"
33 " -c <str> : Specify a custom output string\n"
34 " -g : Turn on the gnuplot output to <filename.gnu>\n"
35 " -h/--help : Print this information\n"
36 " -hc : Print information about custom output\n"
37 " -l <len> : Manually specify a length scale to configure the internal\n"
38 " computational grid\n"
39 " -m <mem> : Manually choose the memory allocation per grid block\n"
40 " (default 8)\n"
41 " -n [3] : Manually specify the internal grid size\n"
42 " -o : Ensure that the output file has the same order as the input\n"
43 " file\n"
44 " -p : Make container periodic in all three directions\n"
45 " -px : Make container periodic in the x direction\n"
46 " -py : Make container periodic in the y direction\n"
47 " -pz : Make container periodic in the z direction\n"
48 " -r : Assume the input file has an extra coordinate for radii\n"
49 " -v : Verbose output\n"
50 " --version : Print version information\n"
51 " -wb [6] : Add six plane wall objects to make rectangular box containing\n"
52 " the space x1<x<x2, x3<y<x4, x5<z<x6\n"
53 " -wc [7] : Add a cylinder wall object, centered on (x1,x2,x3),\n"
54 " pointing in (x4,x5,x6), radius x7\n"
55 " -wo [7] : Add a conical wall object, apex at (x1,x2,x3), axis\n"
56 " along (x4,x5,x6), angle x7 in radians\n"
57 " -ws [4] : Add a sphere wall object, centered on (x1,x2,x3),\n"
58 " with radius x4\n"
59 " -wp [4] : Add a plane wall object, with normal (x1,x2,x3),\n"
60 " and displacement x4\n"
61 " -y : Save POV-Ray particles to <filename_p.pov> and POV-Ray Voronoi\n"
62 " cells to <filename_v.pov>\n"
63 " -yp : Save only POV-Ray particles to <filename_p.pov>\n"
64 " -yv : Save only POV-Ray Voronoi cells to <filename_v.pov>");
65 }
66
67 // This message gets displayed if the user requests information about doing
68 // custom output
custom_output_message()69 void custom_output_message() {
70 puts("The \"-c\" option allows a string to be specified that will customize the output\n"
71 "file to contain a variety of statistics about each computed Voronoi cell. The\n"
72 "string is similar to the standard C printf() function, made up of text with\n"
73 "additional control sequences that begin with percentage signs that are expanded\n"
74 "to different statistics. See http://math.lbl.gov/voro++/doc/custom.html for more\n"
75 "information.\n"
76 "\nParticle-related:\n"
77 " %i The particle ID number\n"
78 " %x The x coordinate of the particle\n"
79 " %y The y coordinate of the particle\n"
80 " %z The z coordinate of the particle\n"
81 " %q The position vector of the particle, short for \"%x %y %z\"\n"
82 " %r The radius of the particle (only printed if -p enabled)\n"
83 "\nVertex-related:\n"
84 " %w The number of vertices in the Voronoi cell\n"
85 " %p A list of the vertices of the Voronoi cell in the format (x,y,z),\n"
86 " relative to the particle center\n"
87 " %P A list of the vertices of the Voronoi cell in the format (x,y,z),\n"
88 " relative to the global coordinate system\n"
89 " %o A list of the orders of each vertex\n"
90 " %m The maximum radius squared of a vertex position, relative to the\n"
91 " particle center\n"
92 "\nEdge-related:\n"
93 " %g The number of edges of the Voronoi cell\n"
94 " %E The total edge distance\n"
95 " %e A list of perimeters of each face\n"
96 "\nFace-related:\n"
97 " %s The number of faces of the Voronoi cell\n"
98 " %F The total surface area of the Voronoi cell\n"
99 " %A A frequency table of the number of edges for each face\n"
100 " %a A list of the number of edges for each face\n"
101 " %f A list of areas of each face\n"
102 " %t A list of bracketed sequences of vertices that make up each face\n"
103 " %l A list of normal vectors for each face\n"
104 " %n A list of neighboring particle or wall IDs corresponding to each face\n"
105 "\nVolume-related:\n"
106 " %v The volume of the Voronoi cell\n"
107 " %c The centroid of the Voronoi cell, relative to the particle center\n"
108 " %C The centroid of the Voronoi cell, in the global coordinate system");
109 }
110
111 // Ths message is displayed if the user requests version information
version_message()112 void version_message() {
113 puts("Voro++ version 0.4.4 (January 17th 2012)");
114 }
115
116 // Prints an error message. This is called when the program is unable to make
117 // sense of the command-line options.
error_message()118 void error_message() {
119 fputs("voro++: Unrecognized command-line options; type \"voro++ -h\" for more\ninformation.\n",stderr);
120 }
121
122 // Carries out the Voronoi computation and outputs the results to the requested
123 // files
124 template<class c_loop,class c_class>
cmd_line_output(c_loop & vl,c_class & con,const char * format,FILE * outfile,FILE * gnu_file,FILE * povp_file,FILE * povv_file,bool verbose,double & vol,int & vcc,int & tp)125 void cmd_line_output(c_loop &vl,c_class &con,const char* format,FILE* outfile,FILE* gnu_file,FILE* povp_file,FILE* povv_file,bool verbose,double &vol,int &vcc,int &tp) {
126 int pid,ps=con.ps;double x,y,z,r;
127 if(con.contains_neighbor(format)) {
128 voronoicell_neighbor c;
129 if(vl.start()) do if(con.compute_cell(c,vl)) {
130 vl.pos(pid,x,y,z,r);
131 if(outfile!=NULL) c.output_custom(format,pid,x,y,z,r,outfile);
132 if(gnu_file!=NULL) c.draw_gnuplot(x,y,z,gnu_file);
133 if(povp_file!=NULL) {
134 fprintf(povp_file,"// id %d\n",pid);
135 if(ps==4) fprintf(povp_file,"sphere{<%g,%g,%g>,%g}\n",x,y,z,r);
136 else fprintf(povp_file,"sphere{<%g,%g,%g>,s}\n",x,y,z);
137 }
138 if(povv_file!=NULL) {
139 fprintf(povv_file,"// cell %d\n",pid);
140 c.draw_pov(x,y,z,povv_file);
141 }
142 if(verbose) {vol+=c.volume();vcc++;}
143 } while(vl.inc());
144 } else {
145 voronoicell c;
146 if(vl.start()) do if(con.compute_cell(c,vl)) {
147 vl.pos(pid,x,y,z,r);
148 if(outfile!=NULL) c.output_custom(format,pid,x,y,z,r,outfile);
149 if(gnu_file!=NULL) c.draw_gnuplot(x,y,z,gnu_file);
150 if(povp_file!=NULL) {
151 fprintf(povp_file,"// id %d\n",pid);
152 if(ps==4) fprintf(povp_file,"sphere{<%g,%g,%g>,%g}\n",x,y,z,r);
153 else fprintf(povp_file,"sphere{<%g,%g,%g>,s}\n",x,y,z);
154 }
155 if(povv_file!=NULL) {
156 fprintf(povv_file,"// cell %d\n",pid);
157 c.draw_pov(x,y,z,povv_file);
158 }
159 if(verbose) {vol+=c.volume();vcc++;}
160 } while(vl.inc());
161 }
162 if(verbose) tp=con.total_particles();
163 }
164
main(int argc,char ** argv)165 int main(int argc,char **argv) {
166 int i=1,j=-7,custom_output=0,nx,ny,nz,init_mem(8);
167 double ls=0;
168 blocks_mode bm=none;
169 bool gnuplot_output=false,povp_output=false,povv_output=false,polydisperse=false;
170 bool xperiodic=false,yperiodic=false,zperiodic=false,ordered=false,verbose=false;
171 pre_container *pcon=NULL;pre_container_poly *pconp=NULL;
172 wall_list wl;
173
174 // If there's one argument, check to see if it's requesting help.
175 // Otherwise, bail out with an error.
176 if(argc==2) {
177 if(strcmp(argv[1],"-h")==0||strcmp(argv[1],"--help")==0) {
178 help_message();return 0;
179 } else if(strcmp(argv[1],"-hc")==0) {
180 custom_output_message();return 0;
181 } else if(strcmp(argv[1],"--version")==0) {
182 version_message();return 0;
183 } else {
184 error_message();
185 return VOROPP_CMD_LINE_ERROR;
186 }
187 }
188
189 // If there aren't enough command-line arguments, then bail out
190 // with an error.
191 if(argc<7) {
192 error_message();
193 return VOROPP_CMD_LINE_ERROR;
194 }
195
196 // We have enough arguments. Now start searching for command-line
197 // options.
198 while(i<argc-7) {
199 if(strcmp(argv[i],"-c")==0) {
200 if(i>=argc-8) {error_message();wl.deallocate();return VOROPP_CMD_LINE_ERROR;}
201 if(custom_output==0) {
202 custom_output=++i;
203 } else {
204 fputs("voro++: multiple custom output strings detected\n",stderr);
205 wl.deallocate();
206 return VOROPP_CMD_LINE_ERROR;
207 }
208 } else if(strcmp(argv[i],"-g")==0) {
209 gnuplot_output=true;
210 } else if(strcmp(argv[i],"-h")==0||strcmp(argv[i],"--help")==0) {
211 help_message();wl.deallocate();return 0;
212 } else if(strcmp(argv[i],"-hc")==0) {
213 custom_output_message();wl.deallocate();return 0;
214 } else if(strcmp(argv[i],"-l")==0) {
215 if(i>=argc-8) {error_message();wl.deallocate();return VOROPP_CMD_LINE_ERROR;}
216 if(bm!=none) {
217 fputs("voro++: Conflicting options about grid setup (-l/-n)\n",stderr);
218 wl.deallocate();
219 return VOROPP_CMD_LINE_ERROR;
220 }
221 bm=length_scale;
222 i++;ls=atof(argv[i]);
223 } else if(strcmp(argv[i],"-m")==0) {
224 i++;init_mem=atoi(argv[i]);
225 } else if(strcmp(argv[i],"-n")==0) {
226 if(i>=argc-10) {error_message();wl.deallocate();return VOROPP_CMD_LINE_ERROR;}
227 if(bm!=none) {
228 fputs("voro++: Conflicting options about grid setup (-l/-n)\n",stderr);
229 wl.deallocate();
230 return VOROPP_CMD_LINE_ERROR;
231 }
232 bm=specified;
233 i++;
234 nx=atoi(argv[i++]);
235 ny=atoi(argv[i++]);
236 nz=atoi(argv[i]);
237 if(nx<=0||ny<=0||nz<=0) {
238 fputs("voro++: Computational grid specified with -n must be greater than one\n"
239 "in each direction\n",stderr);
240 wl.deallocate();
241 return VOROPP_CMD_LINE_ERROR;
242 }
243 } else if(strcmp(argv[i],"-o")==0) {
244 ordered=true;
245 } else if(strcmp(argv[i],"-p")==0) {
246 xperiodic=yperiodic=zperiodic=true;
247 } else if(strcmp(argv[i],"-px")==0) {
248 xperiodic=true;
249 } else if(strcmp(argv[i],"-py")==0) {
250 yperiodic=true;
251 } else if(strcmp(argv[i],"-pz")==0) {
252 zperiodic=true;
253 } else if(strcmp(argv[i],"-r")==0) {
254 polydisperse=true;
255 } else if(strcmp(argv[i],"-v")==0) {
256 verbose=true;
257 } else if(strcmp(argv[i],"--version")==0) {
258 version_message();
259 wl.deallocate();
260 return 0;
261 } else if(strcmp(argv[i],"-wb")==0) {
262 if(i>=argc-13) {error_message();wl.deallocate();return VOROPP_CMD_LINE_ERROR;}
263 i++;
264 double w0=atof(argv[i++]),w1=atof(argv[i++]);
265 double w2=atof(argv[i++]),w3=atof(argv[i++]);
266 double w4=atof(argv[i++]),w5=atof(argv[i]);
267 wl.add_wall(new wall_plane(-1,0,0,-w0,j));j--;
268 wl.add_wall(new wall_plane(1,0,0,w1,j));j--;
269 wl.add_wall(new wall_plane(0,-1,0,-w2,j));j--;
270 wl.add_wall(new wall_plane(0,1,0,w3,j));j--;
271 wl.add_wall(new wall_plane(0,0,-1,-w4,j));j--;
272 wl.add_wall(new wall_plane(0,0,1,w5,j));j--;
273 } else if(strcmp(argv[i],"-ws")==0) {
274 if(i>=argc-11) {error_message();wl.deallocate();return VOROPP_CMD_LINE_ERROR;}
275 i++;
276 double w0=atof(argv[i++]),w1=atof(argv[i++]);
277 double w2=atof(argv[i++]),w3=atof(argv[i]);
278 wl.add_wall(new wall_sphere(w0,w1,w2,w3,j));
279 j--;
280 } else if(strcmp(argv[i],"-wp")==0) {
281 if(i>=argc-11) {error_message();wl.deallocate();return VOROPP_CMD_LINE_ERROR;}
282 i++;
283 double w0=atof(argv[i++]),w1=atof(argv[i++]);
284 double w2=atof(argv[i++]),w3=atof(argv[i]);
285 wl.add_wall(new wall_plane(w0,w1,w2,w3,j));
286 j--;
287 } else if(strcmp(argv[i],"-wc")==0) {
288 if(i>=argc-14) {error_message();wl.deallocate();return VOROPP_CMD_LINE_ERROR;}
289 i++;
290 double w0=atof(argv[i++]),w1=atof(argv[i++]);
291 double w2=atof(argv[i++]),w3=atof(argv[i++]);
292 double w4=atof(argv[i++]),w5=atof(argv[i++]);
293 double w6=atof(argv[i]);
294 wl.add_wall(new wall_cylinder(w0,w1,w2,w3,w4,w5,w6,j));
295 j--;
296 } else if(strcmp(argv[i],"-wo")==0) {
297 if(i>=argc-14) {error_message();wl.deallocate();return VOROPP_CMD_LINE_ERROR;}
298 i++;
299 double w0=atof(argv[i++]),w1=atof(argv[i++]);
300 double w2=atof(argv[i++]),w3=atof(argv[i++]);
301 double w4=atof(argv[i++]),w5=atof(argv[i++]);
302 double w6=atof(argv[i]);
303 wl.add_wall(new wall_cone(w0,w1,w2,w3,w4,w5,w6,j));
304 j--;
305 } else if(strcmp(argv[i],"-y")==0) {
306 povp_output=povv_output=true;
307 } else if(strcmp(argv[i],"-yp")==0) {
308 povp_output=true;
309 } else if(strcmp(argv[i],"-yv")==0) {
310 povv_output=true;
311 } else {
312 wl.deallocate();
313 error_message();
314 return VOROPP_CMD_LINE_ERROR;
315 }
316 i++;
317 }
318
319 // Check the memory guess is positive
320 if(init_mem<=0) {
321 fputs("voro++: The memory allocation must be positive\n",stderr);
322 wl.deallocate();
323 return VOROPP_CMD_LINE_ERROR;
324 }
325
326 // Read in the dimensions of the test box, and estimate the number of
327 // boxes to divide the region up into
328 double ax=atof(argv[i]),bx=atof(argv[i+1]);
329 double ay=atof(argv[i+2]),by=atof(argv[i+3]);
330 double az=atof(argv[i+4]),bz=atof(argv[i+5]);
331
332 // Check that for each coordinate, the minimum value is smaller
333 // than the maximum value
334 if(bx<ax) {
335 fputs("voro++: Minimum x coordinate exceeds maximum x coordinate\n",stderr);
336 wl.deallocate();
337 return VOROPP_CMD_LINE_ERROR;
338 }
339 if(by<ay) {
340 fputs("voro++: Minimum y coordinate exceeds maximum y coordinate\n",stderr);
341 wl.deallocate();
342 return VOROPP_CMD_LINE_ERROR;
343 }
344 if(bz<az) {
345 fputs("voro++: Minimum z coordinate exceeds maximum z coordinate\n",stderr);
346 wl.deallocate();
347 return VOROPP_CMD_LINE_ERROR;
348 }
349
350 if(bm==none) {
351 if(polydisperse) {
352 pconp=new pre_container_poly(ax,bx,ay,by,az,bz,xperiodic,yperiodic,zperiodic);
353 pconp->import(argv[i+6]);
354 pconp->guess_optimal(nx,ny,nz);
355 } else {
356 pcon=new pre_container(ax,bx,ay,by,az,bz,xperiodic,yperiodic,zperiodic);
357 pcon->import(argv[i+6]);
358 pcon->guess_optimal(nx,ny,nz);
359 }
360 } else {
361 double nxf,nyf,nzf;
362 if(bm==length_scale) {
363
364 // Check that the length scale is positive and
365 // reasonably large
366 if(ls<tolerance) {
367 fputs("voro++: ",stderr);
368 if(ls<0) {
369 fputs("The length scale must be positive\n",stderr);
370 } else {
371 fprintf(stderr,"The length scale is smaller than the safe limit of %g. Either\nincrease the particle length scale, or recompile with a different limit.\n",tolerance);
372 }
373 wl.deallocate();
374 return VOROPP_CMD_LINE_ERROR;
375 }
376 ls=0.6/ls;
377 nxf=(bx-ax)*ls+1;
378 nyf=(by-ay)*ls+1;
379 nzf=(bz-az)*ls+1;
380
381 nx=int(nxf);ny=int(nyf);nz=int(nzf);
382 } else {
383 nxf=nx;nyf=ny;nzf=nz;
384 }
385
386 // Compute the number regions based on the length scale
387 // provided. If the total number exceeds a cutoff then bail
388 // out, to prevent making a massive memory allocation. Do this
389 // test using floating point numbers, since huge integers could
390 // potentially wrap around to negative values.
391 if(nxf*nyf*nzf>max_regions) {
392 fprintf(stderr,"voro++: Number of computational blocks exceeds the maximum allowed of %d.\n"
393 "Either increase the particle length scale, or recompile with an increased\nmaximum.",max_regions);
394 wl.deallocate();
395 return VOROPP_MEMORY_ERROR;
396 }
397 }
398
399 // Check that the output filename is a sensible length
400 int flen=strlen(argv[i+6]);
401 if(flen>4096) {
402 fputs("voro++: Filename too long\n",stderr);
403 wl.deallocate();
404 return VOROPP_CMD_LINE_ERROR;
405 }
406
407 // Open files for output
408 char *buffer=new char[flen+7];
409 sprintf(buffer,"%s.vol",argv[i+6]);
410 FILE *outfile=safe_fopen(buffer,"w"),*gnu_file,*povp_file,*povv_file;
411 if(gnuplot_output) {
412 sprintf(buffer,"%s.gnu",argv[i+6]);
413 gnu_file=safe_fopen(buffer,"w");
414 } else gnu_file=NULL;
415 if(povp_output) {
416 sprintf(buffer,"%s_p.pov",argv[i+6]);
417 povp_file=safe_fopen(buffer,"w");
418 } else povp_file=NULL;
419 if(povv_output) {
420 sprintf(buffer,"%s_v.pov",argv[i+6]);
421 povv_file=safe_fopen(buffer,"w");
422 } else povv_file=NULL;
423 delete [] buffer;
424
425 const char *c_str=(custom_output==0?(polydisperse?"%i %q %v %r":"%i %q %v"):argv[custom_output]);
426
427 // Now switch depending on whether polydispersity was enabled, and
428 // whether output ordering is requested
429 double vol=0;int tp=0,vcc=0;
430 if(polydisperse) {
431 if(ordered) {
432 particle_order vo;
433 container_poly con(ax,bx,ay,by,az,bz,nx,ny,nz,xperiodic,yperiodic,zperiodic,init_mem);
434 con.add_wall(wl);
435 if(bm==none) {
436 pconp->setup(vo,con);delete pconp;
437 } else con.import(vo,argv[i+6]);
438
439 c_loop_order vlo(con,vo);
440 cmd_line_output(vlo,con,c_str,outfile,gnu_file,povp_file,povv_file,verbose,vol,vcc,tp);
441 } else {
442 container_poly con(ax,bx,ay,by,az,bz,nx,ny,nz,xperiodic,yperiodic,zperiodic,init_mem);
443 con.add_wall(wl);
444
445 if(bm==none) {
446 pconp->setup(con);delete pconp;
447 } else con.import(argv[i+6]);
448
449 c_loop_all vla(con);
450 cmd_line_output(vla,con,c_str,outfile,gnu_file,povp_file,povv_file,verbose,vol,vcc,tp);
451 }
452 } else {
453 if(ordered) {
454 particle_order vo;
455 container con(ax,bx,ay,by,az,bz,nx,ny,nz,xperiodic,yperiodic,zperiodic,init_mem);
456 con.add_wall(wl);
457 if(bm==none) {
458 pcon->setup(vo,con);delete pcon;
459 } else con.import(vo,argv[i+6]);
460
461 c_loop_order vlo(con,vo);
462 cmd_line_output(vlo,con,c_str,outfile,gnu_file,povp_file,povv_file,verbose,vol,vcc,tp);
463 } else {
464 container con(ax,bx,ay,by,az,bz,nx,ny,nz,xperiodic,yperiodic,zperiodic,init_mem);
465 con.add_wall(wl);
466 if(bm==none) {
467 pcon->setup(con);delete pcon;
468 } else con.import(argv[i+6]);
469 c_loop_all vla(con);
470 cmd_line_output(vla,con,c_str,outfile,gnu_file,povp_file,povv_file,verbose,vol,vcc,tp);
471 }
472 }
473
474 // Print information if verbose output requested
475 if(verbose) {
476 printf("Container geometry : [%g:%g] [%g:%g] [%g:%g]\n"
477 "Computational grid size : %d by %d by %d (%s)\n"
478 "Filename : %s\n"
479 "Output string : %s%s\n",ax,bx,ay,by,az,bz,nx,ny,nz,
480 bm==none?"estimated from file":(bm==length_scale?
481 "estimated using length scale":"directly specified"),
482 argv[i+6],c_str,custom_output==0?" (default)":"");
483 printf("Total imported particles : %d (%.2g per grid block)\n"
484 "Total V. cells computed : %d\n"
485 "Total container volume : %g\n"
486 "Total V. cell volume : %g\n",tp,((double) tp)/(nx*ny*nz),
487 vcc,(bx-ax)*(by-ay)*(bz-az),vol);
488 }
489
490 // Close output files
491 fclose(outfile);
492 if(gnu_file!=NULL) fclose(gnu_file);
493 if(povp_file!=NULL) fclose(povp_file);
494 if(povv_file!=NULL) fclose(povv_file);
495 return 0;
496 }
497
498