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_prd.hh 8 * \brief Header file for the container_periodic_base and related classes. */ 9 10 #ifndef VOROPP_CONTAINER_PRD_HH 11 #define VOROPP_CONTAINER_PRD_HH 12 13 #include <cstdio> 14 #include <vector> 15 16 #include "config.hh" 17 #include "common.hh" 18 #include "v_base.hh" 19 #include "cell.hh" 20 #include "c_loops.hh" 21 #include "v_compute.hh" 22 #include "unitcell.hh" 23 #include "rad_option.hh" 24 25 namespace voro { 26 27 /** \brief Class for representing a particle system in a 3D periodic 28 * non-orthogonal periodic domain. 29 * 30 * This class represents a particle system in a three-dimensional 31 * non-orthogonal periodic domain. The domain is defined by three periodicity 32 * vectors (bx,0,0), (bxy,by,0), and (bxz,byz,bz) that represent a 33 * parallelepiped. Internally, the class stores particles in the box 0<x<bx, 34 * 0<y<by, 0<z<bz, and constructs periodic images of particles that displaced 35 * by the three periodicity vectors when they are necessary for the 36 * computation. The internal memory structure for this class is significantly 37 * different from the container_base class in order to handle the dynamic 38 * construction of these periodic images. 39 * 40 * The class is derived from the unitcell class, which encapsulates information 41 * about the domain geometry, and the voro_base class, which encapsulates 42 * information about the underlying computational grid. */ 43 class container_periodic_base : public unitcell, public voro_base { 44 public: 45 /** The lower y index (inclusive) of the primary domain within 46 * the block structure. */ 47 int ey; 48 /** The lower z index (inclusive) of the primary domain within 49 * the block structure. */ 50 int ez; 51 /** The upper y index (exclusive) of the primary domain within 52 * the block structure. */ 53 int wy; 54 /** The upper z index (exclusive) of the primary domain within 55 * the block structure. */ 56 int wz; 57 /** The total size of the block structure (including images) in 58 * the y direction. */ 59 int oy; 60 /** The total size of the block structure (including images) in 61 * the z direction. */ 62 int oz; 63 /** The total number of blocks. */ 64 int oxyz; 65 /** This array holds the numerical IDs of each particle in each 66 * computational box. */ 67 int **id; 68 /** A two dimensional array holding particle positions. For the 69 * derived container_poly class, this also holds particle 70 * radii. */ 71 double **p; 72 /** This array holds the number of particles within each 73 * computational box of the container. */ 74 int *co; 75 /** This array holds the maximum amount of particle memory for 76 * each computational box of the container. If the number of 77 * particles in a particular box ever approaches this limit, 78 * more is allocated using the add_particle_memory() function. 79 */ 80 int *mem; 81 /** An array holding information about periodic image 82 * construction at a given location. */ 83 char *img; 84 /** The initial amount of memory to allocate for particles 85 * for each block. */ 86 const int init_mem; 87 /** The amount of memory in the array structure for each 88 * particle. This is set to 3 when the basic class is 89 * initialized, so that the array holds (x,y,z) positions. If 90 * the container class is initialized as part of the derived 91 * class container_poly, then this is set to 4, to also hold 92 * the particle radii. */ 93 const int ps; 94 container_periodic_base(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_, 95 int nx_,int ny_,int nz_,int init_mem_,int ps); 96 ~container_periodic_base(); 97 /** Prints all particles in the container, including those that 98 * have been constructed in image blocks. */ print_all_particles()99 inline void print_all_particles() { 100 int ijk,q; 101 for(ijk=0;ijk<oxyz;ijk++) for(q=0;q<co[ijk];q++) 102 printf("%d %g %g %g\n",id[ijk][q],p[ijk][ps*q],p[ijk][ps*q+1],p[ijk][ps*q+2]); 103 } 104 void region_count(); 105 /** Initializes the Voronoi cell prior to a compute_cell 106 * operation for a specific particle being carried out by a 107 * voro_compute class. The cell is initialized to be the 108 * pre-computed unit Voronoi cell based on planes formed by 109 * periodic images of the particle. 110 * \param[in,out] c a reference to a voronoicell object. 111 * \param[in] ijk the block that the particle is within. 112 * \param[in] q the index of the particle within its block. 113 * \param[in] (ci,cj,ck) the coordinates of the block in the 114 * container coordinate system. 115 * \param[out] (i,j,k) the coordinates of the test block 116 * relative to the voro_compute 117 * coordinate system. 118 * \param[out] (x,y,z) the position of the particle. 119 * \param[out] disp a block displacement used internally by the 120 * compute_cell routine. 121 * \return False if the plane cuts applied by walls completely 122 * removed the cell, true otherwise. */ 123 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)124 inline bool 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) { 125 c=unit_voro; 126 double *pp=p[ijk]+ps*q; 127 x=*(pp++);y=*(pp++);z=*pp; 128 i=nx;j=ey;k=ez; 129 return true; 130 } 131 /** Initializes parameters for a find_voronoi_cell call within 132 * the voro_compute template. 133 * \param[in] (ci,cj,ck) the coordinates of the test block in 134 * the container coordinate system. 135 * \param[in] ijk the index of the test block 136 * \param[out] (i,j,k) the coordinates of the test block 137 * relative to the voro_compute 138 * coordinate system. 139 * \param[out] disp a block displacement used internally by the 140 * find_voronoi_cell routine (but not needed 141 * in this instance.) */ initialize_search(int ci,int cj,int ck,int ijk,int & i,int & j,int & k,int & disp)142 inline void initialize_search(int ci,int cj,int ck,int ijk,int &i,int &j,int &k,int &disp) { 143 i=nx;j=ey;k=ez; 144 } 145 /** Returns the position of a particle currently being computed 146 * relative to the computational block that it is within. It is 147 * used to select the optimal worklist entry to use. 148 * \param[in] (x,y,z) the position of the particle. 149 * \param[in] (ci,cj,ck) the block that the particle is within. 150 * \param[out] (fx,fy,fz) the position relative to the block. 151 */ frac_pos(double x,double y,double z,double ci,double cj,double ck,double & fx,double & fy,double & fz)152 inline void frac_pos(double x,double y,double z,double ci,double cj,double ck,double &fx,double &fy,double &fz) { 153 fx=x-boxx*ci; 154 fy=y-boxy*(cj-ey); 155 fz=z-boxz*(ck-ez); 156 } 157 /** Calculates the index of block in the container structure 158 * corresponding to given coordinates. 159 * \param[in] (ci,cj,ck) the coordinates of the original block 160 * in the current computation, relative 161 * to the container coordinate system. 162 * \param[in] (ei,ej,ek) the displacement of the current block 163 * from the original block. 164 * \param[in,out] (qx,qy,qz) the periodic displacement that 165 * must be added to the particles 166 * within the computed block. 167 * \param[in] disp a block displacement used internally by the 168 * find_voronoi_cell and compute_cell routines 169 * (but not needed in this instance.) 170 * \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)171 inline int region_index(int ci,int cj,int ck,int ei,int ej,int ek,double &qx,double &qy,double &qz,int &disp) { 172 int qi=ci+(ei-nx),qj=cj+(ej-ey),qk=ck+(ek-ez); 173 int iv(step_div(qi,nx));if(iv!=0) {qx=iv*bx;qi-=nx*iv;} else qx=0; 174 create_periodic_image(qi,qj,qk); 175 return qi+nx*(qj+oy*qk); 176 } 177 void create_all_images(); 178 void check_compartmentalized(); 179 protected: 180 void add_particle_memory(int i); 181 void put_locate_block(int &ijk,double &x,double &y,double &z); 182 void put_locate_block(int &ijk,double &x,double &y,double &z,int &ai,int &aj,int &ak); 183 /** Creates particles within an image block by copying them 184 * from the primary domain and shifting them. If the given 185 * block is aligned with the primary domain in the z-direction, 186 * the routine calls the simpler create_side_image routine 187 * where the image block may comprise of particles from up to 188 * two primary blocks. Otherwise is calls the more complex 189 * create_vertical_image where the image block may comprise of 190 * particles from up to four primary blocks. 191 * \param[in] (di,dj,dk) the coordinates of the image block to 192 * create. */ create_periodic_image(int di,int dj,int dk)193 inline void create_periodic_image(int di,int dj,int dk) { 194 if(di<0||di>=nx||dj<0||dj>=oy||dk<0||dk>=oz) 195 voro_fatal_error("Constructing periodic image for nonexistent point",VOROPP_INTERNAL_ERROR); 196 if(dk>=ez&&dk<wz) { 197 if(dj<ey||dj>=wy) create_side_image(di,dj,dk); 198 } else create_vertical_image(di,dj,dk); 199 } 200 void create_side_image(int di,int dj,int dk); 201 void create_vertical_image(int di,int dj,int dk); 202 void put_image(int reg,int fijk,int l,double dx,double dy,double dz); 203 inline void remap(int &ai,int &aj,int &ak,int &ci,int &cj,int &ck,double &x,double &y,double &z,int &ijk); 204 }; 205 206 /** \brief Extension of the container_periodic_base class for computing regular 207 * Voronoi tessellations. 208 * 209 * This class is an extension of the container_periodic_base that has routines 210 * specifically for computing the regular Voronoi tessellation with no 211 * dependence on particle radii. */ 212 class container_periodic : public container_periodic_base, public radius_mono { 213 public: 214 container_periodic(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_, 215 int nx_,int ny_,int nz_,int init_mem_); 216 void clear(); 217 void put(int n,double x,double y,double z); 218 void put(int n,double x,double y,double z,int &ai,int &aj,int &ak); 219 void put(particle_order &vo,int n,double x,double y,double z); 220 void import(FILE *fp=stdin); 221 void import(particle_order &vo,FILE *fp=stdin); 222 /** Imports a list of particles from an open file stream into 223 * the container. Entries of four numbers (Particle ID, x 224 * position, y position, z position) are searched for. If the 225 * file cannot be successfully read, then the routine causes a 226 * fatal error. 227 * \param[in] filename the name of the file to open and read 228 * from. */ import(const char * filename)229 inline void import(const char* filename) { 230 FILE *fp=safe_fopen(filename,"r"); 231 import(fp); 232 fclose(fp); 233 } 234 /** Imports a list of particles from an open file stream into 235 * the container. Entries of four numbers (Particle ID, x 236 * position, y position, z position) are searched for. In 237 * addition, the order in which particles are read is saved 238 * into an ordering class. If the file cannot be successfully 239 * read, then the routine causes a fatal error. 240 * \param[in,out] vo the ordering class to use. 241 * \param[in] filename the name of the file to open and read 242 * from. */ import(particle_order & vo,const char * filename)243 inline void import(particle_order &vo,const char* filename) { 244 FILE *fp=safe_fopen(filename,"r"); 245 import(vo,fp); 246 fclose(fp); 247 } 248 void compute_all_cells(); 249 double sum_cell_volumes(); 250 /** Dumps particle IDs and positions to a file. 251 * \param[in] vl the loop class to use. 252 * \param[in] fp a file handle to write to. */ 253 template<class c_loop> draw_particles(c_loop & vl,FILE * fp)254 void draw_particles(c_loop &vl,FILE *fp) { 255 double *pp; 256 if(vl.start()) do { 257 pp=p[vl.ijk]+3*vl.q; 258 fprintf(fp,"%d %g %g %g\n",id[vl.ijk][vl.q],*pp,pp[1],pp[2]); 259 } while(vl.inc()); 260 } 261 /** Dumps all of the particle IDs and positions to a file. 262 * \param[in] fp a file handle to write to. */ draw_particles(FILE * fp=stdout)263 inline void draw_particles(FILE *fp=stdout) { 264 c_loop_all_periodic vl(*this); 265 draw_particles(vl,fp); 266 } 267 /** Dumps all of the particle IDs and positions to a file. 268 * \param[in] filename the name of the file to write to. */ draw_particles(const char * filename)269 inline void draw_particles(const char *filename) { 270 FILE *fp=safe_fopen(filename,"w"); 271 draw_particles(fp); 272 fclose(fp); 273 } 274 /** Dumps particle positions in POV-Ray format. 275 * \param[in] vl the loop class to use. 276 * \param[in] fp a file handle to write to. */ 277 template<class c_loop> draw_particles_pov(c_loop & vl,FILE * fp)278 void draw_particles_pov(c_loop &vl,FILE *fp) { 279 double *pp; 280 if(vl.start()) do { 281 pp=p[vl.ijk]+3*vl.q; 282 fprintf(fp,"// id %d\nsphere{<%g,%g,%g>,s}\n", 283 id[vl.ijk][vl.q],*pp,pp[1],pp[2]); 284 } while(vl.inc()); 285 } 286 /** Dumps all particle positions in POV-Ray format. 287 * \param[in] fp a file handle to write to. */ draw_particles_pov(FILE * fp=stdout)288 inline void draw_particles_pov(FILE *fp=stdout) { 289 c_loop_all_periodic vl(*this); 290 draw_particles_pov(vl,fp); 291 } 292 /** Dumps all particle positions in POV-Ray format. 293 * \param[in] filename the name of the file to write to. */ draw_particles_pov(const char * filename)294 inline void draw_particles_pov(const char *filename) { 295 FILE *fp=safe_fopen(filename,"w"); 296 draw_particles_pov(fp); 297 fclose(fp); 298 } 299 /** Computes Voronoi cells and saves the output in gnuplot 300 * format. 301 * \param[in] vl the loop class to use. 302 * \param[in] fp a file handle to write to. */ 303 template<class c_loop> draw_cells_gnuplot(c_loop & vl,FILE * fp)304 void draw_cells_gnuplot(c_loop &vl,FILE *fp) { 305 voronoicell c;double *pp; 306 if(vl.start()) do if(compute_cell(c,vl)) { 307 pp=p[vl.ijk]+ps*vl.q; 308 c.draw_gnuplot(*pp,pp[1],pp[2],fp); 309 } while(vl.inc()); 310 } 311 /** Computes all Voronoi cells and saves the output in gnuplot 312 * format. 313 * \param[in] fp a file handle to write to. */ draw_cells_gnuplot(FILE * fp=stdout)314 inline void draw_cells_gnuplot(FILE *fp=stdout) { 315 c_loop_all_periodic vl(*this); 316 draw_cells_gnuplot(vl,fp); 317 } 318 /** Compute all Voronoi cells and saves the output in gnuplot 319 * format. 320 * \param[in] filename the name of the file to write to. */ draw_cells_gnuplot(const char * filename)321 inline void draw_cells_gnuplot(const char *filename) { 322 FILE *fp=safe_fopen(filename,"w"); 323 draw_cells_gnuplot(fp); 324 fclose(fp); 325 } 326 /** Computes Voronoi cells and saves the output in POV-Ray 327 * format. 328 * \param[in] vl the loop class to use. 329 * \param[in] fp a file handle to write to. */ 330 template<class c_loop> draw_cells_pov(c_loop & vl,FILE * fp)331 void draw_cells_pov(c_loop &vl,FILE *fp) { 332 voronoicell c;double *pp; 333 if(vl.start()) do if(compute_cell(c,vl)) { 334 fprintf(fp,"// cell %d\n",id[vl.ijk][vl.q]); 335 pp=p[vl.ijk]+ps*vl.q; 336 c.draw_pov(*pp,pp[1],pp[2],fp); 337 } while(vl.inc()); 338 } 339 /** Computes all Voronoi cells and saves the output in POV-Ray 340 * format. 341 * \param[in] fp a file handle to write to. */ draw_cells_pov(FILE * fp=stdout)342 inline void draw_cells_pov(FILE *fp=stdout) { 343 c_loop_all_periodic vl(*this); 344 draw_cells_pov(vl,fp); 345 } 346 /** Computes all Voronoi cells and saves the output in POV-Ray 347 * format. 348 * \param[in] filename the name of the file to write to. */ draw_cells_pov(const char * filename)349 inline void draw_cells_pov(const char *filename) { 350 FILE *fp=safe_fopen(filename,"w"); 351 draw_cells_pov(fp); 352 fclose(fp); 353 } 354 /** Computes the Voronoi cells and saves customized information 355 * about them. 356 * \param[in] vl the loop class to use. 357 * \param[in] format the custom output string to use. 358 * \param[in] fp a file handle to write to. */ 359 template<class c_loop> print_custom(c_loop & vl,const char * format,FILE * fp)360 void print_custom(c_loop &vl,const char *format,FILE *fp) { 361 int ijk,q;double *pp; 362 if(contains_neighbor(format)) { 363 voronoicell_neighbor c; 364 if(vl.start()) do if(compute_cell(c,vl)) { 365 ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q; 366 c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],default_radius,fp); 367 } while(vl.inc()); 368 } else { 369 voronoicell c; 370 if(vl.start()) do if(compute_cell(c,vl)) { 371 ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q; 372 c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],default_radius,fp); 373 } while(vl.inc()); 374 } 375 } 376 void print_custom(const char *format,FILE *fp=stdout); 377 void print_custom(const char *format,const char *filename); 378 bool find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid); 379 /** Computes the Voronoi cell for a particle currently being 380 * referenced by a loop class. 381 * \param[out] c a Voronoi cell class in which to store the 382 * computed cell. 383 * \param[in] vl the loop class to use. 384 * \return True if the cell was computed. If the cell cannot be 385 * computed because it was removed entirely for some reason, 386 * then the routine returns false. */ 387 template<class v_cell,class c_loop> compute_cell(v_cell & c,c_loop & vl)388 inline bool compute_cell(v_cell &c,c_loop &vl) { 389 return vc.compute_cell(c,vl.ijk,vl.q,vl.i,vl.j,vl.k); 390 } 391 /** Computes the Voronoi cell for given particle. 392 * \param[out] c a Voronoi cell class in which to store the 393 * computed cell. 394 * \param[in] ijk the block that the particle is within. 395 * \param[in] q the index of the particle within the block. 396 * \return True if the cell was computed. If the cell cannot be 397 * computed because it was removed entirely for some reason, 398 * then the routine returns false. */ 399 template<class v_cell> compute_cell(v_cell & c,int ijk,int q)400 inline bool compute_cell(v_cell &c,int ijk,int q) { 401 int k(ijk/(nx*oy)),ijkt(ijk-(nx*oy)*k),j(ijkt/nx),i(ijkt-j*nx); 402 return vc.compute_cell(c,ijk,q,i,j,k); 403 } 404 /** Computes the Voronoi cell for a ghost particle at a given 405 * location. 406 * \param[out] c a Voronoi cell class in which to store the 407 * computed cell. 408 * \param[in] (x,y,z) the location of the ghost particle. 409 * \return True if the cell was computed. If the cell cannot be 410 * computed, if it is removed entirely by a wall or boundary 411 * condition, then the routine returns false. */ 412 template<class v_cell> compute_ghost_cell(v_cell & c,double x,double y,double z)413 inline bool compute_ghost_cell(v_cell &c,double x,double y,double z) { 414 int ijk; 415 put_locate_block(ijk,x,y,z); 416 double *pp=p[ijk]+3*co[ijk]++; 417 *(pp++)=x;*(pp++)=y;*(pp++)=z; 418 bool q=compute_cell(c,ijk,co[ijk]-1); 419 co[ijk]--; 420 return q; 421 } 422 private: 423 voro_compute<container_periodic> vc; 424 friend class voro_compute<container_periodic>; 425 }; 426 427 /** \brief Extension of the container_periodic_base class for computing radical 428 * Voronoi tessellations. 429 * 430 * This class is an extension of container_periodic_base that has routines 431 * specifically for computing the radical Voronoi tessellation that depends 432 * on the particle radii. */ 433 class container_periodic_poly : public container_periodic_base, public radius_poly { 434 public: 435 container_periodic_poly(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_, 436 int nx_,int ny_,int nz_,int init_mem_); 437 void clear(); 438 void put(int n,double x,double y,double z,double r); 439 void put(int n,double x,double y,double z,double r,int &ai,int &aj,int &ak); 440 void put(particle_order &vo,int n,double x,double y,double z,double r); 441 void import(FILE *fp=stdin); 442 void import(particle_order &vo,FILE *fp=stdin); 443 /** Imports a list of particles from an open file stream into 444 * the container_poly class. Entries of five numbers (Particle 445 * ID, x position, y position, z position, radius) are searched 446 * for. If the file cannot be successfully read, then the 447 * routine causes a fatal error. 448 * \param[in] filename the name of the file to open and read 449 * from. */ import(const char * filename)450 inline void import(const char* filename) { 451 FILE *fp=safe_fopen(filename,"r"); 452 import(fp); 453 fclose(fp); 454 } 455 /** Imports a list of particles from an open file stream into 456 * the container_poly class. Entries of five numbers (Particle 457 * ID, x position, y position, z position, radius) are searched 458 * for. In addition, the order in which particles are read is 459 * saved into an ordering class. If the file cannot be 460 * successfully read, then the routine causes a fatal error. 461 * \param[in,out] vo the ordering class to use. 462 * \param[in] filename the name of the file to open and read 463 * from. */ import(particle_order & vo,const char * filename)464 inline void import(particle_order &vo,const char* filename) { 465 FILE *fp=safe_fopen(filename,"r"); 466 import(vo,fp); 467 fclose(fp); 468 } 469 void compute_all_cells(); 470 double sum_cell_volumes(); 471 /** Dumps particle IDs, positions and radii to a file. 472 * \param[in] vl the loop class to use. 473 * \param[in] fp a file handle to write to. */ 474 template<class c_loop> draw_particles(c_loop & vl,FILE * fp)475 void draw_particles(c_loop &vl,FILE *fp) { 476 double *pp; 477 if(vl.start()) do { 478 pp=p[vl.ijk]+4*vl.q; 479 fprintf(fp,"%d %g %g %g %g\n",id[vl.ijk][vl.q],*pp,pp[1],pp[2],pp[3]); 480 } while(vl.inc()); 481 } 482 /** Dumps all of the particle IDs, positions and radii to a 483 * file. 484 * \param[in] fp a file handle to write to. */ draw_particles(FILE * fp=stdout)485 inline void draw_particles(FILE *fp=stdout) { 486 c_loop_all_periodic vl(*this); 487 draw_particles(vl,fp); 488 } 489 /** Dumps all of the particle IDs, positions and radii to a 490 * file. 491 * \param[in] filename the name of the file to write to. */ draw_particles(const char * filename)492 inline void draw_particles(const char *filename) { 493 FILE *fp=safe_fopen(filename,"w"); 494 draw_particles(fp); 495 fclose(fp); 496 } 497 /** Dumps particle positions in POV-Ray format. 498 * \param[in] vl the loop class to use. 499 * \param[in] fp a file handle to write to. */ 500 template<class c_loop> draw_particles_pov(c_loop & vl,FILE * fp)501 void draw_particles_pov(c_loop &vl,FILE *fp) { 502 double *pp; 503 if(vl.start()) do { 504 pp=p[vl.ijk]+4*vl.q; 505 fprintf(fp,"// id %d\nsphere{<%g,%g,%g>,%g}\n", 506 id[vl.ijk][vl.q],*pp,pp[1],pp[2],pp[3]); 507 } while(vl.inc()); 508 } 509 /** Dumps all the particle positions in POV-Ray format. 510 * \param[in] fp a file handle to write to. */ draw_particles_pov(FILE * fp=stdout)511 inline void draw_particles_pov(FILE *fp=stdout) { 512 c_loop_all_periodic vl(*this); 513 draw_particles_pov(vl,fp); 514 } 515 /** Dumps all the particle positions in POV-Ray format. 516 * \param[in] filename the name of the file to write to. */ draw_particles_pov(const char * filename)517 inline void draw_particles_pov(const char *filename) { 518 FILE *fp(safe_fopen(filename,"w")); 519 draw_particles_pov(fp); 520 fclose(fp); 521 } 522 /** Computes Voronoi cells and saves the output in gnuplot 523 * format. 524 * \param[in] vl the loop class to use. 525 * \param[in] fp a file handle to write to. */ 526 template<class c_loop> draw_cells_gnuplot(c_loop & vl,FILE * fp)527 void draw_cells_gnuplot(c_loop &vl,FILE *fp) { 528 voronoicell c;double *pp; 529 if(vl.start()) do if(compute_cell(c,vl)) { 530 pp=p[vl.ijk]+ps*vl.q; 531 c.draw_gnuplot(*pp,pp[1],pp[2],fp); 532 } while(vl.inc()); 533 } 534 /** Compute all Voronoi cells and saves the output in gnuplot 535 * format. 536 * \param[in] fp a file handle to write to. */ draw_cells_gnuplot(FILE * fp=stdout)537 inline void draw_cells_gnuplot(FILE *fp=stdout) { 538 c_loop_all_periodic vl(*this); 539 draw_cells_gnuplot(vl,fp); 540 } 541 /** Compute all Voronoi cells and saves the output in gnuplot 542 * format. 543 * \param[in] filename the name of the file to write to. */ draw_cells_gnuplot(const char * filename)544 inline void draw_cells_gnuplot(const char *filename) { 545 FILE *fp(safe_fopen(filename,"w")); 546 draw_cells_gnuplot(fp); 547 fclose(fp); 548 } 549 /** Computes Voronoi cells and saves the output in POV-Ray 550 * format. 551 * \param[in] vl the loop class to use. 552 * \param[in] fp a file handle to write to. */ 553 template<class c_loop> draw_cells_pov(c_loop & vl,FILE * fp)554 void draw_cells_pov(c_loop &vl,FILE *fp) { 555 voronoicell c;double *pp; 556 if(vl.start()) do if(compute_cell(c,vl)) { 557 fprintf(fp,"// cell %d\n",id[vl.ijk][vl.q]); 558 pp=p[vl.ijk]+ps*vl.q; 559 c.draw_pov(*pp,pp[1],pp[2],fp); 560 } while(vl.inc()); 561 } 562 /** Computes all Voronoi cells and saves the output in POV-Ray 563 * format. 564 * \param[in] fp a file handle to write to. */ draw_cells_pov(FILE * fp=stdout)565 inline void draw_cells_pov(FILE *fp=stdout) { 566 c_loop_all_periodic vl(*this); 567 draw_cells_pov(vl,fp); 568 } 569 /** Computes all Voronoi cells and saves the output in POV-Ray 570 * format. 571 * \param[in] filename the name of the file to write to. */ draw_cells_pov(const char * filename)572 inline void draw_cells_pov(const char *filename) { 573 FILE *fp(safe_fopen(filename,"w")); 574 draw_cells_pov(fp); 575 fclose(fp); 576 } 577 /** Computes the Voronoi cells and saves customized information 578 * about them. 579 * \param[in] vl the loop class to use. 580 * \param[in] format the custom output string to use. 581 * \param[in] fp a file handle to write to. */ 582 template<class c_loop> print_custom(c_loop & vl,const char * format,FILE * fp)583 void print_custom(c_loop &vl,const char *format,FILE *fp) { 584 int ijk,q;double *pp; 585 if(contains_neighbor(format)) { 586 voronoicell_neighbor c; 587 if(vl.start()) do if(compute_cell(c,vl)) { 588 ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q; 589 c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],pp[3],fp); 590 } while(vl.inc()); 591 } else { 592 voronoicell c; 593 if(vl.start()) do if(compute_cell(c,vl)) { 594 ijk=vl.ijk;q=vl.q;pp=p[ijk]+ps*q; 595 c.output_custom(format,id[ijk][q],*pp,pp[1],pp[2],pp[3],fp); 596 } while(vl.inc()); 597 } 598 } 599 /** Computes the Voronoi cell for a particle currently being 600 * referenced by a loop class. 601 * \param[out] c a Voronoi cell class in which to store the 602 * computed cell. 603 * \param[in] vl the loop class to use. 604 * \return True if the cell was computed. If the cell cannot be 605 * computed because it was removed entirely for some reason, 606 * then the routine returns false. */ 607 template<class v_cell,class c_loop> compute_cell(v_cell & c,c_loop & vl)608 inline bool compute_cell(v_cell &c,c_loop &vl) { 609 return vc.compute_cell(c,vl.ijk,vl.q,vl.i,vl.j,vl.k); 610 } 611 /** Computes the Voronoi cell for given particle. 612 * \param[out] c a Voronoi cell class in which to store the 613 * computed cell. 614 * \param[in] ijk the block that the particle is within. 615 * \param[in] q the index of the particle within the block. 616 * \return True if the cell was computed. If the cell cannot be 617 * computed because it was removed entirely for some reason, 618 * then the routine returns false. */ 619 template<class v_cell> compute_cell(v_cell & c,int ijk,int q)620 inline bool compute_cell(v_cell &c,int ijk,int q) { 621 int k(ijk/(nx*oy)),ijkt(ijk-(nx*oy)*k),j(ijkt/nx),i(ijkt-j*nx); 622 return vc.compute_cell(c,ijk,q,i,j,k); 623 } 624 /** Computes the Voronoi cell for a ghost particle at a given 625 * location. 626 * \param[out] c a Voronoi cell class in which to store the 627 * computed cell. 628 * \param[in] (x,y,z) the location of the ghost particle. 629 * \param[in] r the radius of the ghost particle. 630 * \return True if the cell was computed. If the cell cannot be 631 * computed, if it is removed entirely by a wall or boundary 632 * condition, then the routine returns false. */ 633 template<class v_cell> compute_ghost_cell(v_cell & c,double x,double y,double z,double r)634 inline bool compute_ghost_cell(v_cell &c,double x,double y,double z,double r) { 635 int ijk; 636 put_locate_block(ijk,x,y,z); 637 double *pp=p[ijk]+4*co[ijk]++,tm=max_radius; 638 *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r; 639 if(r>max_radius) max_radius=r; 640 bool q=compute_cell(c,ijk,co[ijk]-1); 641 co[ijk]--;max_radius=tm; 642 return q; 643 } 644 void print_custom(const char *format,FILE *fp=stdout); 645 void print_custom(const char *format,const char *filename); 646 bool find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid); 647 private: 648 voro_compute<container_periodic_poly> vc; 649 friend class voro_compute<container_periodic_poly>; 650 }; 651 652 } 653 654 #endif 655