branches/exact/src/container.hh

   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:
  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. */
  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. */
  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. */
  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>
  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>
 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. */
 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                  */
 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. */
 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. */
 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. */
 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. */
 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. */
 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. */
 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>
 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. */
 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. */
 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>
 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. */
 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. */
 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>
 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. */
 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. */
 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>
 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. */
 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. */
 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>
 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>
 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>
 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. */
 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. */
 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>
 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. */
 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. */
 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>
 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. */
 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. */
 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>
 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. */
 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. */
 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>
 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. */
 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. */
 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>
 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>
 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>
 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