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 container.hh 8 * \brief Header file for the container_base and related classes. */ 9 10 #ifndef VOROPP_CONTAINER_HH 11 #define VOROPP_CONTAINER_HH 12 13 #include <cstdio> 14 #include <cstdlib> 15 #include <cmath> 16 #include <vector> 17 using namespace std; 18 19 #include "config.hh" 20 #include "common.hh" 21 #include "v_base.hh" 22 #include "cell.hh" 23 #include "c_loops.hh" 24 #include "v_compute.hh" 25 #include "rad_option.hh" 26 27 namespace voro { 28 29 /** \brief Pure virtual class from which wall objects are derived. 30 * 31 * This is a pure virtual class for a generic wall object. A wall object 32 * can be specified by deriving a new class from this and specifying the 33 * functions.*/ 34 class wall { 35 public: ~wall()36 virtual ~wall() {} 37 /** A pure virtual function for testing whether a point is 38 * inside the wall object. */ 39 virtual bool point_inside(double x,double y,double z) = 0; 40 /** A pure virtual function for cutting a cell without 41 * neighbor-tracking with a wall. */ 42 virtual bool cut_cell(voronoicell &c,double x,double y,double z) = 0; 43 /** A pure virtual function for cutting a cell with 44 * neighbor-tracking enabled with a wall. */ 45 virtual bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) = 0; 46 }; 47 48 /** \brief A class for storing a list of pointers to walls. 49 * 50 * This class stores a list of pointers to wall classes. It contains several 51 * simple routines that make use of the wall classes (such as telling whether a 52 * given position is inside all of the walls or not). It can be used by itself, 53 * but also forms part of container_base, for associating walls with this 54 * class. */ 55 class wall_list { 56 public: 57 /** An array holding pointers to wall objects. */ 58 wall **walls; 59 /** A pointer to the next free position to add a wall pointer. 60 */ 61 wall **wep; 62 wall_list(); 63 ~wall_list(); 64 /** Adds a wall to the list. 65 * \param[in] w the wall to add. */ add_wall(wall * w)66 inline void add_wall(wall *w) { 67 if(wep==wel) increase_wall_memory(); 68 *(wep++)=w; 69 } 70 /** Adds a wall to the list. 71 * \param[in] w a reference to the wall to add. */ add_wall(wall & w)72 inline void add_wall(wall &w) {add_wall(&w);} 73 void add_wall(wall_list &wl); 74 /** Determines whether a given position is inside all of the 75 * walls on the list. 76 * \param[in] (x,y,z) the position to test. 77 * \return True if it is inside, false if it is outside. */ point_inside_walls(double x,double y,double z)78 inline bool point_inside_walls(double x,double y,double z) { 79 for(wall **wp=walls;wp<wep;wp++) if(!((*wp)->point_inside(x,y,z))) return false; 80 return true; 81 } 82 /** Cuts a Voronoi cell by all of the walls currently on 83 * the list. 84 * \param[in] c a reference to the Voronoi cell class. 85 * \param[in] (x,y,z) the position of the cell. 86 * \return True if the cell still exists, false if the cell is 87 * deleted. */ 88 template<class c_class> apply_walls(c_class & c,double x,double y,double z)89 bool apply_walls(c_class &c,double x,double y,double z) { 90 for(wall **wp=walls;wp<wep;wp++) if(!((*wp)->cut_cell(c,x,y,z))) return false; 91 return true; 92 } 93 void deallocate(); 94 protected: 95 void increase_wall_memory(); 96 /** A pointer to the limit of the walls array, used to 97 * determine when array is full. */ 98 wall **wel; 99 /** The current amount of memory allocated for walls. */ 100 int current_wall_size; 101 }; 102 103 /** \brief Class for representing a particle system in a three-dimensional 104 * rectangular box. 105 * 106 * This class represents a system of particles in a three-dimensional 107 * rectangular box. Any combination of non-periodic and periodic coordinates 108 * can be used in the three coordinate directions. The class is not intended 109 * for direct use, but instead forms the base of the container and 110 * container_poly classes that add specialized routines for computing the 111 * regular and radical Voronoi tessellations respectively. It contains routines 112 * that are commonly between these two classes, such as those for drawing the 113 * domain, and placing particles within the internal data structure. 114 * 115 * The class is derived from the wall_list class, which encapsulates routines 116 * for associating walls with the container, and the voro_base class, which 117 * encapsulates routines about the underlying computational grid. */ 118 class container_base : public voro_base, public wall_list { 119 public: 120 /** The minimum x coordinate of the container. */ 121 const double ax; 122 /** The maximum x coordinate of the container. */ 123 const double bx; 124 /** The minimum y coordinate of the container. */ 125 const double ay; 126 /** The maximum y coordinate of the container. */ 127 const double by; 128 /** The minimum z coordinate of the container. */ 129 const double az; 130 /** The maximum z coordinate of the container. */ 131 const double bz; 132 /** A boolean value that determines if the x coordinate in 133 * periodic or not. */ 134 const bool xperiodic; 135 /** A boolean value that determines if the y coordinate in 136 * periodic or not. */ 137 const bool yperiodic; 138 /** A boolean value that determines if the z coordinate in 139 * periodic or not. */ 140 const bool zperiodic; 141 /** This array holds the numerical IDs of each particle in each 142 * computational box. */ 143 int **id; 144 /** A two dimensional array holding particle positions. For the 145 * derived container_poly class, this also holds particle 146 * radii. */ 147 double **p; 148 /** This array holds the number of particles within each 149 * computational box of the container. */ 150 int *co; 151 /** This array holds the maximum amount of particle memory for 152 * each computational box of the container. If the number of 153 * particles in a particular box ever approaches this limit, 154 * more is allocated using the add_particle_memory() function. 155 */ 156 int *mem; 157 /** The amount of memory in the array structure for each 158 * particle. This is set to 3 when the basic class is 159 * initialized, so that the array holds (x,y,z) positions. If 160 * the container class is initialized as part of the derived 161 * class container_poly, then this is set to 4, to also hold 162 * the particle radii. */ 163 const int ps; 164 container_base(double ax_,double bx_,double ay_,double by_,double az_,double bz_, 165 int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_, 166 int init_mem,int ps_); 167 ~container_base(); 168 bool point_inside(double x,double y,double z); 169 void region_count(); 170 /** Initializes the Voronoi cell prior to a compute_cell 171 * operation for a specific particle being carried out by a 172 * voro_compute class. The cell is initialized to fill the 173 * entire container. For non-periodic coordinates, this is set 174 * by the position of the walls. For periodic coordinates, the 175 * space is equally divided in either direction from the 176 * particle's initial position. Plane cuts made by any walls 177 * that have been added are then applied to the cell. 178 * \param[in,out] c a reference to a voronoicell object. 179 * \param[in] ijk the block that the particle is within. 180 * \param[in] q the index of the particle within its block. 181 * \param[in] (ci,cj,ck) the coordinates of the block in the 182 * container coordinate system. 183 * \param[out] (i,j,k) the coordinates of the test block 184 * relative to the voro_compute 185 * coordinate system. 186 * \param[out] (x,y,z) the position of the particle. 187 * \param[out] disp a block displacement used internally by the 188 * compute_cell routine. 189 * \return False if the plane cuts applied by walls completely 190 * removed the cell, true otherwise. */ 191 template<class v_cell> initialize_voronoicell(v_cell & c,int ijk,int q,int ci,int cj,int ck,int & i,int & j,int & k,double & x,double & y,double & z,int & disp)192 inline bool initialize_voronoicell(v_cell &c,int ijk,int q,int ci,int cj,int ck, 193 int &i,int &j,int &k,double &x,double &y,double &z,int &disp) { 194 double x1,x2,y1,y2,z1,z2,*pp=p[ijk]+ps*q; 195 x=*(pp++);y=*(pp++);z=*pp; 196 if(xperiodic) {x1=-(x2=0.5*(bx-ax));i=nx;} else {x1=ax-x;x2=bx-x;i=ci;} 197 if(yperiodic) {y1=-(y2=0.5*(by-ay));j=ny;} else {y1=ay-y;y2=by-y;j=cj;} 198 if(zperiodic) {z1=-(z2=0.5*(bz-az));k=nz;} else {z1=az-z;z2=bz-z;k=ck;} 199 c.init(x1,x2,y1,y2,z1,z2); 200 if(!apply_walls(c,x,y,z)) return false; 201 disp=ijk-i-nx*(j+ny*k); 202 return true; 203 } 204 /** Initializes parameters for a find_voronoi_cell call within 205 * the voro_compute template. 206 * \param[in] (ci,cj,ck) the coordinates of the test block in 207 * the container coordinate system. 208 * \param[in] ijk the index of the test block 209 * \param[out] (i,j,k) the coordinates of the test block 210 * relative to the voro_compute 211 * coordinate system. 212 * \param[out] disp a block displacement used internally by the 213 * find_voronoi_cell routine. */ initialize_search(int ci,int cj,int ck,int ijk,int & i,int & j,int & k,int & disp)214 inline void initialize_search(int ci,int cj,int ck,int ijk,int &i,int &j,int &k,int &disp) { 215 i=xperiodic?nx:ci; 216 j=yperiodic?ny:cj; 217 k=zperiodic?nz:ck; 218 disp=ijk-i-nx*(j+ny*k); 219 } 220 /** Returns the position of a particle currently being computed 221 * relative to the computational block that it is within. It is 222 * used to select the optimal worklist entry to use. 223 * \param[in] (x,y,z) the position of the particle. 224 * \param[in] (ci,cj,ck) the block that the particle is within. 225 * \param[out] (fx,fy,fz) the position relative to the block. 226 */ frac_pos(double x,double y,double z,double ci,double cj,double ck,double & fx,double & fy,double & fz)227 inline void frac_pos(double x,double y,double z,double ci,double cj,double ck, 228 double &fx,double &fy,double &fz) { 229 fx=x-ax-boxx*ci; 230 fy=y-ay-boxy*cj; 231 fz=z-az-boxz*ck; 232 } 233 /** Calculates the index of block in the container structure 234 * corresponding to given coordinates. 235 * \param[in] (ci,cj,ck) the coordinates of the original block 236 * in the current computation, relative 237 * to the container coordinate system. 238 * \param[in] (ei,ej,ek) the displacement of the current block 239 * from the original block. 240 * \param[in,out] (qx,qy,qz) the periodic displacement that 241 * must be added to the particles 242 * within the computed block. 243 * \param[in] disp a block displacement used internally by the 244 * find_voronoi_cell and compute_cell routines. 245 * \return The block index. */ region_index(int ci,int cj,int ck,int ei,int ej,int ek,double & qx,double & qy,double & qz,int & disp)246 inline int region_index(int ci,int cj,int ck,int ei,int ej,int ek,double &qx,double &qy,double &qz,int &disp) { 247 if(xperiodic) {if(ci+ei<nx) {ei+=nx;qx=-(bx-ax);} else if(ci+ei>=(nx<<1)) {ei-=nx;qx=bx-ax;} else qx=0;} 248 if(yperiodic) {if(cj+ej<ny) {ej+=ny;qy=-(by-ay);} else if(cj+ej>=(ny<<1)) {ej-=ny;qy=by-ay;} else qy=0;} 249 if(zperiodic) {if(ck+ek<nz) {ek+=nz;qz=-(bz-az);} else if(ck+ek>=(nz<<1)) {ek-=nz;qz=bz-az;} else qz=0;} 250 return disp+ei+nx*(ej+ny*ek); 251 } 252 void draw_domain_gnuplot(FILE *fp=stdout); 253 /** Draws an outline of the domain in Gnuplot format. 254 * \param[in] filename the filename to write to. */ draw_domain_gnuplot(const char * filename)255 inline void draw_domain_gnuplot(const char* filename) { 256 FILE *fp=safe_fopen(filename,"w"); 257 draw_domain_gnuplot(fp); 258 fclose(fp); 259 } 260 void draw_domain_pov(FILE *fp=stdout); 261 /** Draws an outline of the domain in Gnuplot format. 262 * \param[in] filename the filename to write to. */ draw_domain_pov(const char * filename)263 inline void draw_domain_pov(const char* filename) { 264 FILE *fp=safe_fopen(filename,"w"); 265 draw_domain_pov(fp); 266 fclose(fp); 267 } 268 /** Sums up the total number of stored particles. 269 * \return The number of particles. */ total_particles()270 inline int total_particles() { 271 int tp=*co; 272 for(int *cop=co+1;cop<co+nxyz;cop++) tp+=*cop; 273 return tp; 274 } 275 protected: 276 void add_particle_memory(int i); 277 inline bool put_locate_block(int &ijk,double &x,double &y,double &z); 278 inline bool put_remap(int &ijk,double &x,double &y,double &z); 279 inline bool remap(int &ai,int &aj,int &ak,int &ci,int &cj,int &ck,double &x,double &y,double &z,int &ijk); 280 }; 281 282 /** \brief Extension of the container_base class for computing regular Voronoi 283 * tessellations. 284 * 285 * This class is an extension of the container_base class that has routines 286 * specifically for computing the regular Voronoi tessellation with no 287 * dependence on particle radii. */ 288 class container : public container_base, public radius_mono { 289 public: 290 container(double ax_,double bx_,double ay_,double by_,double az_,double bz_, 291 int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,int init_mem); 292 void clear(); 293 void put(int n,double x,double y,double z); 294 void put(particle_order &vo,int n,double x,double y,double z); 295 void import(FILE *fp=stdin); 296 void import(particle_order &vo,FILE *fp=stdin); 297 /** Imports a list of particles from an open file stream into 298 * the container. Entries of four numbers (Particle ID, x 299 * position, y position, z position) are searched for. If the 300 * file cannot be successfully read, then the routine causes a 301 * fatal error. 302 * \param[in] filename the name of the file to open and read 303 * from. */ import(const char * filename)304 inline void import(const char* filename) { 305 FILE *fp=safe_fopen(filename,"r"); 306 import(fp); 307 fclose(fp); 308 } 309 /** Imports a list of particles from an open file stream into 310 * the container. Entries of four numbers (Particle ID, x 311 * position, y position, z position) are searched for. In 312 * addition, the order in which particles are read is saved 313 * into an ordering class. If the file cannot be successfully 314 * read, then the routine causes a fatal error. 315 * \param[in,out] vo the ordering class to use. 316 * \param[in] filename the name of the file to open and read 317 * from. */ import(particle_order & vo,const char * filename)318 inline void import(particle_order &vo,const char* filename) { 319 FILE *fp=safe_fopen(filename,"r"); 320 import(vo,fp); 321 fclose(fp); 322 } 323 void compute_all_cells(); 324 double sum_cell_volumes(); 325 /** Dumps particle IDs and positions to a file. 326 * \param[in] vl the loop class to use. 327 * \param[in] fp a file handle to write to. */ 328 template<class c_loop> draw_particles(c_loop & vl,FILE * fp)329 void draw_particles(c_loop &vl,FILE *fp) { 330 double *pp; 331 if(vl.start()) do { 332 pp=p[vl.ijk]+3*vl.q; 333 fprintf(fp,"%d %g %g %g\n",id[vl.ijk][vl.q],*pp,pp[1],pp[2]); 334 } while(vl.inc()); 335 } 336 /** Dumps all of the particle IDs and positions to a file. 337 * \param[in] fp a file handle to write to. */ draw_particles(FILE * fp=stdout)338 inline void draw_particles(FILE *fp=stdout) { 339 c_loop_all vl(*this); 340 draw_particles(vl,fp); 341 } 342 /** Dumps all of the particle IDs and positions to a file. 343 * \param[in] filename the name of the file to write to. */ draw_particles(const char * filename)344 inline void draw_particles(const char *filename) { 345 FILE *fp=safe_fopen(filename,"w"); 346 draw_particles(fp); 347 fclose(fp); 348 } 349 /** Dumps particle positions in POV-Ray format. 350 * \param[in] vl the loop class to use. 351 * \param[in] fp a file handle to write to. */ 352 template<class c_loop> draw_particles_pov(c_loop & vl,FILE * fp)353 void draw_particles_pov(c_loop &vl,FILE *fp) { 354 double *pp; 355 if(vl.start()) do { 356 pp=p[vl.ijk]+3*vl.q; 357 fprintf(fp,"// id %d\nsphere{<%g,%g,%g>,s}\n", 358 id[vl.ijk][vl.q],*pp,pp[1],pp[2]); 359 } while(vl.inc()); 360 } 361 /** Dumps all particle positions in POV-Ray format. 362 * \param[in] fp a file handle to write to. */ draw_particles_pov(FILE * fp=stdout)363 inline void draw_particles_pov(FILE *fp=stdout) { 364 c_loop_all vl(*this); 365 draw_particles_pov(vl,fp); 366 } 367 /** Dumps all particle positions in POV-Ray format. 368 * \param[in] filename the name of the file to write to. */ draw_particles_pov(const char * filename)369 inline void draw_particles_pov(const char *filename) { 370 FILE *fp=safe_fopen(filename,"w"); 371 draw_particles_pov(fp); 372 fclose(fp); 373 } 374 /** Computes Voronoi cells and saves the output in gnuplot 375 * format. 376 * \param[in] vl the loop class to use. 377 * \param[in] fp a file handle to write to. */ 378 template<class c_loop> draw_cells_gnuplot(c_loop & vl,FILE * fp)379 void draw_cells_gnuplot(c_loop &vl,FILE *fp) { 380 voronoicell c;double *pp; 381 if(vl.start()) do if(compute_cell(c,vl)) { 382 pp=p[vl.ijk]+ps*vl.q; 383 c.draw_gnuplot(*pp,pp[1],pp[2],fp); 384 } while(vl.inc()); 385 } 386 /** Computes all Voronoi cells and saves the output in gnuplot 387 * format. 388 * \param[in] fp a file handle to write to. */ draw_cells_gnuplot(FILE * fp=stdout)389 inline void draw_cells_gnuplot(FILE *fp=stdout) { 390 c_loop_all vl(*this); 391 draw_cells_gnuplot(vl,fp); 392 } 393 /** Computes all Voronoi cells and saves the output in gnuplot 394 * format. 395 * \param[in] filename the name of the file to write to. */ draw_cells_gnuplot(const char * filename)396 inline void draw_cells_gnuplot(const char *filename) { 397 FILE *fp=safe_fopen(filename,"w"); 398 draw_cells_gnuplot(fp); 399 fclose(fp); 400 } 401 /** Computes Voronoi cells and saves the output in POV-Ray 402 * format. 403 * \param[in] vl the loop class to use. 404 * \param[in] fp a file handle to write to. */ 405 template<class c_loop> draw_cells_pov(c_loop & vl,FILE * fp)406 void draw_cells_pov(c_loop &vl,FILE *fp) { 407 voronoicell c;double *pp; 408 if(vl.start()) do if(compute_cell(c,vl)) { 409 fprintf(fp,"// cell %d\n",id[vl.ijk][vl.q]); 410 pp=p[vl.ijk]+ps*vl.q; 411 c.draw_pov(*pp,pp[1],pp[2],fp); 412 } while(vl.inc()); 413 } 414 /** Computes all Voronoi cells and saves the output in POV-Ray 415 * format. 416 * \param[in] fp a file handle to write to. */ draw_cells_pov(FILE * fp=stdout)417 inline void draw_cells_pov(FILE *fp=stdout) { 418 c_loop_all vl(*this); 419 draw_cells_pov(vl,fp); 420 } 421 /** Computes all Voronoi cells and saves the output in POV-Ray 422 * format. 423 * \param[in] filename the name of the file to write to. */ draw_cells_pov(const char * filename)424 inline void draw_cells_pov(const char *filename) { 425 FILE *fp=safe_fopen(filename,"w"); 426 draw_cells_pov(fp); 427 fclose(fp); 428 } 429 /** Computes the Voronoi cells and saves customized information 430 * about them. 431 * \param[in] vl the loop class to use. 432 * \param[in] format the custom output string to use. 433 * \param[in] fp a file handle to write to. */ 434 template<class c_loop> print_custom(c_loop & vl,const char * format,FILE * fp)435 void print_custom(c_loop &vl,const char *format,FILE *fp) { 436 int ijk,q;double *pp; 437 if(contains_neighbor(format)) { 438 voronoicell_neighbor c; 439 if(vl.start()) do if(compute_cell(c,vl)) { 440 ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q; 441 c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],default_radius,fp); 442 } while(vl.inc()); 443 } else { 444 voronoicell c; 445 if(vl.start()) do if(compute_cell(c,vl)) { 446 ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q; 447 c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],default_radius,fp); 448 } while(vl.inc()); 449 } 450 } 451 void print_custom(const char *format,FILE *fp=stdout); 452 void print_custom(const char *format,const char *filename); 453 bool find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid); 454 /** Computes the Voronoi cell for a particle currently being 455 * referenced by a loop class. 456 * \param[out] c a Voronoi cell class in which to store the 457 * computed cell. 458 * \param[in] vl the loop class to use. 459 * \return True if the cell was computed. If the cell cannot be 460 * computed, if it is removed entirely by a wall or boundary 461 * condition, then the routine returns false. */ 462 template<class v_cell,class c_loop> compute_cell(v_cell & c,c_loop & vl)463 inline bool compute_cell(v_cell &c,c_loop &vl) { 464 return vc.compute_cell(c,vl.ijk,vl.q,vl.i,vl.j,vl.k); 465 } 466 /** Computes the Voronoi cell for given particle. 467 * \param[out] c a Voronoi cell class in which to store the 468 * computed cell. 469 * \param[in] ijk the block that the particle is within. 470 * \param[in] q the index of the particle within the block. 471 * \return True if the cell was computed. If the cell cannot be 472 * computed, if it is removed entirely by a wall or boundary 473 * condition, then the routine returns false. */ 474 template<class v_cell> compute_cell(v_cell & c,int ijk,int q)475 inline bool compute_cell(v_cell &c,int ijk,int q) { 476 int k=ijk/nxy,ijkt=ijk-nxy*k,j=ijkt/nx,i=ijkt-j*nx; 477 return vc.compute_cell(c,ijk,q,i,j,k); 478 } 479 private: 480 voro_compute<container> vc; 481 friend class voro_compute<container>; 482 }; 483 484 /** \brief Extension of the container_base class for computing radical Voronoi 485 * tessellations. 486 * 487 * This class is an extension of container_base class that has routines 488 * specifically for computing the radical Voronoi tessellation that depends on 489 * the particle radii. */ 490 class container_poly : public container_base, public radius_poly { 491 public: 492 container_poly(double ax_,double bx_,double ay_,double by_,double az_,double bz_, 493 int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,int init_mem); 494 void clear(); 495 void put(int n,double x,double y,double z,double r); 496 void put(particle_order &vo,int n,double x,double y,double z,double r); 497 void import(FILE *fp=stdin); 498 void import(particle_order &vo,FILE *fp=stdin); 499 /** Imports a list of particles from an open file stream into 500 * the container_poly class. Entries of five numbers (Particle 501 * ID, x position, y position, z position, radius) are searched 502 * for. If the file cannot be successfully read, then the 503 * routine causes a fatal error. 504 * \param[in] filename the name of the file to open and read 505 * from. */ import(const char * filename)506 inline void import(const char* filename) { 507 FILE *fp=safe_fopen(filename,"r"); 508 import(fp); 509 fclose(fp); 510 } 511 /** Imports a list of particles from an open file stream into 512 * the container_poly class. Entries of five numbers (Particle 513 * ID, x position, y position, z position, radius) are searched 514 * for. In addition, the order in which particles are read is 515 * saved into an ordering class. If the file cannot be 516 * successfully read, then the routine causes a fatal error. 517 * \param[in,out] vo the ordering class to use. 518 * \param[in] filename the name of the file to open and read 519 * from. */ import(particle_order & vo,const char * filename)520 inline void import(particle_order &vo,const char* filename) { 521 FILE *fp=safe_fopen(filename,"r"); 522 import(vo,fp); 523 fclose(fp); 524 } 525 void compute_all_cells(); 526 double sum_cell_volumes(); 527 /** Dumps particle IDs, positions and radii to a file. 528 * \param[in] vl the loop class to use. 529 * \param[in] fp a file handle to write to. */ 530 template<class c_loop> draw_particles(c_loop & vl,FILE * fp)531 void draw_particles(c_loop &vl,FILE *fp) { 532 double *pp; 533 if(vl.start()) do { 534 pp=p[vl.ijk]+4*vl.q; 535 fprintf(fp,"%d %g %g %g %g\n",id[vl.ijk][vl.q],*pp,pp[1],pp[2],pp[3]); 536 } while(vl.inc()); 537 } 538 /** Dumps all of the particle IDs, positions and radii to a 539 * file. 540 * \param[in] fp a file handle to write to. */ draw_particles(FILE * fp=stdout)541 inline void draw_particles(FILE *fp=stdout) { 542 c_loop_all vl(*this); 543 draw_particles(vl,fp); 544 } 545 /** Dumps all of the particle IDs, positions and radii to a 546 * file. 547 * \param[in] filename the name of the file to write to. */ draw_particles(const char * filename)548 inline void draw_particles(const char *filename) { 549 FILE *fp=safe_fopen(filename,"w"); 550 draw_particles(fp); 551 fclose(fp); 552 } 553 /** Dumps particle positions in POV-Ray format. 554 * \param[in] vl the loop class to use. 555 * \param[in] fp a file handle to write to. */ 556 template<class c_loop> draw_particles_pov(c_loop & vl,FILE * fp)557 void draw_particles_pov(c_loop &vl,FILE *fp) { 558 double *pp; 559 if(vl.start()) do { 560 pp=p[vl.ijk]+4*vl.q; 561 fprintf(fp,"// id %d\nsphere{<%g,%g,%g>,%g}\n", 562 id[vl.ijk][vl.q],*pp,pp[1],pp[2],pp[3]); 563 } while(vl.inc()); 564 } 565 /** Dumps all the particle positions in POV-Ray format. 566 * \param[in] fp a file handle to write to. */ draw_particles_pov(FILE * fp=stdout)567 inline void draw_particles_pov(FILE *fp=stdout) { 568 c_loop_all vl(*this); 569 draw_particles_pov(vl,fp); 570 } 571 /** Dumps all the particle positions in POV-Ray format. 572 * \param[in] filename the name of the file to write to. */ draw_particles_pov(const char * filename)573 inline void draw_particles_pov(const char *filename) { 574 FILE *fp=safe_fopen(filename,"w"); 575 draw_particles_pov(fp); 576 fclose(fp); 577 } 578 /** Computes Voronoi cells and saves the output in gnuplot 579 * format. 580 * \param[in] vl the loop class to use. 581 * \param[in] fp a file handle to write to. */ 582 template<class c_loop> draw_cells_gnuplot(c_loop & vl,FILE * fp)583 void draw_cells_gnuplot(c_loop &vl,FILE *fp) { 584 voronoicell c;double *pp; 585 if(vl.start()) do if(compute_cell(c,vl)) { 586 pp=p[vl.ijk]+ps*vl.q; 587 c.draw_gnuplot(*pp,pp[1],pp[2],fp); 588 } while(vl.inc()); 589 } 590 /** Compute all Voronoi cells and saves the output in gnuplot 591 * format. 592 * \param[in] fp a file handle to write to. */ draw_cells_gnuplot(FILE * fp=stdout)593 inline void draw_cells_gnuplot(FILE *fp=stdout) { 594 c_loop_all vl(*this); 595 draw_cells_gnuplot(vl,fp); 596 } 597 /** Compute all Voronoi cells and saves the output in gnuplot 598 * format. 599 * \param[in] filename the name of the file to write to. */ draw_cells_gnuplot(const char * filename)600 inline void draw_cells_gnuplot(const char *filename) { 601 FILE *fp=safe_fopen(filename,"w"); 602 draw_cells_gnuplot(fp); 603 fclose(fp); 604 } 605 /** Computes Voronoi cells and saves the output in POV-Ray 606 * format. 607 * \param[in] vl the loop class to use. 608 * \param[in] fp a file handle to write to. */ 609 template<class c_loop> draw_cells_pov(c_loop & vl,FILE * fp)610 void draw_cells_pov(c_loop &vl,FILE *fp) { 611 voronoicell c;double *pp; 612 if(vl.start()) do if(compute_cell(c,vl)) { 613 fprintf(fp,"// cell %d\n",id[vl.ijk][vl.q]); 614 pp=p[vl.ijk]+ps*vl.q; 615 c.draw_pov(*pp,pp[1],pp[2],fp); 616 } while(vl.inc()); 617 } 618 /** Computes all Voronoi cells and saves the output in POV-Ray 619 * format. 620 * \param[in] fp a file handle to write to. */ draw_cells_pov(FILE * fp=stdout)621 inline void draw_cells_pov(FILE *fp=stdout) { 622 c_loop_all vl(*this); 623 draw_cells_pov(vl,fp); 624 } 625 /** Computes all Voronoi cells and saves the output in POV-Ray 626 * format. 627 * \param[in] filename the name of the file to write to. */ draw_cells_pov(const char * filename)628 inline void draw_cells_pov(const char *filename) { 629 FILE *fp=safe_fopen(filename,"w"); 630 draw_cells_pov(fp); 631 fclose(fp); 632 } 633 /** Computes the Voronoi cells and saves customized information 634 * about them. 635 * \param[in] vl the loop class to use. 636 * \param[in] format the custom output string to use. 637 * \param[in] fp a file handle to write to. */ 638 template<class c_loop> print_custom(c_loop & vl,const char * format,FILE * fp)639 void print_custom(c_loop &vl,const char *format,FILE *fp) { 640 int ijk,q;double *pp; 641 if(contains_neighbor(format)) { 642 voronoicell_neighbor c; 643 if(vl.start()) do if(compute_cell(c,vl)) { 644 ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q; 645 c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],pp[3],fp); 646 } while(vl.inc()); 647 } else { 648 voronoicell c; 649 if(vl.start()) do if(compute_cell(c,vl)) { 650 ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q; 651 c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],pp[3],fp); 652 } while(vl.inc()); 653 } 654 } 655 /** Computes the Voronoi cell for a particle currently being 656 * referenced by a loop class. 657 * \param[out] c a Voronoi cell class in which to store the 658 * computed cell. 659 * \param[in] vl the loop class to use. 660 * \return True if the cell was computed. If the cell cannot be 661 * computed, if it is removed entirely by a wall or boundary 662 * condition, then the routine returns false. */ 663 template<class v_cell,class c_loop> compute_cell(v_cell & c,c_loop & vl)664 inline bool compute_cell(v_cell &c,c_loop &vl) { 665 return vc.compute_cell(c,vl.ijk,vl.q,vl.i,vl.j,vl.k); 666 } 667 /** Computes the Voronoi cell for given particle. 668 * \param[out] c a Voronoi cell class in which to store the 669 * computed cell. 670 * \param[in] ijk the block that the particle is within. 671 * \param[in] q the index of the particle within the block. 672 * \return True if the cell was computed. If the cell cannot be 673 * computed, if it is removed entirely by a wall or boundary 674 * condition, then the routine returns false. */ 675 template<class v_cell> compute_cell(v_cell & c,int ijk,int q)676 inline bool compute_cell(v_cell &c,int ijk,int q) { 677 int k=ijk/nxy,ijkt=ijk-nxy*k,j=ijkt/nx,i=ijkt-j*nx; 678 return vc.compute_cell(c,ijk,q,i,j,k); 679 } 680 void print_custom(const char *format,FILE *fp=stdout); 681 void print_custom(const char *format,const char *filename); 682 bool find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid); 683 private: 684 voro_compute<container_poly> vc; 685 friend class voro_compute<container_poly>; 686 }; 687 688 } 689 690 #endif 691