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